Added protection against failed fit, for automatic QA trending.
[u/mrichter/AliRoot.git] / PWGPP / TOF / trending / MakeTrendingTOFQA.C
1 /*
2   fbellini@cern.ch - last update on 21/02/2014
3   Macro to run the TOF QA trending by accessing the std QA output, 
4   to be mainly used with the automatic scripts to fill the QA repository.
5   Launch with 
6   aliroot -l -b -q "MakeTrendingTOFQA.C(\"${fullpath}/QAresults.root\", ${run}, ...) 
7   The macro produces a file containing the tree of trending variables and the main plots.
8   A feature that displays the plots in canvases must be enable when needed.
9 */
10
11 Int_t MakeTrendingTOFQA(TString qafilename,       //full path of the QA output; set IsOnGrid to prepend "alien://"
12                         Int_t runNumber,          // run number
13                         Bool_t isMC=kFALSE,       //MC flag, to disable meaningless checks
14                         Bool_t canvasE = kFALSE,  //enable display plots on canvas and save png
15                         Bool_t IsOnGrid = kFALSE) //set to kTRUE to access files on the grid
16 {   
17   // macro to generate tree with TOF QA trending variables
18   // access qa PWGPP output files  
19   if (!qafilename) {
20     Printf("Error - Invalid input file");
21     return 1;
22   }
23
24   /*set graphic style*/
25   gStyle->SetCanvasColor(kWhite);
26   gStyle->SetFrameFillColor(kWhite);
27   gStyle->SetFrameBorderMode(0);
28   gStyle->SetCanvasBorderMode(0);
29   gStyle->SetTitleFillColor(kWhite);
30   gStyle->SetTitleBorderSize(0)  ;
31   gStyle->SetTitleFont(42);
32   gStyle->SetTextFont(42);
33   gStyle->SetStatColor(kWhite); 
34   gStyle->SetStatBorderSize(1);
35   TGaxis::SetMaxDigits(3);
36   gStyle->SetOptStat(10);
37
38   char defaultQAoutput[30]="QAresults.root";
39   char * treePostFileName="trending.root";
40   
41   if (IsOnGrid) TGrid::Connect("alien://");
42   TFile * fin = TFile::Open(qafilename,"r");
43   if (!fin) {
44     Printf("ERROR: QA output not found. Exiting...\n");
45     return -1;
46   } else {
47     Printf("INFO: QA output file %s open. \n",fin->GetName());
48   }
49   
50   //access histograms lists
51   char tofQAdirName[15]="TOF_Performance";
52   char genListName[15]="cGeneralTOFqa";
53   char t0ListName[15]="cTimeZeroTOFqa";
54   char pidListName[15]="cPIDTOFqa"; 
55   char posListName[15]="cPositiveTOFqa";
56   char negListName[15]="cNegativeTOFqa";
57   
58   TDirectoryFile * tofQAdir=(TDirectoryFile*)fin->Get(tofQAdirName);
59   if (!tofQAdir) {
60     Printf("ERROR: TOF QA directory not present in input file.\n");
61     return -1;
62   }
63   TList * generalList=(TList*)tofQAdir->Get(genListName);
64   TList  *timeZeroList=(TList*)tofQAdir->Get(t0ListName);
65   TList  *pidList=(TList*)tofQAdir->Get(pidListName);
66   TList  *posList=(TList*)tofQAdir->Get(posListName);
67   TList  *negList=(TList*)tofQAdir->Get(negListName);
68   
69   if (!generalList) Printf("WARNING: general QA histograms absent or not accessible\n");
70   if (!timeZeroList) Printf("WARNING: timeZero QA histograms absent or not accessible\n");
71   if (!pidList) Printf("WARNING: PID QA histograms absent or not accessible\n");
72   if (!posList) Printf("WARNING: general QA histograms for positive tracks absent or not accessible\n");
73   if (!negList) Printf("WARNING: general QA histograms for negative tracks absent or not accessible\n");
74   
75   if ( (!generalList) && (!timeZeroList) && (!pidList) ){
76     printf("ERROR: no QA available \n");
77     return -1;
78   }
79   
80   Printf(":::: Getting post-analysis info for run %i",runNumber);
81   TFile * trendFile = new TFile(treePostFileName,"recreate");
82
83   Double_t avTime=-9999., peakTime=-9999., spreadTime=-9999., peakTimeErr=-9999., spreadTimeErr=-9999., negTimeRatio=-9999.,
84     avRawTime=-9999., peakRawTime=-9999., spreadRawTime=-9999., peakRawTimeErr=-9999., spreadRawTimeErr=-9999., avTot=-9999., peakTot=-9999.,spreadTot=-9999.,  peakTotErr=-9999.,spreadTotErr=-9999.,
85     orphansRatio=-9999., avL=-9999., negLratio=-9999.,
86     effPt1=-9999., effPt2=-9999., matchEffLinFit1Gev=-9999.,matchEffLinFit1GevErr=-9999.;
87   
88   Double_t avPiDiffTime=-9999.,peakPiDiffTime=-9999., spreadPiDiffTime=-9999.,peakPiDiffTimeErr=-9999., spreadPiDiffTimeErr=-9999.;
89   
90   Double_t avT0A=-9999.,peakT0A=-9999., spreadT0A=-9999.,peakT0AErr=-9999., spreadT0AErr=-9999.;
91   Double_t avT0C=-9999.,peakT0C=-9999., spreadT0C=-9999.,peakT0CErr=-9999., spreadT0CErr=-9999.;
92   Double_t avT0AC=-9999.,peakT0AC=-9999., spreadT0AC=-9999.,peakT0ACErr=-9999., spreadT0ACErr=-9999.;
93   Double_t avT0res=-9999.,peakT0res=-9999., spreadT0res=-9999.,peakT0resErr=-9999., spreadT0resErr=-9999.;
94   Double_t avT0fillRes=-9999.;
95
96   Float_t avMulti=0;
97   Float_t fractionEventsWHits=-9999.;
98   /*number of good (HW ok && efficient && !noisy) TOF channels from OCDB*/
99   Double_t goodChannelRatio=0.0;
100
101   TTree * ttree=new TTree("trending","tree of trending variables");
102   ttree->Branch("run",&runNumber,"run/I");
103   ttree->Branch("avMulti",&avMulti,"avMulti/F");
104   ttree->Branch("fractionEventsWHits",&fractionEventsWHits,"fractionEventsWHits/F");
105   ttree->Branch("goodChannelsRatio",&goodChannelRatio,"goodChannelRatio/D");
106   ttree->Branch("avTime",&avTime,"avTime/D"); //mean time
107   ttree->Branch("peakTime",&peakTime,"peakTime/D"); //main peak time after fit
108   ttree->Branch("spreadTime",&spreadTime,"spreadTime/D"); //spread of main peak of time after fit
109   ttree->Branch("peakTimeErr",&peakTimeErr,"peakTimeErr/D"); //main peak time after fit error
110   ttree->Branch("spreadTimeErr",&spreadTimeErr,"spreadTimeErr/D"); //spread of main peak of time after fit error
111   ttree->Branch("negTimeRatio",&negTimeRatio,"negTimeRatio/D"); //negative time ratio
112   
113   ttree->Branch("avRawTime",&avRawTime,"avRawTime/D"); //mean raw time
114   ttree->Branch("peakRawTime",&peakRawTime,"peakRawTime/D"); //mean peak of RAW TIME after fit
115   ttree->Branch("spreadRawTime",&spreadRawTime,"spreadRawTime/D"); //spread of main peak of raw time after fit
116   ttree->Branch("peakRawTimeErr",&peakRawTimeErr,"peakRawTimeErr/D"); //main peak raw  time after fit error
117   ttree->Branch("spreadRawTimeErr",&spreadRawTimeErr,"spreadRawTimeErr/D"); //spread of  raw main peak of time after fit error
118   
119   ttree->Branch("avTot",&avTot,"avTot/D"); //main peak tot
120   ttree->Branch("peakTot",&peakTot,"peakTot/D"); // main peak of tot after fit
121   ttree->Branch("spreadTot",&spreadTot,"spreadTot/D"); //spread of main peak of tot after fit
122   ttree->Branch("peakTotErr",&peakTotErr,"peakTotErr/D"); // main peak of tot after fit
123   ttree->Branch("spreadTotErr",&spreadTotErr,"spreadTotErr/D"); //spread of main peak of tot after fit
124   
125   ttree->Branch("orphansRatio",&orphansRatio,"orphansRatio/D"); //orphans ratio
126
127   ttree->Branch("avL",&avL,"avL/D"); //mean track length
128   ttree->Branch("negLratio",&negLratio,"negLratio/D");//ratio of tracks with track length <350 cm
129   ttree->Branch("effPt1",&effPt1,"effPt1/D");//matching eff at 1 GeV/c
130   ttree->Branch("effPt2",&effPt2,"effPt2/D"); //matching eff at 2 GeV/c
131   ttree->Branch("matchEffLinFit1Gev",&matchEffLinFit1Gev,"matchEffLinFit1Gev/D");//matching eff fit param 
132   ttree->Branch("matchEffLinFit1GevErr",&matchEffLinFit1GevErr,"matchEffLinFit1GevErr/D");////matching eff fit param error
133   
134   ttree->Branch("avPiDiffTime",&avPiDiffTime,"avPiDiffTime/D"); //mean t-texp
135   ttree->Branch("peakPiDiffTime",&peakPiDiffTime,"peakPiDiffTime/D"); //main peak t-texp after fit
136   ttree->Branch("spreadPiDiffTime",&spreadPiDiffTime,"spreadPiDiffTime/D"); //spread of main peak t-texp after fit
137   ttree->Branch("peakPiDiffTimeErr",&peakPiDiffTimeErr,"peakPiDiffTimeErr/D"); //main peak t-texp after fit error
138   ttree->Branch("spreadPiDiffTimeErr",&spreadPiDiffTimeErr,"spreadPiDiffTimeErr/D"); //spread of main peak of t-texp after fit error
139
140   ttree->Branch("avT0A",&avT0A,"avT0A/D"); //main peak t0A
141   ttree->Branch("peakT0A",&peakT0A,"peakT0A/D"); // main peak of t0A after fit
142   ttree->Branch("spreadT0A",&spreadT0A,"spreadTot/D"); //spread of main peak of t0A after fit
143   ttree->Branch("peakT0AErr",&peakT0AErr,"peakT0AErr/D"); // main peak of t0A after fit
144   ttree->Branch("spreadT0AErr",&spreadT0AErr,"spreadT0AErr/D"); //spread of main peak of t0A after fit
145
146   ttree->Branch("avT0C",&avT0C,"avT0C/D"); //main peak t0C
147   ttree->Branch("peakT0C",&peakT0C,"peakT0C/D"); // main peak of t0C after fit
148   ttree->Branch("spreadT0C",&spreadT0C,"spreadT0C/D"); //spread of main peak of t0C after fit
149   ttree->Branch("peakT0CErr",&peakT0CErr,"peakT0CErr/D"); // main peak of t0C after fit
150   ttree->Branch("spreadT0CErr",&spreadT0CErr,"spreadT0CErr/D"); //spread of main peak of t0C after fit
151  
152   ttree->Branch("avT0AC",&avT0AC,"avT0AC/D"); //main peak t0AC
153   ttree->Branch("peakT0AC",&peakT0AC,"peakT0AC/D"); // main peak of t0AC after fit
154   ttree->Branch("spreadT0AC",&spreadT0AC,"spreadT0AC/D"); //spread of main peak of t0AC after fit
155   ttree->Branch("peakT0ACErr",&peakT0ACErr,"peakT0ACErr/D"); // main peak of t0AC after fit
156   ttree->Branch("spreadT0ACErr",&spreadT0ACErr,"spreadT0ACErr/D"); //spread of main peak of t0AC after fit
157  
158   ttree->Branch("avT0res",&avT0res,"avT0res/D"); //main peak t0AC
159   ttree->Branch("peakT0res",&peakT0res,"peakT0res/D"); // main peak of t0AC after fit
160   ttree->Branch("spreadT0res",&spreadT0res,"spreadT0res/D"); //spread of main peak of t0AC after fit
161   ttree->Branch("peakT0resErr",&peakT0resErr,"peakT0resErr/D"); // main peak of t0AC after fit
162   ttree->Branch("spreadT0resErr",&spreadT0resErr,"spreadT0resErr/D"); //spread of main peak of t0AC after fit
163   ttree->Branch("avT0fillRes",&avT0fillRes,"avT0fillRes/D"); //t0 fill res
164
165   //save quantities for trending
166   goodChannelRatio=(Double_t)GetGoodTOFChannelsRatio(runNumber);
167         
168   //--------------------------------- Multiplicity ----------------------------------//
169
170   TH1F * hMulti = (TH1F*) generalList->FindObject("hTOFmatchedPerEvt");
171   TH1F* hFractionEventsWhits = new TH1F("hFractionEventsWhits","hFractionEventsWhits;fraction of events with hits (%)",200,0.,100.);
172   Float_t fraction=0.0;
173   if (hMulti->GetEntries()>0.0) {
174     fraction = ((Float_t) hMulti->GetBinContent(1))/((Float_t) hMulti->GetEntries());
175     avMulti = hMulti->GetMean();
176   } else fraction=0.0;
177   hFractionEventsWhits->Fill(fraction*100.);
178   
179   //--------------------------------- T0F signal ----------------------------------//
180   TH1F * hRawTime = (TH1F*)generalList->FindObject("hTOFmatchedESDrawTime");
181   if ((hRawTime)&&(hRawTime->GetEntries()>0)){
182     avRawTime=hRawTime->GetMean();
183     if (!isMC){
184       hRawTime->Fit("landau","RQ0","",200.,250.);
185       if (hRawTime->GetFunction("landau")) {
186         peakRawTime=(hRawTime->GetFunction("landau"))->GetParameter(1);
187         spreadRawTime=(hRawTime->GetFunction("landau"))->GetParameter(2);
188         peakRawTimeErr=(hRawTime->GetFunction("landau"))->GetParError(1);
189         spreadRawTimeErr=(hRawTime->GetFunction("landau"))->GetParError(2);
190       //        printf("Main peak raw time (landau): mean = %f +- %f\n",peakTime,peakTimeErr );
191       // printf("Main peak raw time (landau): spread = %f +- %f\n",spreadRawTime,spreadRawTimeErr );
192       }
193     } else {
194       printf("Reminder: Raw time not available in MC simulated data.");
195     }
196   }
197   MakeUpHisto(hRawTime, "matched tracks", 21, kGreen+2);
198   hRawTime->Rebin(2);
199   
200   TH1F * hTime = (TH1F*)generalList->FindObject("hTOFmatchedESDtime");
201   if ((hTime)&&(hTime->GetEntries()>0)) {
202     avTime=hTime->GetMean();
203     hTime->Fit("landau","RQ0","",0.,50.);
204     if (hTime->GetFunction("landau")) {
205       peakTime=(hTime->GetFunction("landau"))->GetParameter(1);
206       spreadTime=(hTime->GetFunction("landau"))->GetParameter(2);
207       peakTimeErr=(hTime->GetFunction("landau"))->GetParError(1);
208       spreadTimeErr=(hTime->GetFunction("landau"))->GetParError(2);
209       negTimeRatio=((Double_t)hTime->Integral(1,3)*100.)/((Double_t)hTime->Integral());
210     // printf("Main peak time (landau): mean = %f +- %f\n",peakTime,peakTimeErr );
211     // printf("Main peak time (landau): spread = %f +- %f\n",spreadTime,spreadTimeErr );
212     // printf("Ratio of tracks with time<12.5 ns / total = %f\n",negTimeRatio );
213     }
214     MakeUpHisto(hTime, "matched tracks", 20, kBlue+2);
215     hTime->Rebin(2);
216     
217     TLegend *lTime = new TLegend(0.7125881,0.6052519,0.979435,0.7408306,NULL,"brNDC");
218     lTime->SetTextSize(0.04281433);
219     lTime->AddEntry(hRawTime, "raw","L");
220     lTime->AddEntry(hTime, "ESD","L"); 
221     lTime->SetFillColor(kWhite);
222     lTime->SetShadowColor(0);
223   }
224       
225   TH1F * hTot = (TH1F*)generalList->FindObject("hTOFmatchedESDToT");
226   if ((hTot)&&(hTot->GetEntries()>0)){
227     avTot=hTot->GetMean();
228     hTot->Fit("gaus","","",0.,50.);
229     if (hTot->GetFunction("gaus")) {
230       peakTot=(hTot->GetFunction("gaus"))->GetParameter(1);
231       spreadTot=(hTot->GetFunction("gaus"))->GetParameter(2);
232       peakTotErr=(hTot->GetFunction("gaus"))->GetParError(1);
233       spreadTotErr=(hTot->GetFunction("gaus"))->GetParError(2);
234     // printf("Main peak ToT (gaus): mean = %f +- %f\n",peakTot,peakTotErr );
235     // printf("Main peak ToT (gaus): spread = %f +- %f\n",spreadTot,spreadTotErr );     
236     }
237   }      
238   MakeUpHisto(hTot, "matched tracks", 8, kViolet-3);
239   
240   char orphansTxt[200];
241   if (hTot->GetEntries()>1){
242     orphansRatio=((Float_t) hTot->GetBinContent(1))/((Float_t) hTot->GetEntries()) ;
243   }
244   sprintf(orphansTxt,"orphans/matched = %4.2f%%",orphansRatio*100.);
245   TPaveText *tOrphans = new TPaveText(0.38,0.63,0.88,0.7, "NDC");
246   tOrphans->SetBorderSize(0);
247   tOrphans->SetTextSize(0.045);
248   tOrphans->SetFillColor(0); //white background
249   tOrphans->SetTextAlign(12);
250   tOrphans->SetTextColor(kViolet-3);
251   tOrphans->AddText(orphansTxt);
252   
253   TH1F * hL=(TH1F*)generalList->FindObject("hTOFmatchedESDtrkLength");
254   char negLengthTxt[200];
255   if (hL->GetEntries()>0){
256     avL=hL->GetMean();
257     negLratio=(hL->Integral(1,750))/((Float_t) hL->GetEntries()) ;
258   }
259   MakeUpHisto(hL, "matched tracks", 1, kBlue+2);
260   sprintf(negLengthTxt,"trk with L<350cm /matched = %4.2f%%", negLratio*100.);
261   TPaveText *tLength = new TPaveText(0.15,0.83,0.65,0.87, "NDC");
262   tLength->SetBorderSize(0);
263   tLength->SetTextSize(0.04);
264   tLength->SetFillColor(0); //white background
265   tLength->SetTextAlign(11);
266   tLength->SetTextColor(kOrange-3);
267   tLength->AddText(negLengthTxt);
268
269   //--------------------------------- residuals -------------------------------------//
270   TH2F* hDxPos4profile = (TH2F*) generalList->FindObject("hTOFmatchedDxVsPtPos");
271   TH2F* hDxNeg4profile = (TH2F*) generalList->FindObject("hTOFmatchedDxVsPtNeg");    
272   const Int_t ybinMin = 0;
273   const Int_t ybinMax = hDxPos4profile->GetYaxis()->GetNbins() ;
274   TProfile * profDxPos = (TProfile*)hDxPos4profile->ProfileX("profDxPos", ybinMin,ybinMax); 
275   profDxPos->SetLineWidth(2);
276   TProfile * profDxNeg = (TProfile*)hDxNeg4profile->ProfileX("profDxNeg", ybinMin, ybinMax); 
277   profDxNeg->SetLineWidth(2);  
278  
279   TH1 *profRatioPosOverNegDx = (TH1*) profDxPos->Clone();
280   profRatioPosOverNegDx->SetName("profRatioPosOverNegDx");
281   profRatioPosOverNegDx->Divide((TH1*) profDxNeg);
282   profRatioPosOverNegDx->GetYaxis()->SetRangeUser(-5.,5.);
283   profRatioPosOverNegDx->GetXaxis()->SetRangeUser(0.,2.);
284
285   TH2F* hTOFmatchedDzVsStrip = (TH2F*)generalList->FindObject("hTOFmatchedDzVsStrip");
286
287   // fout->cd();
288   // hTOFmatchedDzVsStrip->Write();
289   // hDxPos4profile->Write();
290   // hDxNeg4profile->Write();
291   // profDxPos->Write();
292   // profDxNeg->Write();
293   // profRatioPosOverNegDx->Write();
294  
295   //--------------------------------- matching eff ----------------------------------//
296   //matching as function of pT
297   TH1F * hMatchingVsPt = new TH1F("hMatchingVsPt","Matching probability vs. Pt; Pt(GeV/c); matching probability", 50, 0., 5. );
298   TH1F * hDenom = (TH1F*)generalList->FindObject("hESDprimaryTrackPt"); 
299   if (hDenom) {  
300     hMatchingVsPt=(TH1F*) generalList->FindObject("hTOFmatchedESDPt")->Clone(); 
301     hMatchingVsPt->Rebin(5);
302     hDenom->Rebin(5);
303     hMatchingVsPt->Divide(hDenom);
304     hMatchingVsPt->GetYaxis()->SetTitle("matching efficiency");
305     hMatchingVsPt->SetTitle("TOF matching efficiency as function of transverse momentum");
306     hMatchingVsPt->GetYaxis()->SetRangeUser(0,1.2); 
307   }
308
309   if (hMatchingVsPt->GetEntries()>0){
310     hMatchingVsPt->Fit("pol0","","",1.0,10.);
311     hMatchingVsPt->Draw();
312     if (hMatchingVsPt->GetFunction("pol0")){
313       matchEffLinFit1Gev=(hMatchingVsPt->GetFunction("pol0"))->GetParameter(0);
314       matchEffLinFit1GevErr=(hMatchingVsPt->GetFunction("pol0"))->GetParError(0);       
315       //printf("Matching efficiency fit param is %f +- %f\n",matchEffLinFit1Gev,matchEffLinFit1GevErr );
316     }
317   } else {
318     printf("WARNING: matching efficiency plot has 0 entries. Skipped!\n");
319   }
320   MakeUpHisto(hMatchingVsPt, "efficiency", 1, kBlue+2);
321
322   //matching as function of eta
323   TH1F * hMatchingVsEta =new TH1F("hMatchingVsEta","Matching probability vs. #\Eta; #\Eta; matching probability", 20, -1., 1.);
324   hDenom->Clear();
325   hDenom=(TH1F*)generalList->FindObject("hTOFprimaryESDeta"); 
326   if (hDenom) {  
327     hMatchingVsEta=(TH1F*) generalList->FindObject("hTOFmatchedESDeta")->Clone(); 
328     hMatchingVsEta->Rebin(5);
329     hDenom->Rebin(5);
330     hMatchingVsEta->Divide(hDenom);
331     hMatchingVsEta->GetXaxis()->SetRangeUser(-1,1);
332     hMatchingVsEta->GetYaxis()->SetTitle("matching efficiency");
333     hMatchingVsEta->GetYaxis()->SetRangeUser(0,1.2);
334     hMatchingVsEta->SetTitle("TOF matching efficiency as function of pseudorapidity");
335   }
336   MakeUpHisto(hMatchingVsEta, "efficiency", 1, kBlue+2);
337
338   
339   //matching as function of phi
340   TH1F * hMatchingVsPhi = new TH1F("hMatchingVsPhi","Matching probability vs. Phi; Phi(rad); matching probability", 628, 0., 6.28);
341   hDenom->Clear();
342   hDenom=(TH1F*)generalList->FindObject("hTOFprimaryESDphi");  
343   if (hDenom) {  
344     hMatchingVsPhi=(TH1F*) generalList->FindObject("hTOFmatchedESDphi")->Clone(); 
345     hMatchingVsPhi->Rebin(2);
346     hDenom->Rebin(2);    
347     hMatchingVsPhi->Divide(hDenom);
348     hMatchingVsPhi->GetYaxis()->SetTitle("matching efficiency");
349     hMatchingVsPhi->SetTitle("TOF matching efficiency as function of phi");
350     hMatchingVsPhi->GetYaxis()->SetRangeUser(0,1.2);
351   }
352   MakeUpHisto(hMatchingVsPhi, "efficiency", 1, kBlue+2);
353
354   // fout->cd();
355   // hMatchingVsPt->Write();
356   // hMatchingVsEta->Write();
357   // hMatchingVsPhi->Write();
358   
359
360   //--------------------------------- t-texp ----------------------------------//
361   TH2F * hBetaP=(TH2F*)pidList->FindObject("hTOFmatchedESDpVsBeta");
362   if (hBetaP) hBetaP->GetYaxis()->SetRangeUser(0.,1.2);
363   
364   TH1F * hMass=(TH1F*)pidList->FindObject("hTOFmatchedMass");
365   MakeUpHisto(hMass, "tracks", 1, kBlue+2);
366   // hMass->SetFillColor(kAzure+10);
367   // hMass->SetFillStyle(1001);
368   hMass->Rebin(2);
369   
370   //pions
371   TH2F * hDiffTimeT0TOFPion1GeV=(TH2F*)pidList->FindObject("hTOFmatchedTimePion1GeV"); 
372   TH1F * hPionDiff=(TH1F*)pidList->FindObject("hTOFmatchedExpTimePi"); 
373   TH1F * hKaonDiff=(TH1F*)pidList->FindObject("hTOFmatchedExpTimeKa"); 
374   TH1F * hProtonDiff=(TH1F*)pidList->FindObject("hTOFmatchedExpTimePro"); 
375   if ((hPionDiff)&&(hPionDiff->GetEntries()>0)) {
376     avPiDiffTime=hPionDiff->GetMean();
377     hPionDiff->Fit("gaus","","",-1000.,500.);
378     if (hPionDiff->GetFunction("gaus")){
379       peakPiDiffTime=(hPionDiff->GetFunction("gaus"))->GetParameter(1);
380       spreadPiDiffTime=(hPionDiff->GetFunction("gaus"))->GetParameter(2);
381       peakPiDiffTimeErr=(hPionDiff->GetFunction("gaus"))->GetParError(1);
382       spreadPiDiffTimeErr=(hPionDiff->GetFunction("gaus"))->GetParError(2);
383       // printf("Main peak t-t_exp (gaus): mean = %f +- %f\n",peakPiDiffTime,peakPiDiffTimeErr );
384       // printf("Main peak t-t_exp (gaus): spread = %f +- %f\n",spreadPiDiffTime,spreadPiDiffTimeErr );
385     }
386   }
387   
388   TH2F * hDiffTimePi=(TH2F*)pidList->FindObject("hTOFmatchedExpTimePiVsP"); 
389   hDiffTimePi->GetYaxis()->SetRangeUser(-5000.,5000.);
390   hDiffTimePi->SetTitle("PIONS t-t_{exp,#pi} (from tracking) vs. P");
391
392   // TProfile * profDiffTimePi = (TProfile*)hDiffTimePi->ProfileX("profDiffTimePi", 490, 510); 
393   // if (profDiffTimePi){
394   //   profDiffTimePi->SetLineWidth(2);
395   //   profDiffTimePi->SetLineColor(kRed+2); 
396   // }
397
398   //Kaon
399   TH2F * hDiffTimeKa=(TH2F*)pidList->FindObject("hTOFmatchedExpTimeKaVsP");  
400   hDiffTimeKa->SetTitle("KAONS t-t_{exp,K} (from tracking) vs. P");
401   hDiffTimeKa->GetYaxis()->SetRangeUser(-5000.,5000.);
402   
403   // TProfile * profDiffTimeKa = (TProfile*)hDiffTimeKa->ProfileX("profDiffTimeKa", 490, 510); 
404   // if (profDiffTimeKa) {
405   //   profDiffTimeKa->SetLineWidth(2);
406   //   profDiffTimeKa->SetLineColor(kBlue);  
407   // }
408
409   //Protons
410   TH2F * hDiffTimePro=(TH2F*)pidList->FindObject("hTOFmatchedExpTimeProVsP"); 
411   hDiffTimePro->SetTitle("PROTONS t-t_{exp,p} (from tracking) vs. P");
412   hDiffTimePro->GetYaxis()->SetRangeUser(-5000.,5000.);
413   
414   // TProfile * profDiffTimePro = (TProfile*)hDiffTimePro->ProfileX("profDiffTimePro", 490, 510); 
415   // if (profDiffTimePro) {
416   //   profDiffTimePro->SetLineWidth(2);
417   //   profDiffTimePro->SetLineColor(kGreen+2);  
418   // }
419
420   /*
421   TH2F * hDiffTimePiTh=(TH2F*)pidList->FindObject("hTOFtheoreticalExpTimePiVsP");  
422   hDiffTimePiTh->GetYaxis()->SetRangeUser(-2000.,2000.); 
423
424   TProfile * profDiffTimePiTh = (TProfile*)hDiffTimePiTh->ProfileX("profDiffTimePiTh", 490, 510);
425   if (profDiffTimePiTh) {
426     profDiffTimePiTh->SetLineWidth(2);
427     profDiffTimePiTh->SetLineColor(kRed+2);
428   }
429
430   TH2F * hDiffTimeKaTh=(TH2F*)pidList->FindObject("hTOFtheoreticalExpTimeKaVsP"); 
431   hDiffTimeKaTh->GetYaxis()->SetRangeUser(-2000.,2000.);
432
433   TProfile * profDiffTimeKaTh = (TProfile*)hDiffTimeKaTh->ProfileX("profDiffTimeKaTh", 490, 510); 
434   if (profDiffTimeKaTh) {
435     profDiffTimeKaTh->SetLineWidth(2);
436     profDiffTimeKaTh->SetLineColor(kBlue);
437   }
438
439   TH2F * hDiffTimeProTh=(TH2F*)pidList->FindObject("hTOFtheoreticalExpTimePro"); 
440   hDiffTimeProTh->GetYaxis()->SetRangeUser(-2000.,2000.);
441  
442   TProfile * profDiffTimeProTh = (TProfile*)hDiffTimeProTh->ProfileX("profDiffTimeProTh", 490, 510);
443   if (profDiffTimeProTh) {
444     profDiffTimeProTh->SetLineWidth(2);
445     profDiffTimeProTh->SetLineColor(kGreen+2);
446   }
447   TLegend * lPid=new TLegend(0.75,0.75,0.95,0.95,"PID");
448   lPid->AddEntry(profDiffTimePi,"#pi^{#pm}","l");
449   lPid->AddEntry(profDiffTimeKa,"K^{#pm}","l");
450   lPid->AddEntry(profDiffTimePro,"p^{#pm}","l");
451   
452      if (canvasE){
453      TCanvas *cPidPerformanceTh= new TCanvas("cPidPerformanceTh","summary of pid performance - theoretical times",700,700);
454      cPidPerformanceTh->Divide(2,2);
455      cPidPerformanceTh->cd(1);
456      gPad->SetLogz();
457      gPad->SetGridx();
458      gPad->SetGridy();
459      hDiffTimePiTh->Draw("colz");
460      profDiffTimePiTh->Draw("same");
461      cPidPerformanceTh->cd(2);
462      gPad->SetLogz();
463      gPad->SetGridx();
464      gPad->SetGridy();
465      hDiffTimeKaTh->Draw("colz");
466      profDiffTimeKaTh->Draw("same");
467      cPidPerformanceTh->cd(3);
468      gPad->SetLogz();
469      gPad->SetGridx();
470      gPad->SetGridy();
471      hDiffTimeProTh->Draw("colz");
472      profDiffTimeProTh->Draw("same");
473      }
474     fout->cd();
475     hPionDiff->Write();
476     hKaonDiff->Write();
477     hProtonDiff->Write();
478     hBetaP->Write();
479     hMass->Write();
480     hDiffTimeT0TOFPion1GeV->Write();
481     hDiffTimePi->Write();
482     profDiffTimePi->Write();
483     hDiffTimeKa->Write();
484     profDiffTimeKa->Write();
485     hDiffTimePro->Write();
486     profDiffTimePro->Write();
487     //lPid->Draw();
488     hDiffTimePiTh->Write();
489     profDiffTimePiTh->Write();
490     hDiffTimeKaTh->Write();
491     profDiffTimeKaTh->Write();
492     hDiffTimeProTh->Write();
493     profDiffTimeProTh->Write();
494   */
495
496   TH2F * hSigmaPi=(TH2F*)pidList->FindObject("hTOFExpSigmaPi"); 
497   hSigmaPi->GetYaxis()->SetRangeUser(-5.,5.);
498   TProfile * profSigmaPi = (TProfile*)hSigmaPi->ProfileX("profSigmaPi"); 
499   profSigmaPi->SetLineWidth(2);
500   profSigmaPi->SetLineColor(kRed+2); 
501
502   TH2F * hSigmaKa=(TH2F*)pidList->FindObject("hTOFExpSigmaKa"); 
503   hSigmaKa->GetYaxis()->SetRangeUser(-5.,5.);
504   TProfile * profSigmaKa = (TProfile*)hSigmaKa->ProfileX("profSigmaKa"); 
505   profSigmaKa->SetLineWidth(2);
506   profSigmaKa->SetLineColor(kBlue);  
507
508   TH2F * hSigmaPro=(TH2F*)pidList->FindObject("hTOFExpSigmaPro"); 
509   hSigmaPro->GetYaxis()->SetRangeUser(-5.,5.);
510   TProfile * profSigmaPro = (TProfile*)hSigmaPro->ProfileX("profSigmaPro"); 
511   profSigmaPro->SetLineWidth(2);
512   profSigmaPro->SetLineColor(kGreen+2);  
513
514   // fout->cd();
515   // hSigmaPi->Write();
516   // profSigmaPi->Write();
517   // hSigmaKa->Write();
518   // profSigmaKa->Write();
519   // hSigmaPro->Write();
520   // profSigmaPro->Write();
521   
522   //--------------------------------- T0 detector ----------------------------------//
523         
524   TH1F*hT0A=(TH1F*)timeZeroList->FindObject("hEventT0DetA");
525   if ((hT0A)&&(hT0A->GetEntries()>0)) {
526     avT0A = hT0A->GetMean();
527     hT0A->Fit("gaus","RQ0", "", -1000., 1000.);
528     if (hT0A->GetFunction("gaus")) {
529       peakT0A=(hT0A->GetFunction("gaus"))->GetParameter(1);
530       spreadT0A=(hT0A->GetFunction("gaus"))->GetParameter(2);
531       peakT0AErr=(hT0A->GetFunction("gaus"))->GetParError(1);
532       spreadT0AErr=(hT0A->GetFunction("gaus"))->GetParError(2); 
533       // printf("Main peak T0A(gaus): mean = %f +- %f\n",peakT0A,peakT0AErr );
534       // printf("Main peak T0A (gaus): spread = %f +- %f\n",spreadT0A,spreadT0AErr );
535       //add integral of main peak over total
536     }
537   }
538   MakeUpHisto(hT0A, "events", 8, kBlue);
539   hT0A->Rebin(2);
540
541   TH1F*hT0C=(TH1F*)timeZeroList->FindObject("hEventT0DetC");
542   if ((hT0C)&&(hT0C->GetEntries()>0)) {
543     avT0C=hT0C->GetMean();
544     hT0C->Fit("gaus","RQ0","", -1000., 1000.);
545     if (hT0C->GetFunction("gaus")) {
546       peakT0C=(hT0C->GetFunction("gaus"))->GetParameter(1);
547       spreadT0C=(hT0C->GetFunction("gaus"))->GetParameter(2);
548       peakT0CErr=(hT0C->GetFunction("gaus"))->GetParError(1);
549       spreadT0CErr=(hT0C->GetFunction("gaus"))->GetParError(2); 
550       // printf("Main peak T0C(gaus): mean = %f +- %f\n",peakT0C,peakT0CErr );
551       // printf("Main peak T0C (gaus): spread = %f +- %f\n",spreadT0C,spreadT0CErr );
552       //add integral of main peak over total
553     }
554   }
555   MakeUpHisto(hT0C, "events", 8, kGreen+1);
556   hT0C->Rebin(2);
557         
558   TH1F*hT0AC=(TH1F*)timeZeroList->FindObject("hEventT0DetAND");
559   if ((hT0AC)&&(hT0AC->GetEntries()>0)) {
560     avT0AC=hT0AC->GetMean();
561     hT0AC->Fit("gaus","RQ0", "",-1000., 1000.);
562     if (hT0AC->GetFunction("gaus")) {
563       peakT0AC=(hT0AC->GetFunction("gaus"))->GetParameter(1);
564       spreadT0AC=(hT0AC->GetFunction("gaus"))->GetParameter(2);
565       peakT0ACErr=(hT0AC->GetFunction("gaus"))->GetParError(1);
566       spreadT0ACErr=(hT0AC->GetFunction("gaus"))->GetParError(2);       
567       // printf("Main peak T0AC(gaus): mean = %f +- %f\n",peakT0AC,peakT0ACErr );
568       // printf("Main peak T0AC (gaus): spread = %f +- %f\n",spreadT0AC,spreadT0ACErr );         
569     }
570   }
571   MakeUpHisto(hT0AC, "events", 8, kRed+1);
572   hT0AC->Rebin(2);
573   
574   TLegend *lT0 = new TLegend(0.7125881,0.6052519,0.979435,0.7408306,NULL,"brNDC");
575   lT0->SetTextSize(0.041);
576   lT0->AddEntry(hT0AC, "T0 A&C","L");
577   lT0->AddEntry(hT0A, "T0 A","L"); 
578   lT0->AddEntry(hT0C, "T0 C","L");
579   lT0->SetFillColor(kWhite);
580   lT0->SetShadowColor(0);
581   
582   TH1F*hT0res=(TH1F*)timeZeroList->FindObject("hT0DetRes");
583   if ((hT0res)&&(hT0res->GetEntries()>0)) {
584     avT0res=hT0res->GetMean();
585     hT0res->Fit("gaus");
586     if (hT0res->GetFunction("gaus")) {
587       peakT0res=(hT0res->GetFunction("gaus"))->GetParameter(1);
588       spreadT0res=(hT0res->GetFunction("gaus"))->GetParameter(2);
589       peakT0resErr=(hT0res->GetFunction("gaus"))->GetParError(1);
590       spreadT0resErr=(hT0res->GetFunction("gaus"))->GetParError(2);     
591       // printf("Main peak T0res(gaus): mean = %f +- %f\n",peakT0res,peakT0resErr );
592       // printf("Main peak T0res (gaus): spread = %f +- %f\n",spreadT0res,spreadT0resErr );      
593     //add integral of main peak over total
594     }
595   }
596   TH1F*hT0fillRes=(TH1F*)timeZeroList->FindObject("hT0fillRes");
597   if ((hT0fillRes)&&(hT0fillRes->GetEntries()>0)) {
598     avT0fillRes=hT0fillRes->GetMean();
599   }
600   
601   // fout->cd();
602   // hT0AC->Write();
603   // hT0A->Write();
604   // hT0C->Write();
605   // hT0res->Write();
606   // hT0fillRes->Write();
607   // lT0->Write();
608         
609   //Fill tree and save to file
610   ttree->Fill();
611   printf("==============  Saving trending quantities in tree for run %i ===============\n",runNumber);
612   trendFile->cd();
613   ttree->Write();
614   trendFile->Close();
615   
616   if (canvasE){
617     // TString plotDir(Form("Plots_run%d",runNumber));
618     // gSystem->Exec(Form("mkdir %s",plotDir.Data()));
619     TString plotDir(".");
620
621     TCanvas *cTrackProperties= new TCanvas("cTrackProperties","summary of matched tracks properties", 1200, 500);
622     cTrackProperties->Divide(3,1);
623     cTrackProperties->cd(1);
624     gPad->SetLogy();
625     hTime->Draw("");
626     hRawTime ->Draw("same");
627     lTime->Draw();  
628     cTrackProperties->cd(2);
629     gPad->SetLogy();
630     hTot->Draw("");
631     tOrphans->Draw(); 
632     cTrackProperties->cd(3);
633     gPad->SetLogy();
634     hL->Draw("");
635     tLength->Draw(); 
636
637     TCanvas *cResiduals= new TCanvas("residuals","residuals", 900,500);
638     cResiduals->Divide(2,1);
639     cResiduals->cd(1);
640     gPad->SetLogz();
641     hDxPos4profile->GetYaxis()->SetRangeUser(-5.,5.);
642     hDxPos4profile->Draw("colz");
643     profDxPos->SetLineColor(kRed);
644     profDxPos ->Draw("same");
645     cResiduals->cd(2);
646     gPad->SetLogz();
647     hDxNeg4profile->GetYaxis()->SetRangeUser(-5.,5.); 
648     hDxNeg4profile->Draw("colz");
649     profDxNeg->SetLineColor(kBlue);
650     profDxNeg->Draw("same"); 
651
652     TCanvas* cProfile = new TCanvas("cProfile","cProfile",50,50, 750,550);
653     cProfile->cd();
654     gPad->SetLogz();
655     hTOFmatchedDzVsStrip->Draw("colz");
656     Int_t binmin = hTOFmatchedDzVsStrip->GetYaxis()->FindBin(-3);
657     Int_t binmax = hTOFmatchedDzVsStrip->GetYaxis()->FindBin(3);
658     TProfile* hDzProfile = (TProfile*)hTOFmatchedDzVsStrip->ProfileX("hDzProfile",binmin, binmax);
659     hDzProfile->SetLineWidth(3);
660     hDzProfile->Draw("same");
661     
662     TCanvas *cMatchingPerformance= new TCanvas("cMatchingPerformance","summary of matching performance",1200,500);
663     cMatchingPerformance->Divide(3,1);
664     cMatchingPerformance->cd(1);
665     hMatchingVsPt->Draw();
666     cMatchingPerformance->cd(2);
667     hMatchingVsEta->Draw();
668     cMatchingPerformance->cd(3);
669     hMatchingVsPhi->Draw();
670     
671     TCanvas *cPidPerformance= new TCanvas("cPidPerformance","summary of pid performance", 900,500);
672     cPidPerformance->Divide(2,1);
673     cPidPerformance->cd(1);
674     gPad->SetLogz();
675     hBetaP->Draw("colz");   
676     cPidPerformance->cd(2);
677     gPad->SetLogy();
678     hMass->Draw("HIST ");
679
680     TCanvas *cPidPerformance2= new TCanvas("cPidPerformance2","summary of pid performance - expected times", 1200, 500);
681     cPidPerformance2->Divide(3,1);
682     cPidPerformance2->cd(1);
683     gPad->SetLogz();
684     hDiffTimePi->Draw("colz");
685     //profDiffTimePi->Draw("same");
686     cPidPerformance2->cd(2);
687     gPad->SetLogz();
688     hDiffTimeKa->Draw("colz");
689     //profDiffTimeKa->Draw("same");
690     cPidPerformance2->cd(3);
691     gPad->SetLogz();
692     hDiffTimePro->Draw("colz");
693     //profDiffTimePro->Draw("same");
694   
695     TLegend * lSigmaPid=new TLegend(0.75,0.75,0.95,0.95,"#sigma_{PID}");
696     lSigmaPid->AddEntry(profSigmaPi,"#pi^{#pm}","l");
697     lSigmaPid->AddEntry(profSigmaKa,"K^{#pm}","l");
698     lSigmaPid->AddEntry(profSigmaPro,"p^{#pm}","l");      
699     TCanvas *cPidPerformance3= new TCanvas("cPidPerformance3","summary of pid performance - sigmas",1200,500);
700     cPidPerformance3->Divide(3,1);
701     cPidPerformance3->cd(1);
702     gPad->SetLogz();
703     hSigmaPi->Draw("colz");
704     profSigmaPi->Draw("same");
705     cPidPerformance3->cd(2);
706     gPad->SetLogz();
707     hSigmaKa->Draw("colz");
708     profSigmaKa->Draw("same");
709     cPidPerformance3->cd(3);
710     gPad->SetLogz();
711     hSigmaPro->Draw("colz");
712     profSigmaPro->Draw("same");
713   
714     TCanvas *cT0detector= new TCanvas("cT0detector","T0 detector",800,600);
715     cT0detector->Divide(2,1);
716     cT0detector->cd(1);
717     gPad->SetGridx();
718     hT0AC->Draw("");
719     hT0AC->SetTitle("timeZero measured by T0 detector");
720     hT0A ->Draw("same");
721     hT0C ->Draw("same");
722     lT0->Draw();  
723     cT0detector->cd(2);
724     hT0res->Draw();
725     
726     cPidPerformance3->Print(Form("%s/%i_PID_sigmas.png", plotDir.Data(), runNumber));
727     cPidPerformance->Print(Form("%s/%i_PID.png",plotDir.Data(), runNumber));
728     //cPidPerformanceTh->Print(Form("%s/PID_theoreticalTimes.png",plotDir.Data()));
729     cPidPerformance2->Print(Form("%s/%i_PID_ExpTimes.png",plotDir.Data(), runNumber));
730     cMatchingPerformance->Print(Form("%s/%i_Matching.png",plotDir.Data(), runNumber));
731     cTrackProperties->Print(Form("%s/%i_TrackProperties.png",plotDir.Data(), runNumber));
732     cResiduals->Print(Form("%s/%i_Residuals.png",plotDir.Data(), runNumber));
733     cProfile->Print(Form("%s/%i_ProfileDZvsStripNumber.png",plotDir.Data(), runNumber));
734     cT0detector->Print(Form("%s/%i_T0Detector.png",plotDir.Data(), runNumber));
735   }
736  
737   return  0;
738 }
739
740
741 //----------------------------------------------------------
742 Double_t GetGoodTOFChannelsRatio(Int_t run = -1, Bool_t saveMap = kFALSE, TString OCDBstorage = "raw://")
743 {
744   /*
745     It retrieves from OCDB the number of good (= efficient && not noisy && HW ok) TOF channels.
746     Optionally is saves the channel map
747   */
748   if (run<=0) {
749     printf("MakeTrendingTOFqa.C - ERROR in CheckCalibStatus(): invalid run number. Please set a run number.\n"); 
750     return 0.0;
751   }
752   
753   AliCDBManager *cdb = AliCDBManager::Instance();
754   cdb->SetDefaultStorage(OCDBstorage.Data());
755   cdb->SetRun(run);
756   
757   AliCDBEntry *cdbe = cdb->Get("TOF/Calib/Status");
758   if (!cdbe) {
759     printf("MakeTrendingTOFqa.C - ERROR in CheckCalibStatus(): OCDB entry not available. Please, try again.\n");
760     return 0.0;
761   }  
762
763   AliTOFChannelOnlineStatusArray *array = (AliTOFChannelOnlineStatusArray *)cdbe->GetObject();
764   TH2F *hOkMap = new TH2F("hOkMap", "Ok map (!noisy & !problematic & efficient);sector;strip", 72, 0., 18., 91, 0., 91.);
765
766   AliTOFcalibHisto calibHisto;
767   calibHisto.LoadCalibHisto();
768   AliTOFcalib calib;
769   calib.Init();
770   Int_t sector, sectorStrip, padx, fea;
771   Float_t hitmapx, hitmapy;
772   for (Int_t i = 0; i <  array->GetSize(); i++) {
773     sector = calibHisto.GetCalibMap(AliTOFcalibHisto::kSector, i);
774     sectorStrip = calibHisto.GetCalibMap(AliTOFcalibHisto::kSectorStrip, i);
775     padx = calibHisto.GetCalibMap(AliTOFcalibHisto::kPadX, i);
776     fea = padx / 12;
777     hitmapx = sector + ((Double_t)(3 - fea) + 0.5) / 4.;
778     hitmapy = sectorStrip;
779     if ( !(array->GetNoiseStatus(i) == AliTOFChannelOnlineStatusArray::kTOFNoiseBad)   &&
780          (calib.IsChannelEnabled(i,kTRUE,kTRUE)))
781       hOkMap->Fill(hitmapx,hitmapy);
782   }
783   Int_t nOk=(Int_t) hOkMap->GetEntries();
784   Double_t ratioOk=nOk/152928.;
785   if (saveMap) hOkMap->SaveAs(Form("run%i_OKChannelsMap.root",run));
786   cout << "###    Run " << run << ": TOF channels ok = " << nOk << "/ total 152928 channels = " << ratioOk*100. << "% of whole TOF" << endl;
787   return ratioOk;
788 }
789
790 //----------------------------------------------------------
791 void MakeUpHisto(TH1* histo, TString titleY, Int_t marker=20, Color_t color=kBlue+2)
792 {
793   if (!histo) return;
794   histo->SetMarkerStyle(marker);
795   histo->SetMarkerSize(0.7);
796   histo->SetMarkerColor(color);
797   histo->SetLineColor(color);
798   histo->SetFillColor(kWhite);
799   histo->SetFillStyle(0); 
800   histo->GetYaxis()->SetTitle(titleY.Data());
801   histo->GetYaxis()->SetTitleOffset(1.35);
802   histo->GetXaxis()->SetLabelSize(0.03);
803   
804   return;
805 }