StRoot: StEEmcPool/StEEmcTimingMaker/macros/plotEEmcTiming.C Source File

00001 TFile  *file  = 0;
00002 TChain *chain = 0;
00003 
00004 #include <vector>
00005 
00006 #include "/afs/rhic.bnl.gov/star/packages/DEV/StRoot/StEEmcUtil/EEfeeRaw/EEdims.h"
00007 
00008 // summary TTree branches created by StEEmcTimingMaker
00009 Int_t   mRunNumber;
00010 Float_t mTowerDelay;
00011 Float_t mMapmtDelay;
00012 Int_t   mTotalYield;
00013 Int_t   mTowerCrateYield[MaxTwCrates];
00014 Int_t   mMapmtCrateYield[MaxMapmtCrates];
00015 Int_t   mTowerChanYield[MaxTwCrates][MaxTwCrateCh];
00016 Int_t   mMapmtChanYield[MaxMapmtCrates][MaxMapmtCrateCh];
00017 Int_t   mTowerMin;
00018 Int_t   mTowerMax;
00019 Int_t   mMapmtMin;
00020 Int_t   mMapmtMax;
00021 
00022 // vectors to hold TGraphs for tower and mapmt crates
00023 Int_t npoints = 0;
00024 
00025 TGraphErrors *towerCrateCurves[MaxTwCrates];
00026 TGraphErrors *mapmtCrateCurves[MaxMapmtCrates];
00027 
00028 TGraphErrors *towerChanCurves[MaxTwCrates][MaxTwCrateCh];
00029 TGraphErrors *mapmtChanCurves[MaxMapmtCrates][MaxMapmtCrateCh];
00030 
00031 // enables printing of output files (gif) for documentation
00032 // figures will appear in same subdirector as input files, 
00033 // specified below.
00034 Bool_t doprint = true;
00035 //Bool_t doprint = false;
00036 
00037 void plotEEmcTiming( const Char_t *input_dir="timing_files/")
00038 {  
00039   
00040   // chain files
00041   chainFiles(input_dir);
00042 
00043   // setup branches
00044   setBranches(input_dir);
00045 
00046   // get total number of points
00047   Long64_t nruns = chain -> GetEntries();
00048   npoints=(Int_t)nruns;
00049 
00050   // setup the graphs for each crate and each channel
00051   initGraphs();
00052 
00053   // loop over all runs
00054   for ( Long64_t i = 0; i < nruns; i++ ) 
00055     {
00056       chain->GetEntry(i);
00057             
00058       fillCrates((int)i);
00059       fillChannels((int)i);
00060 
00061     }
00062 
00063   // draw timing scan curves for tower crates
00064   drawCrates();
00065   for ( Int_t ii=0;ii<MaxTwCrates;ii++ ) towerChannels(ii);
00066   for ( Int_t ii=0;ii<MaxMapmtCrates;ii++ ) mapmtChannels(ii);
00067 
00068   std::cout << "--------------------------------------------------------" << std::endl;
00069   std::cout << "to view timing curves for any crate" << std::endl;
00070   std::cout << std::endl;
00071   std::cout << "towerChannels(icrate) -- icrate = 0-5 for tower crates 1-6"<<std::endl;
00072   std::cout << "mapmtChannels(icrate) -- icrate = 0-47 for mapmt crates 1-48"<<std::endl;
00073 //  std::cout << "print()               -- make gif files for all crates"<<std::endl;
00074 
00075 }
00076 // ----------------------------------------------------------------------------
00077 void print()
00078 {
00079   drawCrates();
00080   for ( Int_t ii=0;ii<MaxTwCrates;ii++ ) drawChannels(ii);
00081   for ( Int_t ii=0;ii<MaxMapmtCrates;ii++ ) drawMapmt(ii);
00082 }
00083 
00084 // ----------------------------------------------------------------------------
00085 void drawCrates()
00086 {
00087 
00088   // tower crates first
00089   TCanvas *towers=new TCanvas("towers","towers",500,400);
00090   const Char_t *opt[]={"ALP","LP","LP","LP","LP","LP"};
00091 
00092   TLegend *legend=new TLegend(0.125,0.6,0.325,0.85);
00093 
00094   Float_t ymax=0.;
00095   for ( Int_t icr=0;icr<MaxTwCrates;icr++ )
00096     {
00097       towerCrateCurves[icr]->Sort();
00098       towerCrateCurves[icr]->Draw(opt[icr]);
00099       TString crname="tw crate ";crname+=icr+1;
00100       legend->AddEntry( towerCrateCurves[icr], crname, "lp" );
00101       if ( towerCrateCurves[icr]->GetYaxis()->GetXmax() > ymax )
00102         ymax=towerCrateCurves[icr]->GetYaxis()->GetXmax();
00103     }
00104   towerCrateCurves[0]->SetTitle("EEMC Tower Crate Timing Curves");
00105   towerCrateCurves[0]->GetXaxis()->SetTitle("TCD phase[ns]");
00106   TString ytitle=Form("Integral [ped+%i,ped+%i] / N_{events}",mTowerMin,mTowerMax);
00107   towerCrateCurves[0]->GetYaxis()->SetTitle(ytitle);
00108   towerCrateCurves[0]->GetYaxis()->SetRangeUser(0.,ymax);
00109   legend->Draw();
00110 
00111   if(doprint)towers->Print("tower_crates.gif");
00112 
00113   // next mapmt crates, 4 crates per canvas
00114   TCanvas *mapmt = 0;
00115 
00116   for ( Int_t icr = 0; icr<MaxMapmtCrates; icr+=4 ) 
00117     {
00118       
00119       TString cname="mapmt";cname+=icr+MinMapmtCrateID;
00120       mapmt=new TCanvas(cname,cname,500,400);
00121       TString title="EEMC Mapmt Crate Timing Curves ";
00122       title+=icr+MinMapmtCrateID;
00123       title+="-";
00124       title+=icr+MinMapmtCrateID+3;
00125       legend=new TLegend(0.625,0.12,0.825,0.375);
00126 
00127       Float_t ymax = 0.0;
00128       for ( Int_t jcr=0;jcr<4;jcr++ ) {
00129         Int_t index=icr+jcr;
00130         Int_t crateid=MinMapmtCrateID+index;
00131 
00132         mapmtCrateCurves[index]->Sort();
00133         mapmtCrateCurves[index]->Draw(opt[jcr]);
00134         if ( mapmtCrateCurves[index]->GetYaxis()->GetXmax() > ymax )
00135           ymax=mapmtCrateCurves[index]->GetYaxis()->GetXmax();
00136         TString crname="mapmt crate ";crname+=crateid;
00137         legend->AddEntry(mapmtCrateCurves[index],crname,"lp");
00138 
00139       }
00140       mapmtCrateCurves[icr]->SetTitle(title);
00141       mapmtCrateCurves[icr]->GetXaxis()->SetTitle("TCD phase[ns]");
00142       mapmtCrateCurves[icr]->GetYaxis()->SetTitle("Integral [ped+25,ped+125] / Nevents");
00143       TString ytitle=Form("Integral [ped+%i,ped+%i] / N_{events}",mMapmtMin,mMapmtMax);
00144       mapmtCrateCurves[icr]->GetYaxis()->SetTitle(ytitle);
00145       mapmtCrateCurves[icr]->GetYaxis()->SetRangeUser(0.,ymax*1.05);
00146       legend->Draw();
00147 
00148       TString oname="mapmt_crates_";
00149       oname+=icr+MinMapmtCrateID;
00150       oname+="_";
00151       oname+=icr+MinMapmtCrateID+4;
00152       if(doprint)mapmt->Print(oname+".gif");
00153 
00154     }
00155 
00156 
00157 }
00158 // ----------------------------------------------------------------------------
00159 void mapmtChannels( Int_t crate )
00160 {
00161   
00162   static const Int_t stride=16;
00163   Int_t crateid = MinMapmtCrateID+crate;
00164 
00165 
00166   TString fname="mapmt-crate-";fname+=crate+MinMapmtCrateID;fname+=".ps";
00167   TCanvas *canvas=new TCanvas("canvas","canvas",850/2,1100/2);
00168   canvas->Divide(1,2);
00169   Int_t icanvas=0;
00170 
00171 
00172   for ( Int_t ich=0;ich<192;ich+=stride ) 
00173     {
00174 
00175       canvas->cd(1+icanvas%2);
00176       icanvas++;
00177 
00178 
00179       TString pname="crate";
00180       pname+=crateid;
00181       pname+="_ch";
00182       pname+=ich;
00183       pname+="-";
00184       pname+=ich+stride-1;
00185 
00186       const Char_t *opts[]={"ALP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP"};
00187 
00188       // normalize
00189       Float_t ymax=0.0;
00190       Double_t sum[stride];for ( Int_t jj=0;jj<stride;jj++ )sum[jj]=0.;
00191       Double_t max=0.;
00192       for ( Int_t jch=0;jch<stride;jch++ ) // loop over channels in this one graph
00193         {
00194           Int_t index=ich+jch;
00195           Double_t *Y=mapmtChanCurves[crate][index]->GetY();
00196           for ( Int_t ipoint=0;ipoint<npoints;ipoint++ ) {
00197             if ( Y[ipoint]>ymax ) ymax=Y[ipoint];
00198             sum[jch]+=Y[ipoint];
00199           }
00200           if ( sum[jch]>max ) max=sum[jch];
00201         }
00202       if ( max <= 0. ) continue; // meh?
00203 
00204       
00205 
00206       TLegend *legend=new TLegend(0.55,0.11,0.85,0.525);
00207       for ( Int_t jch=0;jch<stride;jch++ )
00208         {
00209           Int_t index=ich+jch;
00210           mapmtChanCurves[crate][index]->SetMarkerSize(0.75);
00211           // offset X axis of each of these
00212           Double_t *X=mapmtChanCurves[crate][index]->GetX();
00213           Double_t *Y=mapmtChanCurves[crate][index]->GetY();
00214           Double_t *EY=mapmtChanCurves[crate][index]->GetEY();
00215           if ( sum[jch]<= 0. ) continue;
00216           //      std::cout<<"before"<<std::endl;
00217           for ( Int_t ip=0;ip<npoints;ip++ ){
00218             Float_t shift = 0.5+ ((float)jch) - ((float)stride)/2.0;
00219             Double_t yy=Y[ip];
00220             X[ip]-= 0.1*shift;
00221             Y[ip]*=max/sum[jch];
00222             EY[ip]*=max/sum[jch];
00223             //      std::cout << "ip="<<ip<<" y="<<yy<<" y'="<<Y[ip]<<std::endl;
00224           }
00225           mapmtChanCurves[crate][index]->Sort();
00226           if ( !jch ) 
00227             mapmtChanCurves[crate][index]->GetXaxis()->SetRangeUser(0.,ymax*1.05);
00228           mapmtChanCurves[crate][index]->SetMarkerColor(38+jch);
00229           mapmtChanCurves[crate][index]->SetLineColor(38+jch);    
00230           mapmtChanCurves[crate][index]->SetMinimum(0.);          
00231           mapmtChanCurves[crate][index]->Draw(opts[jch]);
00232           
00233           TString label="crate ";label+=crate+1;label+=" chan ";label+=index;
00234           legend->AddEntry(mapmtChanCurves[crate][index],label,"lp");
00235 
00236         }
00237       legend->Draw();
00238       canvas->Update();
00239 
00240       //      if(doprint)c->Print(pname+".gif");
00241       if ( !(icanvas%2) ){
00242         canvas->Print(fname+"(");
00243         canvas->Clear();
00244         canvas->Divide(1,2);
00245       }
00246       //      if(doprint)c->Print(pname+".gif");
00247 
00248     }
00249   canvas->Print(fname+")");
00250   gSystem->Exec(TString("ps2pdf ")+fname);
00251 
00252 }
00253 // ----------------------------------------------------------------------------
00254 void towerChannels( Int_t crate )
00255 {
00256 
00257   static  const  Int_t stride=12;
00258 
00259   TString fname="tower-crate-";fname+=crate+1;fname+=".ps";
00260   TCanvas *canvas=new TCanvas("canvas","canvas",850/2,1100/2);
00261   canvas->Divide(1,2);
00262   Int_t icanvas=0;
00263 
00264   for ( Int_t ich=0;ich<120;ich+=stride )
00265     {
00266 
00267       canvas->cd(1+icanvas%2);
00268       icanvas++;
00269 
00270       //    TString aname="crate";aname+=crate;aname+=" channels ";aname+=ich;aname+=" to ";aname+=ich+stride-1;
00271       //    TCanvas *c = new TCanvas(aname,aname,400,300);
00272 
00273       TString pname="crate";
00274       pname+=crate+1;
00275       pname+="_ch";
00276       pname+=ich;
00277       pname+="-";
00278       pname+=ich+stride-1;
00279 
00280       const Char_t *opts[]={"ALP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP","LP"};
00281       
00282       // normalize
00283       Double_t sum[stride];for ( Int_t jj=0;jj<stride;jj++ )sum[jj]=0.;
00284       Double_t max=0.;
00285       for ( Int_t jch=0;jch<stride;jch++ ) // loop over channels in this one graph
00286         {
00287 
00288           Int_t index=ich+jch;
00289           Double_t *Y=towerChanCurves[crate][index]->GetY();
00290           for ( Int_t ipoint=0;ipoint<npoints;ipoint++ ) sum[jch]+=Y[ipoint];
00291           if ( sum[jch]>max ) max=sum[jch];
00292         }
00293       if ( max <= 0. ) continue; // meh?
00294 
00295       TLegend *legend=new TLegend(0.125,0.15,0.325,0.45);
00296       for ( Int_t jch=0;jch<stride;jch++ )
00297         {
00298 
00299           Int_t index=ich+jch;
00300           towerChanCurves[crate][index]->SetMarkerSize(0.75);
00301           // offset X axis of each of these
00302           Double_t *X=towerChanCurves[crate][index]->GetX();
00303           Double_t *Y=towerChanCurves[crate][index]->GetY();
00304           Double_t *EY=towerChanCurves[crate][index]->GetEY();
00305           if ( sum[jch]<= 0. ) continue;
00306           //      std::cout<<"before"<<std::endl;
00307           for ( Int_t ip=0;ip<npoints;ip++ ){
00308             Float_t shift = 0.5+ ((float)jch) - ((float)stride)/2.0;
00309             Double_t yy=Y[ip];
00310             X[ip]-= 0.1*shift;
00311             Y[ip]*=max/sum[jch];
00312             EY[ip]*=max/sum[jch];
00313             //      std::cout << "ip="<<ip<<" y="<<yy<<" y'="<<Y[ip]<<std::endl;
00314           }
00315 
00316           towerChanCurves[crate][index]->Sort();
00317           towerChanCurves[crate][index]->SetMinimum(0.);          
00318           towerChanCurves[crate][index]->SetMarkerColor(38+jch);
00319           towerChanCurves[crate][index]->SetLineColor(38+jch);    
00320           towerChanCurves[crate][index]->Draw(opts[jch]);
00321           
00322           TString label="crate ";label+=crate+1;label+=" chan ";label+=index;
00323           legend->AddEntry(towerChanCurves[crate][index],label,"lp");
00324 
00325         }
00326       legend->Draw();
00327       canvas->Update();
00328 
00329       //      if(doprint)c->Print(pname+".gif");
00330       if ( !(icanvas%2) ){
00331         canvas->Print(fname+"(");
00332         canvas->Clear();
00333         canvas->Divide(1,2);
00334       }
00335 
00336     }
00337   canvas->Print(fname+")");
00338   gSystem->Exec(TString("ps2pdf ")+fname);
00339 
00340 }
00341 
00342 
00343 
00344 // ----------------------------------------------------------------------------
00345 void fillCrates(Int_t ipoint)
00346 {
00347 #if 1
00348   // loop over tower crates
00349   for ( Int_t icr=0;icr<MaxTwCrates;icr++ )
00350     {
00351       Float_t yield = (Float_t)mTowerCrateYield[icr];
00352       Float_t total = (Float_t)mTotalYield;
00353       if ( total > 10.0 ) {
00354         Float_t eyield = TMath::Sqrt(yield);
00355         Float_t etotal = TMath::Sqrt(total);
00356         Float_t r = yield / total;
00357         Float_t e1 = (yield>0)? eyield/yield : 0.0;
00358         Float_t e2 = etotal/total;
00359         Float_t er = r * TMath::Sqrt( e1*e1 + e2*e2 );
00360         towerCrateCurves[icr]->SetPoint(ipoint, mTowerDelay, r );
00361         towerCrateCurves[icr]->SetPointError( ipoint, 0., er );
00362       }
00363       else {
00364         towerCrateCurves[icr]->SetPoint(ipoint, mTowerDelay, -1.0 );
00365         towerCrateCurves[icr]->SetPointError( ipoint, 0., 0. );
00366       }
00367     }
00368   // loop over mapmt crates
00369   for ( Int_t icr=0;icr<MaxMapmtCrates;icr++ )
00370     {
00371       Float_t yield = (Float_t)mMapmtCrateYield[icr];
00372       Float_t total = (Float_t)mTotalYield;
00373       if ( total > 10.0 ) {
00374         Float_t eyield = TMath::Sqrt(yield);
00375         Float_t etotal = TMath::Sqrt(total);
00376         Float_t r = yield / total;
00377         Float_t e1 = (yield>0)? eyield/yield : 0.0;
00378         Float_t e2 = etotal/total;
00379         Float_t er = r * TMath::Sqrt( e1*e1 + e2*e2 );
00380         mapmtCrateCurves[icr]->SetPoint(ipoint, mMapmtDelay, r );
00381         mapmtCrateCurves[icr]->SetPointError( ipoint, 0., er );
00382       }
00383       else {
00384         mapmtCrateCurves[icr]->SetPoint(ipoint, mMapmtDelay, -1. );
00385         mapmtCrateCurves[icr]->SetPointError( ipoint, 0., 0. ); 
00386       }
00387     }
00388 #endif
00389 }
00390 // ----------------------------------------------------------------------------
00391 void fillChannels(Int_t ipoint)
00392 {
00393 
00394 #if 1
00395   // loop over tower crates
00396   for ( Int_t icr=0;icr<MaxTwCrates;icr++ )
00397     {
00398       for ( Int_t ich=0;ich<MaxTwCrateCh;ich++ ) 
00399         {
00400 
00401           Float_t yield = (Float_t)mTowerChanYield[icr][ich];
00402           Float_t total = (Float_t)mTotalYield;
00403           if ( total > 10.0 ) {
00404             Float_t eyield = TMath::Sqrt(yield);
00405             Float_t etotal = TMath::Sqrt(total);
00406             Float_t r = yield / total;
00407             Float_t e1 = (yield>0)? eyield/yield : 0.0;
00408             Float_t e2 = etotal/total;
00409             Float_t er = r * TMath::Sqrt( e1*e1 + e2*e2 );
00410             towerChanCurves[icr][ich]->SetPoint(ipoint, mTowerDelay, r );
00411             towerChanCurves[icr][ich]->SetPointError( ipoint, 0., er );
00412           }
00413           else {
00414             towerChanCurves[icr][ich]->SetPoint(ipoint, mTowerDelay, -1.0 );
00415             towerChanCurves[icr][ich]->SetPointError( ipoint, 0., 0. );
00416           }
00417         }
00418     }
00419 #endif
00420 #if 1
00421   // loop over mapmt crates
00422   for ( Int_t icr=0;icr<MaxMapmtCrates;icr++ )
00423     {
00424       for ( Int_t ich=0;ich<MaxMapmtCrateCh;ich++ ) 
00425         {
00426 
00427           Float_t yield = (Float_t)mMapmtChanYield[icr][ich];
00428           Float_t total = (Float_t)mTotalYield;
00429           if ( total > 10.0 ) {
00430             Float_t eyield = TMath::Sqrt(yield);
00431             Float_t etotal = TMath::Sqrt(total);
00432             Float_t r = yield / total;
00433             Float_t e1 = (yield>0)? eyield/yield : 0.0;
00434             Float_t e2 = etotal/total;
00435             Float_t er = r * TMath::Sqrt( e1*e1 + e2*e2 );
00436             mapmtChanCurves[icr][ich]->SetPoint(ipoint, mMapmtDelay, r );
00437             mapmtChanCurves[icr][ich]->SetPointError( ipoint, 0., er );
00438           }
00439           else {
00440             mapmtChanCurves[icr][ich]->SetPoint(ipoint, mMapmtDelay, -1.0 );
00441             mapmtChanCurves[icr][ich]->SetPointError( ipoint, 0., 0. );
00442           }
00443         }
00444     }
00445 #endif
00446 }
00447 // ----------------------------------------------------------------------------
00448 void initGraphs()
00449 {
00450 
00451   for ( Int_t i=0;i<MaxTwCrates;i++ ){
00452     towerCrateCurves[i] =  new TGraphErrors(npoints);
00453     towerCrateCurves[i]->SetMarkerStyle(20+i);
00454     towerCrateCurves[i]->SetMarkerColor(i+1);
00455     towerCrateCurves[i]->SetLineColor(i+1);
00456 
00457     for ( Int_t j=0;j<MaxTwCrateCh;j++ )
00458       towerChanCurves[i][j]=(TGraphErrors*)towerCrateCurves[i]->Clone();
00459 
00460   }
00461   for ( Int_t i=0;i<MaxMapmtCrates;i++ ){
00462     mapmtCrateCurves[i]= new TGraphErrors(npoints);
00463     mapmtCrateCurves[i]->SetMarkerStyle(20+i%4);
00464     mapmtCrateCurves[i]->SetMarkerColor(1+i%4);
00465     mapmtCrateCurves[i]->SetLineColor(1+i%4);
00466 
00467     for ( Int_t j=0;j<MaxMapmtCrateCh;j++ )
00468       mapmtChanCurves[i][j]=(TGraphErrors*)mapmtCrateCurves[i]->Clone();
00469 
00470   }
00471 
00472   
00473 
00474 
00475 }
00476 // ----------------------------------------------------------------------------
00477 void chainFiles(const Char_t *path)
00478 {
00479   chain=new TChain("timing","Timing summary");
00480 
00481   TFile *tfile = 0;
00482   std::cout << "chaining files in " << path << std::endl;
00483   TSystemDirectory *dir = new TSystemDirectory("dir",path);
00484   
00485   TIter next( dir->GetListOfFiles() );
00486   TObject *file = 0;
00487   while ( file = (TObject*)next() )
00488     {
00489       TString name=file->GetName();
00490 
00491       // sum the event counter histogram
00492       if ( name.Contains(".root") ) {
00493         // open the TFile and
00494         std::cout << " + " << name << std::endl;
00495         // tfile = TFile::Open(name);
00496         chain->Add(name);
00497       }
00498     }
00499 
00500 }
00501 // ----------------------------------------------------------------------------
00502 void setBranches(const Char_t *dir)
00503 {
00504   
00505   chain->SetBranchAddress("mRunNumber",  &mRunNumber );
00506   chain->SetBranchAddress("mTowerDelay", &mTowerDelay );
00507   chain->SetBranchAddress("mMapmtDelay", &mMapmtDelay );
00508 
00509   chain->SetBranchAddress("mTotalYield",      &mTotalYield );
00510   chain->SetBranchAddress("mTowerCrateYield", &mTowerCrateYield );
00511   chain->SetBranchAddress("mMapmtCrateYield", &mMapmtCrateYield );
00512   chain->SetBranchAddress("mTowerChanYield",  &mTowerChanYield  );
00513   chain->SetBranchAddress("mMapmtChanYield",  &mMapmtChanYield  );
00514 
00515   chain->SetBranchAddress("mTowerMin",&mTowerMin);
00516   chain->SetBranchAddress("mTowerMax",&mTowerMax);
00517   chain->SetBranchAddress("mMapmtMin",&mMapmtMin);
00518   chain->SetBranchAddress("mMapmtMax",&mMapmtMax);
00519 
00520 }
00521