Working on the electron cut
[u/mrichter/AliRoot.git] / prod / LHC09a1 / CheckESD.C
1 #if !defined( __CINT__) || defined(__MAKECINT__)\r
2 #include <TFile.h>\r
3 #include <TError.h>\r
4 #include <TH1.h>\r
5 #include <TH2.h>\r
6 #include <TF1.h>\r
7 #include <TCanvas.h>\r
8 #include <TVector3.h>\r
9 #include <TPDGCode.h>\r
10 #include <TParticle.h>\r
11 \r
12 #include "AliRunLoader.h"\r
13 #include "AliLoader.h"\r
14 #include "AliESDEvent.h"\r
15 #include "AliRun.h"\r
16 #include "AliStack.h"\r
17 #include "AliHeader.h"\r
18 #include "AliGenEventHeader.h"\r
19 #include "AliPID.h"\r
20 #else\r
21 // const Int_t kXiMinus = 3312;\r
22 // const Int_t kOmegaMinus = 3334;\r
23 #endif\r
24 \r
25 \r
26 TH1F* CreateHisto(const char* name, const char* title, \r
27                   Int_t nBins, Double_t xMin, Double_t xMax,\r
28                   const char* xLabel = NULL, const char* yLabel = NULL)\r
29 {\r
30 // create a histogram\r
31 \r
32   TH1F* result = new TH1F(name, title, nBins, xMin, xMax);\r
33   result->SetOption("E");\r
34   if (xLabel) result->GetXaxis()->SetTitle(xLabel);\r
35   if (yLabel) result->GetYaxis()->SetTitle(yLabel);\r
36   result->SetMarkerStyle(kFullCircle);\r
37   return result;\r
38 }\r
39 \r
40 TH1F* CreateEffHisto(TH1F* hGen, TH1F* hRec)\r
41 {\r
42 // create an efficiency histogram\r
43 \r
44   Int_t nBins = hGen->GetNbinsX();\r
45   TH1F* hEff = (TH1F*) hGen->Clone("hEff");\r
46   hEff->SetTitle("");\r
47   hEff->SetStats(kFALSE);\r
48   hEff->SetMinimum(0.);\r
49   hEff->SetMaximum(110.);\r
50   hEff->GetYaxis()->SetTitle("#epsilon [%]");\r
51 \r
52   for (Int_t iBin = 0; iBin <= nBins; iBin++) {\r
53     Double_t nGen = hGen->GetBinContent(iBin);\r
54     Double_t nRec = hRec->GetBinContent(iBin);\r
55     if (nGen > 0) {\r
56       Double_t eff = nRec/nGen;\r
57       hEff->SetBinContent(iBin, 100. * eff);\r
58       Double_t error = sqrt(eff*(1.-eff) / nGen);\r
59       if (error == 0) error = 0.0001;\r
60       hEff->SetBinError(iBin, 100. * error);\r
61     } else {\r
62       hEff->SetBinContent(iBin, -100.);\r
63       hEff->SetBinError(iBin, 0);\r
64     }\r
65   }\r
66 \r
67   return hEff;\r
68 }\r
69 \r
70 Bool_t FitHisto(TH1* histo, Double_t& res, Double_t& resError)\r
71 {\r
72 // fit a gaussian to a histogram\r
73 \r
74   static TF1* fitFunc = new TF1("fitFunc", "gaus");\r
75   fitFunc->SetLineWidth(2);\r
76   fitFunc->SetFillStyle(0);\r
77   Double_t maxFitRange = 2;\r
78 \r
79   if (histo->Integral() > 50) {\r
80     Float_t mean = histo->GetMean();\r
81     Float_t rms = histo->GetRMS();\r
82     fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms);\r
83     fitFunc->SetParameters(mean, rms);\r
84     histo->Fit(fitFunc, "QRI0");\r
85     histo->GetFunction("fitFunc")->ResetBit(1<<9);\r
86     res = TMath::Abs(fitFunc->GetParameter(2));\r
87     resError = TMath::Abs(fitFunc->GetParError(2));\r
88     return kTRUE;\r
89   }\r
90 \r
91   return kFALSE;\r
92 }\r
93 \r
94 \r
95 Bool_t CheckESD(const char* gAliceFileName = "galice.root", \r
96                 const char* esdFileName = "AliESDs.root")\r
97 {\r
98 // check the content of the ESD\r
99  \r
100   // check values\r
101   Int_t    checkNGenLow = 1;\r
102 \r
103   Double_t checkEffLow = 0.5;\r
104   Double_t checkEffSigma = 3;\r
105   Double_t checkFakeHigh = 0.5;\r
106   Double_t checkFakeSigma = 3;\r
107 \r
108   Double_t checkResPtInvHigh = 5;\r
109   Double_t checkResPtInvSigma = 3;\r
110   Double_t checkResPhiHigh = 10;\r
111   Double_t checkResPhiSigma = 3;\r
112   Double_t checkResThetaHigh = 10;\r
113   Double_t checkResThetaSigma = 3;\r
114 \r
115   Double_t checkPIDEffLow = 0.5;\r
116   Double_t checkPIDEffSigma = 3;\r
117   Double_t checkResTOFHigh = 500;\r
118   Double_t checkResTOFSigma = 3;\r
119 \r
120   Double_t checkPHOSNLow = 5;\r
121   Double_t checkPHOSEnergyLow = 0.3;\r
122   Double_t checkPHOSEnergyHigh = 1.0;\r
123   Double_t checkEMCALNLow = 50;\r
124   Double_t checkEMCALEnergyLow = 0.05;\r
125   Double_t checkEMCALEnergyHigh = 1.0;\r
126 \r
127   Double_t checkMUONNLow = 1;\r
128   Double_t checkMUONPtLow = 0.5;\r
129   Double_t checkMUONPtHigh = 10.;\r
130 \r
131   Double_t cutPtV0 = 0.3;\r
132   Double_t checkV0EffLow = 0.02;\r
133   Double_t checkV0EffSigma = 3;\r
134   Double_t cutPtCascade = 0.5;\r
135   Double_t checkCascadeEffLow = 0.01;\r
136   Double_t checkCascadeEffSigma = 3;\r
137 \r
138   // open run loader and load gAlice, kinematics and header\r
139   AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName);\r
140   if (!runLoader) {\r
141     Error("CheckESD", "getting run loader from file %s failed", \r
142             gAliceFileName);\r
143     return kFALSE;\r
144   }\r
145   runLoader->LoadgAlice();\r
146   gAlice = runLoader->GetAliRun();\r
147   if (!gAlice) {\r
148     Error("CheckESD", "no galice object found");\r
149     return kFALSE;\r
150   }\r
151   runLoader->LoadKinematics();\r
152   runLoader->LoadHeader();\r
153 \r
154   // open the ESD file\r
155   TFile* esdFile = TFile::Open(esdFileName);\r
156   if (!esdFile || !esdFile->IsOpen()) {\r
157     Error("CheckESD", "opening ESD file %s failed", esdFileName);\r
158     return kFALSE;\r
159   }\r
160   AliESDEvent* esd = new AliESDEvent;\r
161   TTree* tree = (TTree*) esdFile->Get("esdTree");\r
162   if (!tree) {\r
163     Error("CheckESD", "no ESD tree found");\r
164     return kFALSE;\r
165   }\r
166   esd->ReadFromTree(tree);\r
167 \r
168   // efficiency and resolution histograms\r
169   Int_t nBinsPt = 15;\r
170   Float_t minPt = 0.1;\r
171   Float_t maxPt = 3.1;\r
172   TH1F* hGen = CreateHisto("hGen", "generated tracks", \r
173                            nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N");\r
174   TH1F* hRec = CreateHisto("hRec", "reconstructed tracks", \r
175                            nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N");\r
176   Int_t nGen = 0;\r
177   Int_t nRec = 0;\r
178   Int_t nFake = 0;\r
179 \r
180   TH1F* hResPtInv = CreateHisto("hResPtInv", "", 100, -10, 10, \r
181            "(p_{t,rec}^{-1}-p_{t,sim}^{-1}) / p_{t,sim}^{-1} [%]", "N");\r
182   TH1F* hResPhi = CreateHisto("hResPhi", "", 100, -20, 20, \r
183                               "#phi_{rec}-#phi_{sim} [mrad]", "N");\r
184   TH1F* hResTheta = CreateHisto("hResTheta", "", 100, -20, 20, \r
185                                 "#theta_{rec}-#theta_{sim} [mrad]", "N");\r
186 \r
187   // PID\r
188   Int_t partCode[] = \r
189     {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};\r
190   const char* partName[] = \r
191     {"electron", "muon", "pion", "kaon", "proton", "other"};\r
192   Double_t partFrac[] = \r
193     {0.01, 0.01, 0.85, 0.10, 0.05};\r
194   Int_t identified[6][5];\r
195   for (Int_t iGen = 0; iGen < 6; iGen++) {\r
196     for (Int_t iRec = 0; iRec < 5; iRec++) {\r
197       identified[iGen][iRec] = 0;\r
198     }\r
199   }\r
200   Int_t nIdentified = 0;\r
201 \r
202   // dE/dx and TOF\r
203   TH2F* hDEdxRight = new TH2F("hDEdxRight", "", 300, 0, 3, 100, 0, 400);\r
204   hDEdxRight->SetStats(kFALSE);\r
205   hDEdxRight->GetXaxis()->SetTitle("p [GeV/c]");\r
206   hDEdxRight->GetYaxis()->SetTitle("dE/dx_{TPC}");\r
207   hDEdxRight->SetMarkerStyle(kFullCircle);\r
208   hDEdxRight->SetMarkerSize(0.4);\r
209   TH2F* hDEdxWrong = new TH2F("hDEdxWrong", "", 300, 0, 3, 100, 0, 400);\r
210   hDEdxWrong->SetStats(kFALSE);\r
211   hDEdxWrong->GetXaxis()->SetTitle("p [GeV/c]");\r
212   hDEdxWrong->GetYaxis()->SetTitle("dE/dx_{TPC}");\r
213   hDEdxWrong->SetMarkerStyle(kFullCircle);\r
214   hDEdxWrong->SetMarkerSize(0.4);\r
215   hDEdxWrong->SetMarkerColor(kRed);\r
216   TH1F* hResTOFRight = CreateHisto("hResTOFRight", "", 100, -1000, 1000, \r
217                                    "t_{TOF}-t_{track} [ps]", "N");\r
218   TH1F* hResTOFWrong = CreateHisto("hResTOFWrong", "", 100, -1000, 1000, \r
219                                    "t_{TOF}-t_{track} [ps]", "N");\r
220   hResTOFWrong->SetLineColor(kRed);\r
221 \r
222   // calorimeters\r
223   TH1F* hEPHOS = CreateHisto("hEPHOS", "PHOS", 100, 0, 5, "E [GeV]", "N");\r
224   TH1F* hEEMCAL = CreateHisto("hEEMCAL", "EMCAL", 100, 0, 50, "E [GeV]", "N");\r
225 \r
226   // muons\r
227   TH1F* hPtMUON = CreateHisto("hPtMUON", "MUON", 100, 0, 20, \r
228                               "p_{t} [GeV/c]", "N");\r
229 \r
230   // V0s and cascades\r
231   TH1F* hMassK0 = CreateHisto("hMassK0", "K^{0}", 100, 0.4, 0.6, \r
232                               "M(#pi^{+}#pi^{-}) [GeV/c^{2}]", "N");\r
233   TH1F* hMassLambda = CreateHisto("hMassLambda", "#Lambda", 100, 1.0, 1.2, \r
234                                   "M(p#pi^{-}) [GeV/c^{2}]", "N");\r
235   TH1F* hMassLambdaBar = CreateHisto("hMassLambdaBar", "#bar{#Lambda}", \r
236                                      100, 1.0, 1.2, \r
237                                      "M(#bar{p}#pi^{+}) [GeV/c^{2}]", "N");\r
238   Int_t nGenV0s = 0;\r
239   Int_t nRecV0s = 0;\r
240   TH1F* hMassXi = CreateHisto("hMassXi", "#Xi", 100, 1.2, 1.5, \r
241                               "M(#Lambda#pi) [GeV/c^{2}]", "N");\r
242   TH1F* hMassOmega = CreateHisto("hMassOmega", "#Omega", 100, 1.5, 1.8, \r
243                                  "M(#LambdaK) [GeV/c^{2}]", "N");\r
244   Int_t nGenCascades = 0;\r
245   Int_t nRecCascades = 0;\r
246 \r
247   // loop over events\r
248   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {\r
249     runLoader->GetEvent(iEvent);\r
250 \r
251     // select simulated primary particles, V0s and cascades\r
252     AliStack* stack = gAlice->Stack();\r
253     Int_t nParticles = stack->GetNtrack();\r
254     TArrayF vertex(3);\r
255     runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);\r
256     TObjArray selParticles;\r
257     TObjArray selV0s;\r
258     TObjArray selCascades;\r
259     for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {\r
260       TParticle* particle = stack->Particle(iParticle);\r
261       if (!particle) continue;\r
262       if (particle->Pt() < 0.001) continue;\r
263       if (TMath::Abs(particle->Eta()) > 0.9) continue;\r
264       TVector3 dVertex(particle->Vx() - vertex[0], \r
265                        particle->Vy() - vertex[1],\r
266                        particle->Vz() - vertex[2]);\r
267       if (dVertex.Mag() > 0.0001) continue;\r
268 \r
269       switch (TMath::Abs(particle->GetPdgCode())) {\r
270       case kElectron:\r
271       case kMuonMinus:\r
272       case kPiPlus:\r
273       case kKPlus:\r
274       case kProton: {\r
275         if (particle->Pt() > minPt) {\r
276           selParticles.Add(particle);\r
277           nGen++;\r
278           hGen->Fill(particle->Pt());\r
279         }\r
280         break;\r
281       }\r
282       case kK0Short:\r
283       case kLambda0: {\r
284         if (particle->Pt() > cutPtV0) {\r
285           nGenV0s++;\r
286           selV0s.Add(particle);\r
287         }\r
288         break;\r
289       }\r
290       case kXiMinus:\r
291       case kOmegaMinus: {\r
292         if (particle->Pt() > cutPtCascade) {\r
293           nGenCascades++;\r
294           selCascades.Add(particle);\r
295         }\r
296         break;\r
297       }\r
298       default: break;\r
299       }\r
300     }\r
301 \r
302     // get the event summary data\r
303     tree->GetEvent(iEvent);\r
304     if (!esd) {\r
305       Error("CheckESD", "no ESD object found for event %d", iEvent);\r
306       return kFALSE;\r
307     }\r
308 \r
309     // loop over tracks\r
310     for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {\r
311       AliESDtrack* track = esd->GetTrack(iTrack);\r
312 \r
313       // select tracks of selected particles\r
314       Int_t label = TMath::Abs(track->GetLabel());\r
315       if (label > stack->GetNtrack()) continue;     // background\r
316       TParticle* particle = stack->Particle(label);\r
317       if (!selParticles.Contains(particle)) continue;\r
318       if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;\r
319       if (track->GetConstrainedChi2() > 1e9) continue;\r
320       selParticles.Remove(particle);   // don't count multiple tracks\r
321 \r
322       nRec++;\r
323       hRec->Fill(particle->Pt());\r
324       if (track->GetLabel() < 0) nFake++;\r
325 \r
326       // resolutions\r
327       Double_t p[3];\r
328       track->GetConstrainedPxPyPz(p);\r
329       TVector3 pTrack(p);\r
330       hResPtInv->Fill(100. * (1./pTrack.Pt() - 1./particle->Pt()) * \r
331                       particle->Pt());\r
332       hResPhi->Fill(1000. * (pTrack.Phi() - particle->Phi()));\r
333       hResTheta->Fill(1000. * (pTrack.Theta() - particle->Theta()));\r
334 \r
335       // PID\r
336       if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue;\r
337       Int_t iGen = 5;\r
338       for (Int_t i = 0; i < 5; i++) {\r
339         if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i;\r
340       }\r
341       Double_t probability[5];\r
342       track->GetESDpid(probability);\r
343       Double_t pMax = 0;\r
344       Int_t iRec = 0;\r
345       for (Int_t i = 0; i < 5; i++) {\r
346         probability[i] *= partFrac[i];\r
347         if (probability[i] > pMax) {\r
348           pMax = probability[i];\r
349           iRec = i;\r
350         }\r
351       }\r
352       identified[iGen][iRec]++;\r
353       if (iGen == iRec) nIdentified++;\r
354 \r
355       // dE/dx and TOF\r
356       Double_t time[5];\r
357       track->GetIntegratedTimes(time);\r
358       if (iGen == iRec) {\r
359         hDEdxRight->Fill(pTrack.Mag(), track->GetTPCsignal());\r
360         if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {\r
361           hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]);\r
362         }\r
363       } else {\r
364         hDEdxWrong->Fill(pTrack.Mag(), track->GetTPCsignal());\r
365         if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {\r
366           hResTOFWrong->Fill(track->GetTOFsignal() - time[iRec]);\r
367         }\r
368       }\r
369     }\r
370 \r
371     // loop over muon tracks\r
372     {\r
373     for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {\r
374       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);\r
375       Double_t ptInv = TMath::Abs(muonTrack->GetInverseBendingMomentum());\r
376       if (ptInv > 0.001) {\r
377         hPtMUON->Fill(1./ptInv);\r
378       }\r
379     }\r
380     }\r
381 \r
382     // loop over V0s\r
383     for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {\r
384       AliESDv0* v0 = esd->GetV0(iV0);\r
385       if (v0->GetOnFlyStatus()) continue;\r
386       v0->ChangeMassHypothesis(kK0Short);\r
387       hMassK0->Fill(v0->GetEffMass());\r
388       v0->ChangeMassHypothesis(kLambda0);\r
389       hMassLambda->Fill(v0->GetEffMass());\r
390       v0->ChangeMassHypothesis(kLambda0Bar);\r
391       hMassLambdaBar->Fill(v0->GetEffMass());\r
392 \r
393       Int_t negLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel());\r
394       if (negLabel > stack->GetNtrack()) continue;     // background\r
395       Int_t negMother = stack->Particle(negLabel)->GetMother(0);\r
396       if (negMother < 0) continue;\r
397       Int_t posLabel = TMath::Abs(esd->GetTrack(v0->GetPindex())->GetLabel());\r
398       if (posLabel > stack->GetNtrack()) continue;     // background\r
399       Int_t posMother = stack->Particle(posLabel)->GetMother(0);\r
400       if (negMother != posMother) continue;\r
401       TParticle* particle = stack->Particle(negMother);\r
402       if (!selV0s.Contains(particle)) continue;\r
403       selV0s.Remove(particle);\r
404       nRecV0s++;\r
405     }\r
406 \r
407     // loop over Cascades\r
408     for (Int_t iCascade = 0; iCascade < esd->GetNumberOfCascades(); \r
409          iCascade++) {\r
410       AliESDcascade* cascade = esd->GetCascade(iCascade);\r
411       Double_t v0q;\r
412       cascade->ChangeMassHypothesis(v0q,kXiMinus);\r
413       hMassXi->Fill(cascade->GetEffMass());\r
414       cascade->ChangeMassHypothesis(v0q,kOmegaMinus);\r
415       hMassOmega->Fill(cascade->GetEffMass());\r
416 \r
417       Int_t negLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex())\r
418                                   ->GetLabel());\r
419       if (negLabel > stack->GetNtrack()) continue;     // background\r
420       Int_t negMother = stack->Particle(negLabel)->GetMother(0);\r
421       if (negMother < 0) continue;\r
422       Int_t posLabel = TMath::Abs(esd->GetTrack(cascade->GetPindex())\r
423                                   ->GetLabel());\r
424       if (posLabel > stack->GetNtrack()) continue;     // background\r
425       Int_t posMother = stack->Particle(posLabel)->GetMother(0);\r
426       if (negMother != posMother) continue;\r
427       Int_t v0Mother = stack->Particle(negMother)->GetMother(0);\r
428       if (v0Mother < 0) continue;\r
429       Int_t bacLabel = TMath::Abs(esd->GetTrack(cascade->GetBindex())\r
430                                   ->GetLabel());\r
431       if (bacLabel > stack->GetNtrack()) continue;     // background\r
432       Int_t bacMother = stack->Particle(bacLabel)->GetMother(0);\r
433       if (v0Mother != bacMother) continue;\r
434       TParticle* particle = stack->Particle(v0Mother);\r
435       if (!selCascades.Contains(particle)) continue;\r
436       selCascades.Remove(particle);\r
437       nRecCascades++;\r
438     }\r
439 \r
440     // loop over the calorimeters clusters\r
441     {\r
442     Int_t nCaloCluster  = esd->GetNumberOfCaloClusters();\r
443     cout<<"CaloClusters "<<nCaloCluster<<endl;\r
444     for (Int_t iCluster=0; iCluster<nCaloCluster; iCluster++){\r
445       if(esd->GetCaloCluster(iCluster)->IsPHOS()) \r
446         hEPHOS->Fill(esd->GetCaloCluster(iCluster)->E());\r
447       if(esd->GetCaloCluster(iCluster)->GetClusterType()==AliESDCaloCluster::kEMCALClusterv1) {\r
448         hEEMCAL->Fill(esd->GetCaloCluster(iCluster)->E());\r
449         //cout<<esd->GetCaloCluster(iCluster)->E()<<endl;\r
450       }\r
451     }\r
452     }\r
453   }\r
454 \r
455   // perform checks\r
456   if (nGen < checkNGenLow) {\r
457     Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen));\r
458   }\r
459 \r
460   TH1F* hEff = CreateEffHisto(hGen, hRec);\r
461 \r
462   Info("CheckESD", "%d out of %d tracks reconstructed including %d "\r
463          "fake tracks", nRec, nGen, nFake);\r
464   if (nGen > 0) {\r
465     // efficiency\r
466     Double_t eff = nRec*1./nGen;\r
467     Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGen);\r
468     Double_t fake = nFake*1./nGen;\r
469     Double_t fakeError = TMath::Sqrt(fake*(1.-fake) / nGen);\r
470     Info("CheckESD", "eff = (%.1f +- %.1f) %%  fake = (%.1f +- %.1f) %%",\r
471          100.*eff, 100.*effError, 100.*fake, 100.*fakeError);\r
472 \r
473     if (eff < checkEffLow - checkEffSigma*effError) {\r
474       Warning("CheckESD", "low efficiency: (%.1f +- %.1f) %%", \r
475               100.*eff, 100.*effError);\r
476     }\r
477     if (fake > checkFakeHigh + checkFakeSigma*fakeError) {\r
478       Warning("CheckESD", "high fake: (%.1f +- %.1f) %%", \r
479               100.*fake, 100.*fakeError);\r
480     }\r
481 \r
482     // resolutions\r
483     Double_t res, resError;\r
484     if (FitHisto(hResPtInv, res, resError)) {\r
485       Info("CheckESD", "relative inverse pt resolution = (%.1f +- %.1f) %%",\r
486            res, resError);\r
487       if (res > checkResPtInvHigh + checkResPtInvSigma*resError) {\r
488         Warning("CheckESD", "bad pt resolution: (%.1f +- %.1f) %%", \r
489                 res, resError);\r
490       }\r
491     }\r
492 \r
493     if (FitHisto(hResPhi, res, resError)) {\r
494       Info("CheckESD", "phi resolution = (%.1f +- %.1f) mrad", res, resError);\r
495       if (res > checkResPhiHigh + checkResPhiSigma*resError) {\r
496         Warning("CheckESD", "bad phi resolution: (%.1f +- %.1f) mrad", \r
497                 res, resError);\r
498       }\r
499     }\r
500 \r
501     if (FitHisto(hResTheta, res, resError)) {\r
502       Info("CheckESD", "theta resolution = (%.1f +- %.1f) mrad", \r
503            res, resError);\r
504       if (res > checkResThetaHigh + checkResThetaSigma*resError) {\r
505         Warning("CheckESD", "bad theta resolution: (%.1f +- %.1f) mrad", \r
506                 res, resError);\r
507       }\r
508     }\r
509 \r
510     // PID\r
511     if (nRec > 0) {\r
512       Double_t eff = nIdentified*1./nRec;\r
513       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nRec);\r
514       Info("CheckESD", "PID eff = (%.1f +- %.1f) %%", \r
515            100.*eff, 100.*effError);\r
516       if (eff < checkPIDEffLow - checkPIDEffSigma*effError) {\r
517         Warning("CheckESD", "low PID efficiency: (%.1f +- %.1f) %%", \r
518                 100.*eff, 100.*effError);\r
519       }\r
520     }\r
521 \r
522     printf("%9s:", "gen\\rec");\r
523     for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {\r
524       printf("%9s", partName[iRec]);\r
525     }\r
526     printf("\n");\r
527     for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) {\r
528       printf("%9s:", partName[iGen]);\r
529       for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {\r
530         printf("%9d", identified[iGen][iRec]);\r
531       }\r
532       printf("\n");\r
533     }\r
534 \r
535     if (FitHisto(hResTOFRight, res, resError)) {\r
536       Info("CheckESD", "TOF resolution = (%.1f +- %.1f) ps", res, resError);\r
537       if (res > checkResTOFHigh + checkResTOFSigma*resError) {\r
538         Warning("CheckESD", "bad TOF resolution: (%.1f +- %.1f) ps", \r
539                 res, resError);\r
540       }\r
541     }\r
542 \r
543     // calorimeters\r
544     if (hEPHOS->Integral() < checkPHOSNLow) {\r
545       Warning("CheckESD", "low number of PHOS particles: %d", \r
546               Int_t(hEPHOS->Integral()));\r
547     } else {\r
548       Double_t mean = hEPHOS->GetMean();\r
549       if (mean < checkPHOSEnergyLow) {\r
550         Warning("CheckESD", "low mean PHOS energy: %.1f GeV", mean);\r
551       } else if (mean > checkPHOSEnergyHigh) {\r
552         Warning("CheckESD", "high mean PHOS energy: %.1f GeV", mean);\r
553       }\r
554     }\r
555 \r
556     if (hEEMCAL->Integral() < checkEMCALNLow) {\r
557       Warning("CheckESD", "low number of EMCAL particles: %d", \r
558               Int_t(hEEMCAL->Integral()));\r
559     } else {\r
560       Double_t mean = hEEMCAL->GetMean();\r
561       if (mean < checkEMCALEnergyLow) {\r
562         Warning("CheckESD", "low mean EMCAL energy: %.1f GeV", mean);\r
563       } else if (mean > checkEMCALEnergyHigh) {\r
564         Warning("CheckESD", "high mean EMCAL energy: %.1f GeV", mean);\r
565       }\r
566     }\r
567 \r
568     // muons\r
569     if (hPtMUON->Integral() < checkMUONNLow) {\r
570       Warning("CheckESD", "low number of MUON particles: %d", \r
571               Int_t(hPtMUON->Integral()));\r
572     } else {\r
573       Double_t mean = hPtMUON->GetMean();\r
574       if (mean < checkMUONPtLow) {\r
575         Warning("CheckESD", "low mean MUON pt: %.1f GeV/c", mean);\r
576       } else if (mean > checkMUONPtHigh) {\r
577         Warning("CheckESD", "high mean MUON pt: %.1f GeV/c", mean);\r
578       }\r
579     }\r
580 \r
581     // V0s\r
582     if (nGenV0s > 0) {\r
583       Double_t eff = nRecV0s*1./nGenV0s;\r
584       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenV0s);\r
585       if (effError == 0) effError = checkV0EffLow / TMath::Sqrt(1.*nGenV0s);\r
586       Info("CheckESD", "V0 eff = (%.1f +- %.1f) %%", \r
587            100.*eff, 100.*effError);\r
588       if (eff < checkV0EffLow - checkV0EffSigma*effError) {\r
589         Warning("CheckESD", "low V0 efficiency: (%.1f +- %.1f) %%", \r
590                 100.*eff, 100.*effError);\r
591       }\r
592     }\r
593 \r
594     // Cascades\r
595     if (nGenCascades > 0) {\r
596       Double_t eff = nRecCascades*1./nGenCascades;\r
597       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenCascades);\r
598       if (effError == 0) effError = checkV0EffLow / \r
599                            TMath::Sqrt(1.*nGenCascades);\r
600       Info("CheckESD", "Cascade eff = (%.1f +- %.1f) %%", \r
601            100.*eff, 100.*effError);\r
602       if (eff < checkCascadeEffLow - checkCascadeEffSigma*effError) {\r
603         Warning("CheckESD", "low Cascade efficiency: (%.1f +- %.1f) %%", \r
604                 100.*eff, 100.*effError);\r
605       }\r
606     }\r
607   }\r
608 \r
609   // draw the histograms if not in batch mode\r
610   if (!gROOT->IsBatch()) {\r
611     new TCanvas;\r
612     hEff->DrawCopy();\r
613     new TCanvas;\r
614     hResPtInv->DrawCopy("E");\r
615     new TCanvas;\r
616     hResPhi->DrawCopy("E");\r
617     new TCanvas;\r
618     hResTheta->DrawCopy("E");\r
619     new TCanvas;\r
620     hDEdxRight->DrawCopy();\r
621     hDEdxWrong->DrawCopy("SAME");\r
622     new TCanvas;\r
623     hResTOFRight->DrawCopy("E");\r
624     hResTOFWrong->DrawCopy("SAME");\r
625     new TCanvas;\r
626     hEPHOS->DrawCopy("E");\r
627     new TCanvas;\r
628     hEEMCAL->DrawCopy("E");\r
629     new TCanvas;\r
630     hPtMUON->DrawCopy("E");\r
631     new TCanvas;\r
632     hMassK0->DrawCopy("E");\r
633     new TCanvas;\r
634     hMassLambda->DrawCopy("E");\r
635     new TCanvas;\r
636     hMassLambdaBar->DrawCopy("E");\r
637     new TCanvas;\r
638     hMassXi->DrawCopy("E");\r
639     new TCanvas;\r
640     hMassOmega->DrawCopy("E");\r
641   }\r
642 \r
643   // write the output histograms to a file\r
644   TFile* outputFile = TFile::Open("check.root", "recreate");\r
645   if (!outputFile || !outputFile->IsOpen()) {\r
646     Error("CheckESD", "opening output file check.root failed");\r
647     return kFALSE;\r
648   }\r
649   hEff->Write();\r
650   hResPtInv->Write();\r
651   hResPhi->Write();\r
652   hResTheta->Write();\r
653   hDEdxRight->Write();\r
654   hDEdxWrong->Write();\r
655   hResTOFRight->Write();\r
656   hResTOFWrong->Write();\r
657   hEPHOS->Write();\r
658   hEEMCAL->Write();\r
659   hPtMUON->Write();\r
660   hMassK0->Write();\r
661   hMassLambda->Write();\r
662   hMassLambdaBar->Write();\r
663   hMassXi->Write();\r
664   hMassOmega->Write();\r
665   outputFile->Close();\r
666   delete outputFile;\r
667 \r
668   // clean up\r
669   delete hGen;\r
670   delete hRec;\r
671   delete hEff;\r
672   delete hResPtInv;\r
673   delete hResPhi;\r
674   delete hResTheta;\r
675   delete hDEdxRight;\r
676   delete hDEdxWrong;\r
677   delete hResTOFRight;\r
678   delete hResTOFWrong;\r
679   delete hEPHOS;\r
680   delete hEEMCAL;\r
681   delete hPtMUON;\r
682   delete hMassK0;\r
683   delete hMassLambda;\r
684   delete hMassLambdaBar;\r
685   delete hMassXi;\r
686   delete hMassOmega;\r
687 \r
688   delete esd;\r
689   esdFile->Close();\r
690   delete esdFile;\r
691 \r
692   runLoader->UnloadHeader();\r
693   runLoader->UnloadKinematics();\r
694   delete runLoader;\r
695 \r
696   // result of check\r
697   Info("CheckESD", "check of ESD was successfull");\r
698   return kTRUE;\r
699 }\r