TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / HighLevelQA / AODQAChecks.C
CommitLineData
22d1b439 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>
22using namespace std;
23
24// global constants
25const Int_t nPart = 3;
26TString Particle[nPart] = {"Pion", "Kaon", "Proton"};
27const Int_t nCharge = 2;
28TString Charge[nCharge] = {"Pos", "Neg"};
29TString Sign[nCharge] = {"Plus", "Minus"};
30Int_t Color[nPart] = {1,2,4};
31Int_t Marker[nPart*nCharge] = {20,21,22,24,25,26};
32TString 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
40Double_t CentCutMin[4] = {0, 30, 30, 30};
41Double_t CentCutMax[4] = {5, 40, 40, 40};
42Double_t QvecCutMin[4] = {0, 0, 0, 1.5};
43Double_t QvecCutMax[4] = {100, 100, 0.4, 100};
44Double_t EtaMin[4] = {-0.8,-0.8,-0.8,-0.8};
45Double_t EtaMax[4] = {0.8, 0.8, 0.8, 0.8};
46Double_t Nsigmapid = 3.;
47UInt_t trkbit = 1024;
48
49// fixed Pt values chosen for efficiency and matching efficiency plots (GeV)
50const Float_t FixedPtEff = 0.4;
51const 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
67void 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");
af472fff 81 gSystem->Load("libTender");
22d1b439 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