StRoot: StEmcCalibrationMaker/StEmcEqualMaker.cxx Source File

00001 #include "StEmcEqualMaker.h"
00002 #include "TFile.h"
00003 #include "TROOT.h"
00004 #include "TF1.h"
00005 #include "StEventTypes.h"
00006 #include "StEvent.h"
00007 #include "StEmcUtil/geometry/StEmcGeom.h"
00008 #include "StDbLib/StDbManager.hh"
00009 #include "StDbLib/StDbTable.h"
00010 #include "StDbLib/StDbConfigNode.hh"
00011 #include "tables/St_emcCalib_Table.h"
00012 
00013 ClassImp(StEmcEqualMaker)
00014 
00015 //_____________________________________________________________________________
00016 StEmcEqualMaker::StEmcEqualMaker(const char *name):StEmcCalibMaker(name)
00017 {
00018   setRange(700);
00019 }
00020 //_____________________________________________________________________________
00021 StEmcEqualMaker::~StEmcEqualMaker()
00022 {
00023 }
00024 //_____________________________________________________________________________
00025 Int_t StEmcEqualMaker::Init()
00026 {
00027   float pi = 3.1415926;
00028         mA = new TH1F("mA","",getNChannel(),0.5,getNChannel()+0.5);
00029   mB = new TH1F("mB","",getNChannel(),0.5,getNChannel()+0.5);
00030   mADistr = new TH1F("mADistr","",200,0,10);
00031   mBDistr = new TH1F("mBDistr","",200,-10,10);
00032   mStatus = new TH1F("mStatus","",getNChannel(),0.5,getNChannel()+0.5);
00033   mRefSlopes = new TH1F("mRefSlopes","",getNChannel(),0.5,getNChannel()+0.5);
00034   mRefAmps = new TH1F("mRefAmps","",getNChannel(),0.5,getNChannel()+0.5);
00035   mSlopes = new TH1F("mSlopes","",getNChannel(),0.5,getNChannel()+0.5);
00036   mAmps = new TH1F("mAmps","",getNChannel(),0.5,getNChannel()+0.5);
00037   mSlopesTheta = new TH2F("mSlopesTheta","",40,45*pi/180,135*pi/180,100,0,100);
00038 
00039   mSpecName ="mSpecEqual";
00040   mAcceptName = "mAcceptEqual";
00041   return StEmcCalibMaker::Init();
00042 }
00043 //_____________________________________________________________________________
00044 void StEmcEqualMaker::Clear(Option_t *option)              
00045 {
00046 }
00047 //_____________________________________________________________________________
00048 Int_t StEmcEqualMaker::Make()
00049 {  
00050   if(!accept()) return kStOk;
00051       
00052   for(int i=0;i<mNChannel;i++)
00053   {
00054     int id = i+1;
00055     float rms = getCalib()->getPedRms(mDetector,id);
00056     float adc = (float) getCalib()->getADCPedSub(mDetector,id);
00057     if(adc!=0 && adc>3*rms) fill(id,adc);
00058   }
00059   
00060   return kStOK;
00061 }
00062 //_____________________________________________________________________________
00063 Int_t StEmcEqualMaker::Finish()
00064 {
00065   saveHist((char*)mFileName.Data());
00066   return StMaker::Finish();
00067 }
00068 //_____________________________________________________________________________
00069 void StEmcEqualMaker::equalize(int mode,int DeltaEta,bool Et)
00070 {
00071         mADistr->Reset();
00072         mBDistr->Reset();
00073         StEmcGeom *g = getGeom();
00074         int Neta = g->NEta();
00075         int Nsub = g->NSub();
00076         int etabin = 0;
00077         float avg[18000];
00078         float MIN = 20;
00079         float MAX = 500;
00080         TF1 *f= new TF1("slopes","22*TMath::Sin(x)");
00081 
00082         for(int eta0 = -Neta; eta0<Neta; eta0+=DeltaEta)
00083         {
00084                 etabin++;
00085                 float sum = 0;
00086                 float n = 0;
00087                 int m1 = 1;
00088                 int m2 = 60;
00089                 if(eta0<0) { m1 = 61; m2 = 120; }
00090                 if(DeltaEta == 2*Neta){ m1=1; m2=120;}
00091                 for(int j = 0;j<DeltaEta;j++)
00092                 {
00093                         int eta = eta0+j;
00094                         if(eta>=0) eta+=1;
00095                         cout <<"etabin = "<<etabin<<" eta = "<<eta<<endl;
00096                         eta = abs(eta);
00097                         //cout <<"Calculating averages ..."<<endl;
00098                         for(int m = m1;m<=m2;m++)
00099                         {
00100                                 for(int sub = 1;sub<=Nsub;sub++)
00101                                 {
00102                                   int id;
00103                                         g->getId(m,eta,sub,id);
00104                                         TH1D * h = getSpec(id);
00105                                         if(h->Integral()>0)
00106                                         {
00107                                                 float m1,r1;
00108                                                 getMeanAndRms(h,MIN,MAX,&m1,&r1);
00109                                                 sum+=m1;
00110                                                 n++;
00111                                                 avg[id-1] = m1;
00112                                         }
00113                                 }
00114                         }
00115                 }
00116                 float AVG = 0;
00117                 if(n>0) AVG=sum/n;
00118                 if(AVG>0)
00119                 {
00120                         int ID = 0;
00121                         float DCA = 99999;
00122                         //cout <<"Looking for reference ..."<<endl;
00123                         for(int m = m1;m<=m2;m++) for(int j = 0;j<DeltaEta;j++) for(int sub = 1;sub<=Nsub;sub++)
00124                         {
00125                                 int id;
00126                                 int eta = eta0+j;
00127                                 if(eta>=0) eta+=1;
00128                           eta = abs(eta);
00129                                 g->getId(m,eta,sub,id);
00130                                 if(avg[id-1]>0 && fabs(avg[id-1]-AVG)<DCA) { DCA = fabs(avg[id-1]-AVG); ID = id;}
00131                         }
00132                         if(ID>0)
00133                         {
00134                                 cout <<"Reference found id = "<<ID<<"  avg = "<<avg[ID-1]<<endl;
00135                           cout <<"Equalizing eta bin ..."<<endl;
00136                                 for(int m = m1;m<=m2;m++) for(int j = 0;j<DeltaEta;j++) for(int sub = 1;sub<=Nsub;sub++) 
00137                                 {
00138                                         int eta = eta0+j;
00139                                         if(eta>=0) eta+=1;
00140                             eta = abs(eta);
00141                                         int id;
00142                                         g->getId(m,eta,sub,id);
00143                                         if(mode!=5) equalizeRelative(ID,id,mode,Et);
00144                                         else equalizeToFunction(id,f);
00145                                         
00146                                 }
00147                         } 
00148                 }
00149         }
00150         delete f;
00151         return;
00152 }
00153 //_____________________________________________________________________________
00154 void StEmcEqualMaker::equalizeRelative(int ref, int channel,int mode, bool Et)
00155 {
00156   // mean and RMS modes  
00157   // 0 -  mean and RMS with liear average
00158   // 1 -  mean with linear average
00159   // 2 -  mean and rms with log average
00160   // 3 -  mean with log average 
00161         // 4 -  exponential fit
00162 
00163   TH1D *h1 = getSpec(ref,"ref");
00164         TH1D *h2 = getSpec(channel,"channel");
00165 
00166   float integral1 = h1->Integral();
00167         float integral2 = h2->Integral();
00168         
00169         mA->SetBinContent(channel,0);
00170         mB->SetBinContent(channel,0);
00171         mStatus->SetBinContent(channel,0);
00172         mRefSlopes->SetBinContent(channel,0);
00173         mRefAmps->SetBinContent(channel,0);
00174         mSlopes->SetBinContent(channel,0);
00175         mAmps->SetBinContent(channel,0);
00176         
00177         float MIN = 20;
00178         float MAX = 500;
00179         
00180         //cout <<"ref = "<<ref<<"  channel = "<<channel<<endl;
00181         //cout <<"h1 = "<<h1<<"  h2 = "<<h2<<endl;
00182         //cout <<"integral 1 = "<<integral1<<"  integral 2 = "<<integral2<<endl;
00183         
00184         if((integral1==0 || integral2==0) && h1!=h2)
00185         {
00186                 delete h1;
00187                 delete h2;
00188                 return;
00189         }   
00190         if((integral1==0 || integral2==0) && h1==h2)
00191         {
00192                 delete h1;
00193                 return;
00194         }   
00195         
00196   bool EqDone=kFALSE;
00197   float a=0,b=0;
00198   
00199   if(mode==-1)
00200   {
00201     a = 1;
00202     b = 0;
00203     EqDone=true;
00204   }
00205   
00206   if(mode>=0 && mode <=3)  
00207   {
00208     float m1,r1,m2,r2;        
00209     if(mode==0 || mode==1)
00210     {
00211       getMeanAndRms(h1,MIN,MAX,&m1,&r1);
00212       getMeanAndRms(h2,MIN,MAX,&m2,&r2);
00213     }
00214     if(mode==2 || mode==3)
00215     {
00216       getLogMeanAndRms(h1,MIN,MAX,&m1,&r1);
00217       getLogMeanAndRms(h2,MIN,MAX,&m2,&r2);
00218     }
00219     if(mode==0 || mode==2)
00220     {
00221       a=r1/r2;
00222       b=m1-a*m2;
00223     }
00224     if(mode==1 || mode==3)
00225     {    
00226       a=m1/m2;
00227       b=0;
00228     }
00229     EqDone=kTRUE;
00230     cout <<"  id = "<<channel<<"  ref = "<<ref<<"  mean = "<<m1<<" , "<<m2
00231          <<"  rms = "<<r1<<" , "<<r2
00232          <<"  a = "<<a<<"  b = "<<b<<endl;
00233   }
00234   
00235         float m1=0,m2=0,A1=0,A2=0;
00236         if(mode==4)
00237         {
00238                 TF1 *f=new TF1("ff","[0]*exp(-x/[1])",MIN,MAX);
00239     float I1 = h1->Integral(h1->FindBin(MIN),h1->FindBin(MAX));
00240                 f->SetParameter(1,10);
00241                 h1->Fit(f,"RQN");
00242                 m1 = f->GetParameter(1);
00243                 A1 = f->GetParameter(0);
00244           mRefSlopes->SetBinContent(channel,m1);
00245           mRefAmps->SetBinContent(channel,A1);
00246         
00247     float I2 = h2->Integral(h2->FindBin(MIN),h2->FindBin(MAX));
00248                 h2->Fit(f,"RQN");
00249                 m2 = f->GetParameter(1);
00250                 A2 = f->GetParameter(0);
00251           mSlopes->SetBinContent(channel,m2);
00252           mAmps->SetBinContent(channel,A2);
00253                 
00254                 a=m1/m2;
00255                 b=-::log((A2*I1)/(A1*I2));
00256                 EqDone=true;
00257     if(!finite(a) || !finite(b) || a<=0 || b>1000) {EqDone = false; a = 0;}
00258     b=0;
00259     cout <<"  id = "<<channel<<"  ref = "<<ref<<"  slopes = "<<m2<<" , "<<m1
00260          <<"  a = "<<a<<"  b = "<<b<<"  EQDONE = "<<(Int_t)EqDone<<endl;
00261                 delete f;
00262         }
00263   
00264   if (EqDone)
00265   {
00266           mA->SetBinContent(channel,a);
00267           mB->SetBinContent(channel,b);
00268           mADistr->Fill(a);
00269           mBDistr->Fill(b);
00270           mStatus->SetBinContent(channel,1);
00271                 
00272                 float theta;
00273                 int mod,eta,sub;
00274                 getGeom()->getBin(channel,mod,eta,sub);
00275                 getGeom()->getTheta(mod,eta,theta);
00276                 mSlopesTheta->Fill(theta,m2);
00277           
00278   }  
00279   else mStatus->SetBinContent(channel,2);
00280         
00281         delete h1;
00282         delete h2;
00283   return;
00284 }
00285 //_____________________________________________________________________________
00286 void StEmcEqualMaker::equalizeToFunction(int channel,TF1 *func)
00287 {
00288         TH1D *h2 = getSpec(channel,"channel");
00289 
00290         float integral2 = h2->Integral();
00291         
00292         mA->SetBinContent(channel,0);
00293         mB->SetBinContent(channel,0);
00294         mStatus->SetBinContent(channel,0);
00295         mRefSlopes->SetBinContent(channel,0);
00296         mRefAmps->SetBinContent(channel,0);
00297         mSlopes->SetBinContent(channel,0);
00298         mAmps->SetBinContent(channel,0);
00299         
00300         float MIN = 20;
00301         float MAX = 500;
00302                 
00303         if(integral2==0) 
00304         { 
00305           delete h2;      
00306           mA->SetBinContent(channel,0);
00307           mB->SetBinContent(channel,0);
00308                 return;
00309   }   
00310         
00311   bool EqDone=kFALSE;
00312   float a=0,b=0;
00313     
00314         float theta;
00315         int mod,eta,sub;
00316         getGeom()->getBin(channel,mod,eta,sub);
00317         getGeom()->getTheta(mod,eta,theta);
00318 
00319         float m1=0,m2=0,A1=0,A2=0;
00320         TF1 *f=new TF1("ff","[0]*exp(-x/[1])",MIN,MAX);
00321         f->SetParameter(1,10);
00322         
00323         m1 = func->Eval(theta);
00324         mRefSlopes->SetBinContent(channel,m1);
00325         mRefAmps->SetBinContent(channel,A1);
00326         
00327   if(channel<=2400)
00328         {
00329           //float I2 = h2->Integral(h2->FindBin(MIN),h2->FindBin(MAX));
00330                 h2->Fit(f,"RQN");
00331                 m2 = f->GetParameter(1);
00332                 A2 = f->GetParameter(0);
00333                 mSlopes->SetBinContent(channel,m2);
00334                 mAmps->SetBinContent(channel,A2);
00335                 
00336                 
00337                 a=m1/m2;
00338         } //else a =1; //protection because crates in the east side are not timed
00339         
00340         EqDone=true;
00341   if(!finite(a) || !finite(b) || a<=0.01 || b>1000 || a > 10) { EqDone = false; a = 0;}
00342   b=0;
00343   cout <<"  id = "<<channel<<"  slopes = "<<m2<<" , "<<m1
00344          <<"  a = "<<a<<"  b = "<<b<<"  EQDONE = "<<(Int_t)EqDone<<endl;
00345         delete f;
00346   
00347   if (EqDone)
00348   {
00349           mA->SetBinContent(channel,a);
00350           mB->SetBinContent(channel,b);
00351           mADistr->Fill(a);
00352           mBDistr->Fill(b);
00353           mStatus->SetBinContent(channel,1);
00354                 
00355                 mSlopesTheta->Fill(theta,m2);
00356           
00357   }  
00358   else mStatus->SetBinContent(channel,2);
00359         
00360         delete h2;
00361   return;
00362 }
00363 //_____________________________________________________________________________
00364 void StEmcEqualMaker::calcSlopes()
00365 {
00366 
00367         float MIN = 20;
00368         float MAX = 500;
00369   mSlopes->Reset();
00370         mAmps->Reset();
00371         mSlopesTheta->Reset();  
00372         TF1 *f=new TF1("ff","[0]*exp(-x/[1])",MIN,MAX);
00373         StEmcGeom *g = getGeom();
00374 
00375         for(int channel = 1;channel<=getNChannel();channel++)
00376         {
00377                 TH1D *h2 = getSpec(channel,"channel");
00378         float I2 = h2->Integral(h2->FindBin(MIN),h2->FindBin(MAX));
00379                 if(I2>0)
00380                 {
00381       float theta;
00382                   int mod,eta,sub;
00383                   g->getBin(channel,mod,eta,sub);
00384                   g->getTheta(mod,eta,theta);
00385 
00386                   float m2=0,A2=0;
00387                 
00388             f->SetParameter(1,10);
00389                   h2->Fit(f,"RQN");
00390                   m2 = f->GetParameter(1);
00391                   A2 = f->GetParameter(0);
00392                   mSlopes->SetBinContent(channel,m2);
00393                   mAmps->SetBinContent(channel,A2);
00394                                 
00395                   mSlopesTheta->Fill(theta,m2);
00396                         cout <<"Slope for channel "<<channel<<"  is "<<m2<<endl;
00397                 }
00398           delete h2;
00399   }  
00400   delete f;
00401  
00402   return;
00403 }
00404 //_____________________________________________________________________________
00405 void StEmcEqualMaker::saveEqual(int date,int time)
00406 {   
00407         char ts[100];
00408   TString n[] = {"bemcEqual","bprsEqual","bsmdeEqual","bsmdpEqual"};
00409         sprintf(ts,"%s.%08d.%06d.root",n[getDetector()-1].Data(),date,time);
00410         TFile *f = new TFile(ts,"RECREATE");
00411         if(getSpec()) getSpec()->Write();
00412         if(mA) mA->Write();
00413         if(mB) mB->Write();
00414         if(mADistr) mADistr->Write();
00415         if(mBDistr) mBDistr->Write();
00416         if(mStatus) mStatus->Write();
00417         if(mRefSlopes) mRefSlopes->Write();
00418         if(mRefAmps) mRefAmps->Write();
00419         if(mSlopes) mSlopes->Write();
00420         if(mAmps) mAmps->Write();
00421         if(mSlopesTheta) mSlopesTheta->Write();
00422         f->Close();
00423         delete f;
00424         return;
00425 }
00426 //_____________________________________________________________________________
00427 void StEmcEqualMaker::loadEqual(char* file)
00428 {   
00429         TFile *f = new TFile(file);
00430         mA->Reset();
00431         mB->Reset();
00432         mADistr->Reset();
00433         mBDistr->Reset();
00434         mStatus->Reset();
00435         mRefSlopes->Reset();
00436         mRefAmps->Reset();
00437         mSlopes->Reset();
00438         mAmps->Reset();
00439         mSlopesTheta->Reset();
00440         getSpec()->Reset();
00441         
00442   if(f->Get("mSpec;1")) getSpec()->Add((TH2F*)f->Get("mSpec;1"),1);
00443   if(f->Get("mA;1")) mA->Add((TH1F*)f->Get("mA;1"),1);
00444   if(f->Get("mB;1")) mB->Add((TH1F*)f->Get("mB;1"),1);
00445   if(f->Get("mADistr;1")) mADistr->Add((TH1F*)f->Get("mADistr;1"),1);
00446   if(f->Get("mBDistr;1")) mBDistr->Add((TH1F*)f->Get("mBDistr;1"),1);
00447   if(f->Get("mStatus;1")) mStatus->Add((TH1F*)f->Get("mStatus;1"),1);
00448   if(f->Get("mRefSlopes;1")) mRefSlopes->Add((TH1F*)f->Get("mRefSlopes;1"),1);
00449   if(f->Get("mRefAmps;1")) mRefAmps->Add((TH1F*)f->Get("mRefAmps;1"),1);
00450   if(f->Get("mSlopes;1")) mSlopes->Add((TH1F*)f->Get("mSlopes;1"),1);
00451   if(f->Get("mAmps;1")) mAmps->Add((TH1F*)f->Get("mAmps;1"),1);
00452   if(f->Get("mSlopesTheta;1")) mSlopesTheta->Add((TH2F*)f->Get("mSlopesTheta;1"),1);
00453         f->Close();
00454         delete f;
00455         return;
00456 }
00457 //_____________________________________________________________________________
00458 TH1F* StEmcEqualMaker::getEtaBinSpec(int eta0, int eta1,TH2F* SPEC)
00459 {   
00460         TH2F *spec = NULL;
00461         if(!SPEC) spec = getSpec(); else spec = SPEC;
00462         if(!spec) return NULL;
00463         int mi = 1;
00464         int mf = 60;
00465         if(eta0<0 && eta1<0)
00466         {
00467                 mi = 61;
00468                 mf = 120;
00469         }
00470         if(eta0*eta1<0)
00471         {
00472                 mi =1;
00473                 mf = 120;
00474         }
00475         int ns = getGeom()->NSub();
00476         int e0 = abs(eta0);
00477         int e1 = abs(eta1);
00478         TH1F *hist = (TH1F*)spec->ProjectionY("etabin",1,1);
00479         hist->Reset();
00480         
00481         for(int m=mi;m<=mf;m++)
00482         {
00483                 for(int e=e0;e<=e1;e++)
00484                 {
00485                         for(int s=1;s<=ns;s++)
00486                         {
00487                                 int id;
00488                                 getGeom()->getId(m,e,s,id);
00489                                 char name[30];
00490                                 sprintf(name,"Tower %d",id);
00491                                 TH1F* tmp = rebin(id,name,spec);
00492                                 if(tmp)
00493                                 {
00494                                   hist->Add(tmp,1);
00495                                   delete tmp; tmp = NULL;
00496                                 }
00497                         }
00498                 }
00499         }
00500         return hist;
00501 }
00502 //_____________________________________________________________________________
00503 TH1F* StEmcEqualMaker::rebin(int id,char *name, TH2F* SPEC)
00504 {   
00505         TH2F *spec = NULL;
00506         if(!SPEC) spec = getSpec(); else spec = SPEC;
00507         if(!spec) return NULL;
00508         
00509         float seg = 50;
00510         float a = mA->GetBinContent(id);
00511         float b = mB->GetBinContent(id);
00512         //cout <<"A = "<<a<<"  B = "<<b<<endl;
00513         
00514   if(a<=0) return NULL; 
00515                 
00516         TH1F *h = (TH1F*)spec->ProjectionY(name,id,id);
00517         if(a==1 && b==0) return h;
00518         
00519         float deltaBin=h->GetBinWidth(1)/seg;
00520         int nbins=h->GetNbinsX();       
00521         
00522         float Y[4096];
00523         for(int i=0;i<4096;i++) Y[i] = 0;
00524         for(int i=1;i<=nbins;i++) Y[i] = h->GetBinContent(i);
00525                         
00526         h->Reset();
00527   
00528   for(int i=1;i<=nbins;i++)
00529   {
00530     float x = h->GetBinLowEdge(i);
00531                 float y = Y[i];
00532                 for(int j=0;j<(int)seg;j++)
00533                 {
00534                         float x1 = x+((float)j+0.5)*deltaBin;
00535                         float x2 = x1*a+b;
00536       //if (a>1) cout <<position<<"  "<<a<<"  "<<x1<<"  "<<x2<<endl;
00537                         h->Fill(x2,y/seg);
00538                 }
00539   }
00540         return h;
00541 }
00542 //_____________________________________________________________________________
00543 // this method leaks memory because it draws a lot of histograms
00544 // use it with care
00545 // AASPSUAIDE
00546 void StEmcEqualMaker::drawEtaBin(int eta0, int eta1,TH2F* SPEC)
00547 {   
00548         TH2F *spec = NULL;
00549         if(!SPEC) spec = getSpec(); else spec = SPEC;
00550         if(!spec) return;
00551         int mi = 1;
00552         int mf = 60;
00553         if(eta0<0 && eta1<0)
00554         {
00555                 mi = 61;
00556                 mf = 120;
00557         }
00558         if(eta0*eta1<0)
00559         {
00560                 mi =1;
00561                 mf = 120;
00562         }
00563         int ns = getGeom()->NSub();
00564         int e0 = abs(eta0);
00565         int e1 = abs(eta1);
00566 
00567   bool first = true;
00568                 
00569         for(int m=mi;m<=mf;m++)
00570         {
00571                 for(int e=e0;e<=e1;e++)
00572                 {
00573                         for(int s=1;s<=ns;s++)
00574                         {
00575                                 int id;
00576                                 getGeom()->getId(m,e,s,id);
00577                                 char name[30];
00578                                 sprintf(name,"Tower %d",id);
00579                                 TH1F* tmp = rebin(id,name,spec);
00580                                 if(tmp)
00581                                 {
00582                                   if(first) { tmp->Draw(); first = false;} else tmp->Draw("same");
00583                                 }
00584                         }
00585                 }
00586         }
00587         return;
00588 }
00589 
00590 
00591 
00592 
00593