TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / HighLevelQA / AODQAChecks.C
1 /////////////////////////////////////////////////////////////////////////////////////////
2 // AODQAChecks.C
3 //
4 // Performs some run-by-run high-level QA checks on AOD data 
5 //    for runs marked with global quality "good" in the RCT
6 //    for the LHC10h period
7 //
8 // Note that the following macros must be in the current
9 //    directory: CheckEfficiencies.C, CheckIntegratedRawYield.C,
10 //    CheckMatchingEfficiency.C, CheckNSigmaStability.C,
11 //    FindOutliers.C, and PlotAndSave.C
12 //
13 // Written by John Groh, summer student (Summer 2012)
14 // Advisors: Michele Floris and Alexander Kalweit
15 /////////////////////////////////////////////////////////////////////////////////////////
16
17 // for cout
18 #include <iostream>
19
20 // for processing .txt file with run numbers
21 #include <fstream>
22 using namespace std;
23
24 // global constants
25 const Int_t nPart = 3;
26 TString Particle[nPart] = {"Pion", "Kaon", "Proton"};
27 const Int_t nCharge = 2;
28 TString Charge[nCharge] = {"Pos", "Neg"};
29 TString Sign[nCharge] = {"Plus", "Minus"};
30 Int_t Color[nPart] = {1,2,4};
31 Int_t Marker[nPart*nCharge] = {20,21,22,24,25,26};
32 TString Names[nPart*nCharge] = {"#pi^{+}",
33                                 "K^{+}",
34                                 "p",
35                                 "#pi^{-}",
36                                 "K^{-}",
37                                 "#bar{p}"};
38
39 // possible cuts to use -   0    1    2    3
40 Double_t CentCutMin[4] =   {0,  30,  30,  30};
41 Double_t CentCutMax[4] =   {5,  40,  40,  40};
42 Double_t QvecCutMin[4] =   {0,   0,   0, 1.5};
43 Double_t QvecCutMax[4] = {100, 100, 0.4, 100};
44 Double_t EtaMin[4] =    {-0.8,-0.8,-0.8,-0.8};
45 Double_t EtaMax[4] =     {0.8, 0.8, 0.8, 0.8};
46 Double_t Nsigmapid = 3.;
47 UInt_t trkbit = 1024;
48
49 // fixed Pt values chosen for efficiency and matching efficiency plots (GeV)
50 const Float_t FixedPtEff = 0.4;
51 const Float_t FixedPtMatchEff = 0.9;
52
53 // inclusions for functions called
54 #include "CheckIntegratedRawYield.C"
55 #include "CheckNSigmaStability.C"
56 #include "CheckEfficiencies.C"
57 #include "CheckMatchingEfficiency.C"
58 #include "PlotAndSave.C"
59 #include "FindOutliers.C"
60
61 //--------------------
62 // Start of function -
63 //--------------------
64 // useMC selects whether to use Monte Carlo or Data - must be run seperately for each
65 // icut selects which cut to use (2 and 3 have low statistics and cause divide-by-zero errors, so don't use)
66 // nSigmaCut is the size of the deviation with which to determine outliers
67 void AODQAChecks(Bool_t useMC = 1, Int_t icut = 1, const Float_t nSigmaCut = 3)
68 {
69   Printf("\n\n--- Running AODQAChecks() with useMC = %i ---\n\n",useMC);
70
71   gSystem->Load("libCore.so");
72   gSystem->Load("libPhysics.so");
73   gSystem->Load("libTree");
74   gSystem->Load("libMatrix");
75   gSystem->Load("libSTEERBase");
76   gSystem->Load("libESD");
77   gSystem->Load("libAOD");
78   gSystem->Load("libANALYSIS");
79   gSystem->Load("libOADB");
80   gSystem->Load("libANALYSISalice");
81   gSystem->Load("libTender");
82   gSystem->Load("libCORRFW");
83   gSystem->Load("libPWGTools");
84   gSystem->Load("libPWGLFspectra");
85   gSystem->Load("libProof.so");
86   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
87
88   TString fold = "AODQAChecks";
89
90   // get number of runs used
91   Int_t nRunsTemp = 0;
92   Int_t dummy;
93   ifstream runList("output/AODQAChecks/RunList.txt");
94   while (!runList.eof())
95     {
96       runList >> dummy;
97       nRunsTemp++;
98     }
99   runList.close();
100
101   // I know this is silly, but you need a const for array sizes.
102   // Also, it seems to want to read one more line than there is in the file
103   const Int_t nRuns = nRunsTemp - 1;
104
105   // fill an array with the run numbers and another with the file names
106   runList.open("output/AODQAChecks/RunList.txt");
107   Int_t runs[nRuns];
108   TString dataFiles[nRuns];
109   
110   for (Int_t irun=0; irun<nRuns; irun++)
111     {
112       runList >> runs[irun];
113       if (useMC) dataFiles[irun] = Form("output/%s/MC/AnalysisResults%i.root",fold.Data(),runs[irun]);
114       else dataFiles[irun] = Form("output/%s/DATA/AnalysisResults%i.root",fold.Data(),runs[irun]);
115     }
116   runList.close();
117
118   // choose which TDirectory in the .root file to use
119   TString sname;
120   if (useMC)
121     sname = Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
122   else
123     sname = Form("OutputAODSpectraTask_Data_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
124
125   // create a .root file for output
126   TFile * fout = new TFile(Form("results/%s/Res_%s_AODQAChecks.root",fold.Data(),sname.Data()),"RECREATE");
127   TFile * file;
128   TDirectoryFile * Dir;
129   AliSpectraAODTrackCuts * tcuts;
130   AliSpectraAODEventCuts * ecuts;
131   AliSpectraAODHistoManager * hman;
132
133   // histograms for stability of nsigma distribution plots
134   TH1F * TPCnsigMeanTrendPion = new TH1F("TPCnsigMeanTrendPion","",nRuns,0,nRuns);
135   TH1F * TPCnsigMeanTrendKaon = new TH1F("TPCnsigMeanTrendKaon","",nRuns,0,nRuns);
136   TH1F * TPCnsigMeanTrendProton = new TH1F("TPCnsigMeanTrendProton","",nRuns,0,nRuns);
137   TH1F * TPCnsigSigmaTrendPion = new TH1F("TPCnsigSigmaTrendPion","",nRuns,0,nRuns);
138   TH1F * TPCnsigSigmaTrendKaon = new TH1F("TPCnsigSigmaTrendKaon","",nRuns,0,nRuns);
139   TH1F * TPCnsigSigmaTrendProton = new TH1F("TPCnsigSigmaTrendProton","",nRuns,0,nRuns);
140   TH1F * TOFnsigMeanTrendPion = new TH1F("TOFnsigMeanTrendPion","",nRuns,0,nRuns);
141   TH1F * TOFnsigMeanTrendKaon = new TH1F("TOFnsigMeanTrendKaon","",nRuns,0,nRuns);
142   TH1F * TOFnsigMeanTrendProton = new TH1F("TOFnsigMeanTrendProton","",nRuns,0,nRuns);
143   TH1F * TOFnsigSigmaTrendPion = new TH1F("TOFnsigSigmaTrendPion","",nRuns,0,nRuns);
144   TH1F * TOFnsigSigmaTrendKaon = new TH1F("TOFnsigSigmaTrendKaon","",nRuns,0,nRuns);
145   TH1F * TOFnsigSigmaTrendProton = new TH1F("TOFnsigSigmaTrendProton","",nRuns,0,nRuns);
146
147   // histograms for integrated raw yield stability plots
148   TH1F * IntegRawYieldAll = new TH1F("IntegRawYieldAll","",nRuns,0,nRuns);
149   TH1F * IntegRawYieldPiPlus = new TH1F("IntegRawYieldPiPlus","",nRuns,0,nRuns);
150   TH1F * IntegRawYieldPiMinus = new TH1F("IntegRawYieldPiMinus","",nRuns,0,nRuns);
151   TH1F * IntegRawYieldKPlus = new TH1F("IntegRawYieldKPlus","",nRuns,0,nRuns);  
152   TH1F * IntegRawYieldKMinus = new TH1F("IntegRawYieldKMinus","",nRuns,0,nRuns);
153   TH1F * IntegRawYieldProton = new TH1F("IntegRawYieldProton","",nRuns,0,nRuns);
154   TH1F * IntegRawYieldAntiproton = new TH1F("IntegRawYieldAntiproton","",nRuns,0,nRuns);
155
156   // objects for efficiency plots
157   if (useMC)
158     {
159       TCanvas * cEfficienciesAllRuns = new TCanvas("cEfficienciesAllRuns","cEfficienciesAllRuns",50,25,700,500);
160       cEfficienciesAllRuns->Divide(3,2);
161       TH1F * EfficiencyPiPlus = new TH1F("EfficiencyPiPlus","",nRuns,0,nRuns);
162       TH1F * EfficiencyKPlus = new TH1F("EfficiencyKPlus","",nRuns,0,nRuns);
163       TH1F * EfficiencyProton = new TH1F("EfficiencyProton","",nRuns,0,nRuns);
164       TH1F * EfficiencyPiMinus = new TH1F("EfficiencyPiMinus","",nRuns,0,nRuns);
165       TH1F * EfficiencyKMinus = new TH1F("EfficiencyKMinus","",nRuns,0,nRuns);
166       TH1F * EfficiencyAntiproton = new TH1F("EfficiencyAntiproton","",nRuns,0,nRuns);
167     }
168   else
169     {
170       TCanvas * cEfficienciesAllRuns = 0x0;
171       TH1F * EfficiencyPiPlus = 0x0;
172       TH1F * EfficiencyKPlus = 0x0;
173       TH1F * EfficiencyProton = 0x0;
174       TH1F * EfficiencyPiMinus = 0x0;
175       TH1F * EfficiencyKMinus = 0x0;
176       TH1F * EfficiencyAntiproton = 0x0;
177     }
178
179
180   // objects for matching efficiency plots
181   TH1F * MatchEffPos = new TH1F("MatchEffPos","",nRuns,0,nRuns);
182   MatchEffPos->Sumw2();
183   TH1F * MatchEffNeg = new TH1F("MatchEffNeg","",nRuns,0,nRuns);
184   MatchEffNeg->Sumw2();
185
186   // overall loop over data files (1 per run)
187   for (Int_t irun=0; irun<nRuns; irun++)
188     {
189       Printf("\n--- Processing data from run %i ---\n", runs[irun]);
190
191       // open file and print all objects in file
192       file = TFile::Open(dataFiles[irun].Data());     
193       if (!file)
194         {
195           Printf("\n\n!!! ERROR: Could not open file %s !!!\n\n",dataFiles[irun].Data());
196           continue;
197         }
198       file->Print();
199       Printf("sname: %s", sname.Data());
200
201       // choose the right directory, event, and histo managers
202       Dir = (TDirectoryFile*)file->Get(Form("%s",sname.Data()));
203       if (!Dir)
204         {
205           Printf("\n\n!!! ERROR: Dir is a null pointer. Skipping to next file !!!\n\n");
206           continue;
207         }
208       ecuts = (AliSpectraAODEventCuts*)Dir->Get("Event Cuts");
209       if (!ecuts)
210         {
211           Printf("\n\n!!! ERROR: ecuts is a null pointer. Skipping to next file !!!\n\n");
212           continue;
213         }      
214       tcuts = (AliSpectraAODTrackCuts*)Dir->Get("Track Cuts");
215       if (!tcuts)
216         {
217           Printf("\n\n!!! ERROR: tcuts is a null pointer. Skipping to next file !!!\n\n");
218           continue;
219         }      
220       hman = (AliSpectraAODHistoManager*)Dir->Get("SpectraHistos");
221       if (!hman)
222         {
223           Printf("\n\n!!! ERROR: hman is a null pointer. Skipping to next file !!!\n\n");
224           continue;
225         }
226
227       //---------------------------------------------------------------------------------------
228       // Find the trends in the means and sigmas of fits to the peaks of fixed-Pt projections -
229       // of the nsigma distributions of the TPC and TOF as a function of the run number       -
230       //---------------------------------------------------------------------------------------
231       CheckNSigmaStability(hman,
232                            TPCnsigMeanTrendPion,
233                            TPCnsigMeanTrendKaon,
234                            TPCnsigMeanTrendProton, 
235                            TPCnsigSigmaTrendPion,
236                            TPCnsigSigmaTrendKaon,
237                            TPCnsigSigmaTrendProton,
238                            TOFnsigMeanTrendPion,
239                            TOFnsigMeanTrendKaon,
240                            TOFnsigMeanTrendProton,
241                            TOFnsigSigmaTrendPion,
242                            TOFnsigSigmaTrendKaon,
243                            TOFnsigSigmaTrendProton,
244                            runs,
245                            nRuns,
246                            irun,
247                            useMC);
248
249       //---------------------------------------------------------------------------------
250       // plot the stability of the integrated raw yield as a function of the run number -
251       //---------------------------------------------------------------------------------
252       CheckIntegratedRawYield(ecuts,
253                               hman,
254                               IntegRawYieldAll,
255                               IntegRawYieldPiPlus,
256                               IntegRawYieldPiMinus,
257                               IntegRawYieldKPlus,
258                               IntegRawYieldKMinus,
259                               IntegRawYieldProton,
260                               IntegRawYieldAntiproton,
261                               runs,
262                               irun,
263                               nRuns,
264                               Names,
265                               useMC);
266       
267       //-----------------------------------------------------------------------------
268       // find the trend in the MC correction factor as a function of the run number -
269       //-----------------------------------------------------------------------------
270       if (useMC) CheckEfficiencies(hman,
271                                    cEfficienciesAllRuns,
272                                    EfficiencyPiPlus,
273                                    EfficiencyKPlus,
274                                    EfficiencyProton,
275                                    EfficiencyPiMinus,
276                                    EfficiencyKMinus,
277                                    EfficiencyAntiproton,
278                                    FixedPtEff,
279                                    runs,
280                                    nRuns,
281                                    irun);
282
283       //-------------------------------------------------------------------------------
284       // for a fixed Pt, plot the matching efficiency as a function of the run number -
285       //-------------------------------------------------------------------------------
286       CheckMatchingEfficiency(tcuts,
287                               MatchEffPos,
288                               MatchEffNeg,
289                               FixedPtMatchEff,
290                               runs,
291                               irun,
292                               nRuns,
293                               useMC);
294
295       //---------------------------------------------------------------------------
296       // save eta-phi distributions for each run on a seperate page of a pdf file -
297       //---------------------------------------------------------------------------
298       // canvas for temporarily drawing
299       TCanvas * cEtaPhi = new TCanvas("cEtaPhi","cEtaPhi");
300       TH2F * hEtaPhi = (TH2F*)((TH2F*)tcuts->GetHistoEtaPhiHighPt()->Clone("hEtaPhi"));
301       hEtaPhi->Rebin2D(4);
302       hEtaPhi->Scale(1./hEtaPhi->GetEntries());
303       gPad->SetGridy();
304       gPad->SetGridx();
305       hEtaPhi->DrawClone("COLZ");
306
307       // legend
308       TLegend * lEtaPhi = new TLegend(0.4,0.85,.6,.95);
309       lEtaPhi->SetFillColor(0);
310       lEtaPhi->AddEntry(hEtaPhi,Form("Run %i",runs[irun]),"");
311       if (useMC) lEtaPhi->AddEntry(hEtaPhi,"MC","");
312       else lEtaPhi->AddEntry(hEtaPhi,"DATA","");
313       lEtaPhi->DrawClone();
314
315       // save to a pdf and close the temporary canvas
316       if (useMC)
317         {
318           if (irun == 0) cEtaPhi->SaveAs("Plots/MC/EtaPhi.pdf(","pdf");
319           else if (irun < nRuns-1) cEtaPhi->SaveAs("Plots/MC/EtaPhi.pdf","pdf");
320           else if (irun == nRuns-1) cEtaPhi->SaveAs("Plots/MC/EtaPhi.pdf)","pdf");
321           cEtaPhi->Close();
322         }
323       else
324         {
325           if (irun == 0) cEtaPhi->SaveAs("Plots/DATA/EtaPhi.pdf(","pdf");
326           else if (irun < nRuns-1) cEtaPhi->SaveAs("Plots/DATA/EtaPhi.pdf","pdf");
327           else if (irun == nRuns-1) cEtaPhi->SaveAs("Plots/DATA/EtaPhi.pdf)","pdf");
328           cEtaPhi->Close();
329         }
330       
331       // close the file and clean up before the next iteration
332       file->Close();
333       //delete file;
334       //delete Dir;
335       //delete hman;
336
337       Printf("\n--- Done Processing data from run %i ---\n", runs[irun]);
338
339     }
340
341   //------------------------------------------------------------
342   // Plot all the results and save them to the output file     -
343   //------------------------------------------------------------
344   PlotAndSave(runs,
345               nRuns,
346               useMC,
347               fout,
348               
349               TPCnsigMeanTrendPion,
350               TPCnsigMeanTrendKaon,
351               TPCnsigMeanTrendProton,
352               TPCnsigSigmaTrendPion,
353               TPCnsigSigmaTrendKaon,
354               TPCnsigSigmaTrendProton,
355               TOFnsigMeanTrendPion,
356               TOFnsigMeanTrendKaon,
357               TOFnsigMeanTrendProton,
358               TOFnsigSigmaTrendPion,
359               TOFnsigSigmaTrendKaon,
360               TOFnsigSigmaTrendProton,
361               
362               IntegRawYieldAll,
363               IntegRawYieldPiPlus,
364               IntegRawYieldKPlus,
365               IntegRawYieldProton,
366               IntegRawYieldPiMinus,
367               IntegRawYieldKMinus,
368               IntegRawYieldAntiproton,
369               
370               EfficiencyPiPlus,
371               EfficiencyKPlus,
372               EfficiencyProton,
373               EfficiencyPiMinus,
374               EfficiencyKMinus,
375               EfficiencyAntiproton,
376               
377               MatchEffPos,
378               MatchEffNeg);
379
380   //------------------------------------------------------------------------------
381   // Find the outliers from the efficiency, raw yield, matching efficiency, and  -
382   // nsigma fit plots and print the corresponding run numbers to the console.    -
383   // Also creates a .txt file in the results/ directory with info about outliers -
384   //------------------------------------------------------------------------------
385   FindOutliers(runs,
386                nRuns,
387                useMC,
388                icut,
389                nSigmaCut,
390                fout,
391                
392                TPCnsigMeanTrendPion,
393                TPCnsigMeanTrendKaon,
394                TPCnsigMeanTrendProton,
395                TPCnsigSigmaTrendPion,
396                TPCnsigSigmaTrendKaon,
397                TPCnsigSigmaTrendProton,
398                TOFnsigMeanTrendPion,
399                TOFnsigMeanTrendKaon,
400                TOFnsigMeanTrendProton,
401                TOFnsigSigmaTrendPion,
402                TOFnsigSigmaTrendKaon,
403                TOFnsigSigmaTrendProton,
404                
405                IntegRawYieldAll,
406                
407                EfficiencyPiPlus,
408                EfficiencyKPlus,
409                EfficiencyProton,
410                EfficiencyPiMinus,
411                EfficiencyKMinus,
412                EfficiencyAntiproton,
413                
414                MatchEffPos,
415                MatchEffNeg);
416
417   
418   // close the output file
419   fout->Close();
420   // also close the unused efficiency canvas
421   if (useMC) cEfficienciesAllRuns->Close();
422   
423   Printf("\n\n--- End of AODQAChecks() ---\n");
424 }
425
426
427
428
429
430
431