16eec53a663ded7501c7b4b4ad6650b64ab46514
[u/mrichter/AliRoot.git] / PWG1 / comparison / AliAnalysisTaskCheckESD.cxx
1 #include "TChain.h"
2 #include "TROOT.h"
3 #include "TFile.h"
4 #include "TError.h"
5 #include "TH1.h"
6 #include "TH2.h"
7 #include "TF1.h"
8 #include "TArrayI.h"
9 #include "TCanvas.h"
10 #include "TVector3.h"
11 #include "TPDGCode.h"
12 #include "TParticle.h"
13 #include "AliESDEvent.h"
14 #include "AliESDv0.h"
15 #include "AliESDcascade.h"
16 #include "AliESDMuonTrack.h"
17 #include "AliESDCaloCluster.h"
18 #include "AliRun.h"
19 #include "AliHeader.h"
20 #include "AliMCEvent.h"
21 #include "AliGenEventHeader.h"
22 #include "AliPID.h"
23
24 #include "AliAnalysisTaskCheckESD.h"
25
26 ClassImp(AliAnalysisTaskCheckESD)
27
28 AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD():
29 AliAnalysisTaskSE("AliAnalysisTaskCheckESD"),
30   fListOfHistos(0),
31   hGen(0),
32   hRec(0),
33   hResPtInv(0),
34   hResPhi(0),
35   hResTheta(0),
36   hDEdxRight(0),
37   hDEdxWrong(0),
38   hResTOFRight(0),
39   hResTOFWrong(0),
40   hEPHOS(0),
41   hEEMCAL(0),
42   hPtMUON(0),
43   hMassK0(0),
44   hMassLambda(0),
45   hMassLambdaBar(0),
46   hMassXi(0),
47   hMassOmega(0),
48   hScalars(0),
49   hArrayHist(0)
50 {
51   // Default constructor
52   // Define input and output slots here
53   // Input slot #0 works with a TChain
54   DefineInput(0, TChain::Class());
55   // Output slot #1 TList
56   DefineOutput(1, TList::Class());
57 }
58
59 AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD(const char* name):
60 AliAnalysisTaskSE(name),
61   fListOfHistos(0),
62   hGen(0),
63   hRec(0),
64   hResPtInv(0),
65   hResPhi(0),
66   hResTheta(0),
67   hDEdxRight(0),
68   hDEdxWrong(0),
69   hResTOFRight(0),
70   hResTOFWrong(0),
71   hEPHOS(0),
72   hEEMCAL(0),
73   hPtMUON(0),
74   hMassK0(0),
75   hMassLambda(0),
76   hMassLambdaBar(0),
77   hMassXi(0),
78   hMassOmega(0),
79   hScalars(0),
80   hArrayHist(0)
81 {
82   // Constructor
83   AliInfo("Constructor AliAnalysisTaskCheckESD");
84   // Define input and output slots here
85   // Input slot #0 works with a TChain
86   DefineInput(0, TChain::Class());
87   // Output slot #1 TList
88   DefineOutput(1, TList::Class());
89 }
90
91 TH1F * AliAnalysisTaskCheckESD::CreateHisto(const char* name, const char* title,Int_t nBins, 
92                                             Double_t xMin, Double_t xMax,
93                                             const char* xLabel, const char* yLabel)
94 {
95   // create a histogram
96   TH1F* result = new TH1F(name, title, nBins, xMin, xMax);
97   result->SetOption("E");
98   if (xLabel) result->GetXaxis()->SetTitle(xLabel);
99   if (yLabel) result->GetYaxis()->SetTitle(yLabel);
100   result->SetMarkerStyle(kFullCircle);
101   return result;
102 }
103
104 TH1F * AliAnalysisTaskCheckESD::CreateEffHisto(TH1F* hGen, TH1F* hRec)
105 {
106   // create an efficiency histogram
107   Int_t nBins = hGen->GetNbinsX();
108   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
109   hEff->SetTitle("");
110   hEff->SetStats(kFALSE);
111   hEff->SetMinimum(0.);
112   hEff->SetMaximum(110.);
113   hEff->GetYaxis()->SetTitle("#epsilon [%]");
114   
115   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
116     Double_t nGen_eff = hGen->GetBinContent(iBin);
117     Double_t nRec_eff = hRec->GetBinContent(iBin);
118     if (nGen_eff > 0) {
119       Double_t eff = nRec_eff/nGen_eff;
120       hEff->SetBinContent(iBin, 100. * eff);
121       Double_t error = sqrt(eff*(1.-eff) / nGen_eff);
122       if (error == 0) error = 0.0001;
123       hEff->SetBinError(iBin, 100. * error);                    
124     }
125     else {
126       hEff->SetBinContent(iBin, -100.);
127       hEff->SetBinError(iBin, 0);
128     }
129   }
130   return hEff;
131 }
132
133
134 Bool_t AliAnalysisTaskCheckESD::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
135 {
136   // fit a gaussian to a histogram
137   static TF1* fitFunc = new TF1("fitFunc", "gaus");
138   fitFunc->SetLineWidth(2);
139   fitFunc->SetFillStyle(0);
140   Double_t maxFitRange = 2;
141   
142   if (histo->Integral() > 50) {
143     Float_t mean = histo->GetMean();
144     Float_t rms = histo->GetRMS();
145     fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms);
146     fitFunc->SetParameters(mean, rms);
147     histo->Fit(fitFunc, "QRI0");
148     histo->GetFunction("fitFunc")->ResetBit(1<<9);
149     res = TMath::Abs(fitFunc->GetParameter(2));
150     resError = TMath::Abs(fitFunc->GetParError(2));
151     return kTRUE;
152   }
153   return kFALSE;
154 }
155
156 void AliAnalysisTaskCheckESD::UserCreateOutputObjects()
157 {
158   // Create histograms
159   AliInfo("AliAnalysisTaskCheckESD::UserCreateOutputObjects");
160   // Create output container
161   fListOfHistos = new TList();
162   
163   // efficiency and resolution histog
164   Int_t nBinsPt = 15;
165   Float_t minPt = 0.1;
166   Float_t maxPt = 3.1;
167   
168   hGen = CreateHisto("hGen", "generated tracks", nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N");
169   hRec = CreateHisto("hRec", "reconstructed tracks", nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N");
170   hResPtInv = CreateHisto("hResPtInv", "", 100, -10, 10, "(p_{t,rec}^{-1}-p_{t,sim}^{-1}) / p_{t,sim}^{-1} [%]", "N");
171   hResPhi = CreateHisto("hResPhi", "", 100, -20, 20, "#phi_{rec}-#phi_{sim} [mrad]", "N");
172   hResTheta = CreateHisto("hResTheta", "", 100, -20, 20, "#theta_{rec}-#theta_{sim} [mrad]", "N");
173   
174   // dE/dx and TOF
175   hDEdxRight = new TH2F("hDEdxRight", "", 300, 0, 3, 100, 0, 400);
176   hDEdxRight->SetStats(kFALSE);
177   hDEdxRight->GetXaxis()->SetTitle("p [GeV/c]");
178   hDEdxRight->GetYaxis()->SetTitle("dE/dx_{TPC}");
179   hDEdxRight->SetMarkerStyle(kFullCircle);
180   hDEdxRight->SetMarkerSize(0.4);
181   hDEdxWrong = new TH2F("hDEdxWrong", "", 300, 0, 3, 100, 0, 400);
182   hDEdxWrong->SetStats(kFALSE);
183   hDEdxWrong->GetXaxis()->SetTitle("p [GeV/c]");
184   hDEdxWrong->GetYaxis()->SetTitle("dE/dx_{TPC}");
185   hDEdxWrong->SetMarkerStyle(kFullCircle);
186   hDEdxWrong->SetMarkerSize(0.4);
187   hDEdxWrong->SetMarkerColor(kRed);
188   
189   hResTOFRight = CreateHisto("hResTOFRight", "", 100, -1000, 1000, "t_{TOF}-t_{track} [ps]", "N");
190   hResTOFWrong = CreateHisto("hResTOFWrong", "", 100, -1000, 1000, "t_{TOF}-t_{track} [ps]", "N");
191   hResTOFWrong->SetLineColor(kRed);
192   
193   // calorimeters
194   hEPHOS = CreateHisto("hEPHOS", "PHOS", 100, 0, 50, "E [GeV]", "N");
195   hEEMCAL = CreateHisto("hEEMCAL", "EMCAL", 100, 0, 50, "E [GeV]", "N");
196   
197   // muons
198   hPtMUON = CreateHisto("hPtMUON", "MUON", 100, 0, 20, "p_{t} [GeV/c]", "N");
199   
200   // V0s and cascades
201   hMassK0 = CreateHisto("hMassK0", "K^{0}", 100, 0.4, 0.6, "M(#pi^{+}#pi^{-}) [GeV/c^{2}]", "N");
202   hMassLambda = CreateHisto("hMassLambda", "#Lambda", 100, 1.0, 1.2, "M(p#pi^{-}) [GeV/c^{2}]", "N");
203   
204   hMassLambdaBar = CreateHisto("hMassLambdaBar", "#bar{#Lambda}", 100, 1.0, 1.2, "M(#bar{p}#pi^{+}) [GeV/c^{2}]", "N");
205   hMassXi = CreateHisto("hMassXi", "#Xi", 100, 1.2, 1.5, "M(#Lambda#pi) [GeV/c^{2}]", "N");
206   hMassOmega = CreateHisto("hMassOmega", "#Omega", 100, 1.5, 1.8, "M(#LambdaK) [GeV/c^{2}]", "N");
207   hScalars = new TH1F("hScalars","Container of scalars",8,0,8);
208   hArrayHist = new TH1F("hArrayHist","Container for Array",
209                         (AliPID::kSPECIES+1)*AliPID::kSPECIES,0,(AliPID::kSPECIES+1)*AliPID::kSPECIES);
210   
211   fListOfHistos->Add(hGen);
212   fListOfHistos->Add(hRec);
213   fListOfHistos->Add(hResPtInv);
214   fListOfHistos->Add(hResPhi);
215   fListOfHistos->Add(hResTheta);
216   fListOfHistos->Add(hDEdxRight);
217   fListOfHistos->Add(hDEdxWrong);
218   fListOfHistos->Add(hResTOFRight);
219   fListOfHistos->Add(hResTOFWrong);
220   fListOfHistos->Add(hEPHOS);
221   fListOfHistos->Add(hEEMCAL);
222   fListOfHistos->Add(hPtMUON);
223   fListOfHistos->Add(hMassK0);
224   fListOfHistos->Add(hMassLambda);
225   fListOfHistos->Add(hMassLambdaBar);
226   fListOfHistos->Add(hMassXi);
227   fListOfHistos->Add(hMassOmega);
228   fListOfHistos->Add(hScalars);
229   fListOfHistos->Add(hArrayHist);
230 }
231
232 void AliAnalysisTaskCheckESD::UserExec(Option_t *option)
233 {
234   const Int_t kXiMinus = 3312;
235   const Int_t kOmegaMinus = 3334;
236   // check the content of the ESD
237   Double_t cutPtV0 = 0.3;
238   Double_t cutPtCascade = 0.5;
239   Float_t minPt = 0.1;
240   
241   // PID
242   Int_t partCode[AliPID::kSPECIES] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
243   Double_t partFrac[AliPID::kSPECIES] = {0.01, 0.01, 0.85, 0.10, 0.05};
244   
245   AliMCEvent* mcEvent = MCEvent();
246   if (!mcEvent) {
247     Printf("ERROR: Could not retrieve MC event");
248     return;
249   }     
250         
251   //Primary vertex needed
252   TArrayF vertex(3);  
253   mcEvent->GenEventHeader()->PrimaryVertex(vertex);
254   
255   TObjArray selParticles; //Use TClonesArray?
256   TObjArray selV0s;
257   TObjArray selCascades;
258   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again
259     AliMCParticle* track = mcEvent->GetTrack(iTracks);
260     TParticle* particle = track->Particle();
261     if (!particle) continue;
262     if (particle->Pt() < 0.001) continue;
263     if (TMath::Abs(particle->Eta()) > 0.9) continue;
264     TVector3 dVertex(particle->Vx() - vertex[0], particle->Vy() - vertex[1], particle->Vz() - vertex[2]);
265     if (dVertex.Mag() > 0.0001) continue;
266     
267     switch (TMath::Abs(particle->GetPdgCode())) {
268     case kElectron:
269     case kMuonMinus:
270     case kPiPlus:
271     case kKPlus:
272     case kProton:
273       {
274         if (particle->Pt() > minPt) {
275           selParticles.Add(particle);
276           hScalars->Fill(0);
277           hGen->Fill(particle->Pt());
278         }
279         break;
280       }
281     case kK0Short:
282     case kLambda0:
283       {
284         if (particle->Pt() > cutPtV0) {
285           hScalars->Fill(3);
286           selV0s.Add(particle);
287         }
288         break;
289       }
290     case kXiMinus:
291     case kOmegaMinus:
292       {
293         if (particle->Pt() > cutPtCascade) {
294           hScalars->Fill(6);
295           selCascades.Add(particle);
296         }
297         break;
298       }
299     default: break;
300     }           
301   }
302         
303   // ESD information
304   AliVEvent* event = InputEvent();
305   if (!event) {
306     Printf("ERROR: Could not retrieve event");
307     return;
308   }
309
310   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
311   // loop over tracks
312   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
313     AliESDtrack* track = esd->GetTrack(iTrack);
314
315     // select tracks of selected particles
316     Int_t label = TMath::Abs(track->GetLabel());
317     if (label > mcEvent->GetNumberOfTracks()) continue;     // background
318     AliMCParticle* mctrack = mcEvent->GetTrack(label);
319     TParticle* particle = mctrack->Particle();
320     if (!selParticles.Contains(particle)) continue;
321     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
322     if (track->GetConstrainedChi2() > 1e9) continue;
323     selParticles.Remove(particle);   // don't count multiple tracks
324     
325     hScalars->Fill(1);
326     hRec->Fill(particle->Pt());
327     if (track->GetLabel() < 0) { 
328       hScalars->Fill(2);
329     }
330
331     // resolutions
332     hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) *particle->Pt());
333     hResPhi->Fill(1000. * (track->Phi() - particle->Phi()));
334     hResTheta->Fill(1000. * (track->Theta() - particle->Theta()));
335     
336     // PID
337     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
338     Int_t iGen = 5;
339     
340     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
341       if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i;
342     }
343     Double_t probability[AliPID::kSPECIES];
344     track->GetESDpid(probability);
345     Double_t pMax = 0;
346     Int_t iRec = 0;
347     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
348       probability[i] *= partFrac[i];
349       if (probability[i] > pMax) {
350         pMax = probability[i];
351         iRec = i;
352       }
353                         
354     }
355     hArrayHist->Fill(AliPID::kSPECIES*iGen + iRec);
356     if (iGen == iRec) {
357       hScalars->Fill(5);
358     }
359
360     // dE/dx and TOF
361     Double_t time[AliPID::kSPECIES];
362     track->GetIntegratedTimes(time);
363     if (iGen == iRec) {
364       hDEdxRight->Fill(particle->P(), track->GetTPCsignal());
365       if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {
366         hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]);
367       }
368     }
369     else {
370       hDEdxWrong->Fill(particle->P(), track->GetTPCsignal());
371       if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {
372         hResTOFWrong->Fill(track->GetTOFsignal() - time[iRec]);
373       }
374     }
375   }
376
377   // loop over muon tracks
378   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
379     AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
380     Double_t ptInv = TMath::Abs(muonTrack->GetInverseBendingMomentum());
381     if (ptInv > 0.001) {
382       hPtMUON->Fill(1./ptInv);
383     }
384   }
385
386   // loop over V0s
387   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
388     AliESDv0* v0 = esd->GetV0(iV0);
389     if (v0->GetOnFlyStatus()) continue;
390     v0->ChangeMassHypothesis(kK0Short);
391     hMassK0->Fill(v0->GetEffMass());
392     v0->ChangeMassHypothesis(kLambda0);
393     hMassLambda->Fill(v0->GetEffMass());
394     v0->ChangeMassHypothesis(kLambda0Bar);
395     hMassLambdaBar->Fill(v0->GetEffMass());
396     
397     Int_t negLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel());
398     if (negLabel > mcEvent->GetNumberOfTracks()) continue;     // background
399     AliMCParticle* track = mcEvent->GetTrack(negLabel);
400     TParticle* particle = track->Particle();
401     Int_t negMother = particle->GetMother(0);
402     if (negMother < 0) continue;
403     Int_t posLabel = TMath::Abs(esd->GetTrack(v0->GetPindex())->GetLabel());
404     if (posLabel > mcEvent->GetNumberOfTracks()) continue;     // background
405     track = mcEvent->GetTrack(posLabel);
406     particle = track->Particle();
407     Int_t posMother = particle->GetMother(0);
408     if (negMother != posMother) continue;
409     track = mcEvent->GetTrack(negMother);
410     particle = track->Particle();
411     if (!selV0s.Contains(particle)) continue;
412     selV0s.Remove(particle);
413     hScalars->Fill(4);
414   }
415
416   // loop over Cascades
417   for (Int_t iCascade = 0; iCascade < esd->GetNumberOfCascades(); iCascade++) {
418     AliESDcascade* cascade = esd->GetCascade(iCascade);
419     Double_t v0q;
420     cascade->ChangeMassHypothesis(v0q,kXiMinus);
421     hMassXi->Fill(cascade->GetEffMass());
422     cascade->ChangeMassHypothesis(v0q,kOmegaMinus);
423     hMassOmega->Fill(cascade->GetEffMass());
424     
425     Int_t negLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex())->GetLabel());
426     if (negLabel > mcEvent->GetNumberOfTracks()) continue;     // background
427     AliMCParticle* track = mcEvent->GetTrack(negLabel);
428     TParticle* particle = track->Particle();
429     Int_t negMother = particle->GetMother(0);
430     if (negMother < 0) continue;
431     Int_t posLabel = TMath::Abs(esd->GetTrack(cascade->GetPindex())->GetLabel());
432     if (posLabel > mcEvent->GetNumberOfTracks()) continue;     // background
433     track = mcEvent->GetTrack(posLabel);
434     particle = track->Particle();
435     Int_t posMother = particle->GetMother(0);
436     if (negMother != posMother) continue;
437     track = mcEvent->GetTrack(negMother);
438     particle = track->Particle();
439     Int_t v0Mother = particle->GetMother(0);
440     if (v0Mother < 0) continue;
441     Int_t bacLabel = TMath::Abs(esd->GetTrack(cascade->GetBindex())->GetLabel());
442     if (bacLabel > mcEvent->GetNumberOfTracks()) continue;     // background
443     track = mcEvent->GetTrack(bacLabel);
444     particle = track->Particle();
445     Int_t bacMother = particle->GetMother(0);
446     if (v0Mother != bacMother) continue;
447     track = mcEvent->GetTrack(v0Mother);
448     particle = track->Particle();
449     if (!selCascades.Contains(particle)) continue;
450     selCascades.Remove(particle);
451     hScalars->Fill(7);
452   }
453   
454   // loop over the clusters
455   for (Int_t iCluster=0; iCluster<esd->GetNumberOfCaloClusters(); iCluster++) {
456     AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster);
457     if (clust->IsPHOS()) hEPHOS->Fill(clust->E());
458     if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E());
459   }
460         
461   // Post output data.
462   PostData(1, fListOfHistos);
463 }
464
465 void AliAnalysisTaskCheckESD::Terminate(Option_t *)
466 {
467   // check values
468   Int_t    checkNGenLow = 1;
469   Double_t checkEffLow = 0.5;
470   Double_t checkEffSigma = 3;
471   Double_t checkFakeHigh = 0.5;
472   Double_t checkFakeSigma = 3;
473   
474   Double_t checkResPtInvHigh = 5;
475   Double_t checkResPtInvSigma = 3;
476   Double_t checkResPhiHigh = 10;
477   Double_t checkResPhiSigma = 3;
478   Double_t checkResThetaHigh = 10;
479   Double_t checkResThetaSigma = 3;
480   
481   Double_t checkPIDEffLow = 0.5;
482   Double_t checkPIDEffSigma = 3;
483   Double_t checkResTOFHigh = 500;
484   Double_t checkResTOFSigma = 3;
485   
486   Double_t checkPHOSNLow = 5;
487   Double_t checkPHOSEnergyLow = 0.3;
488   Double_t checkPHOSEnergyHigh = 1.0;
489   Double_t checkEMCALNLow = 50;
490   Double_t checkEMCALEnergyLow = 0.05;
491   Double_t checkEMCALEnergyHigh = 1.0;
492   
493   Double_t checkMUONNLow = 1;
494   Double_t checkMUONPtLow = 0.5;
495   Double_t checkMUONPtHigh = 10.;
496   
497   Double_t checkV0EffLow = 0.02;
498   Double_t checkV0EffSigma = 3;
499   Double_t checkCascadeEffLow = 0.01;
500   Double_t checkCascadeEffSigma = 3;
501   
502   const char* partName[AliPID::kSPECIES+1] = {"electron", "muon", "pion", "kaon", "proton", "other"};
503   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
504   if (!fListOfHistos) {
505     Printf("ERROR: fListOfHistos not available");
506     return;
507   }
508         
509   hGen = dynamic_cast<TH1F*>(fListOfHistos->At(0));
510   hRec = dynamic_cast<TH1F*>(fListOfHistos->At(1));
511   hResPtInv = dynamic_cast<TH1F*>(fListOfHistos->At(2));
512   hResPhi = dynamic_cast<TH1F*>(fListOfHistos->At(3));
513   hResTheta = dynamic_cast<TH1F*>(fListOfHistos->At(4));
514   hDEdxRight = dynamic_cast<TH2F*>(fListOfHistos->At(5));
515   hDEdxWrong = dynamic_cast<TH2F*>(fListOfHistos->At(6));
516   hResTOFRight = dynamic_cast<TH1F*>(fListOfHistos->At(7));
517   hResTOFWrong = dynamic_cast<TH1F*>(fListOfHistos->At(8));
518   hEPHOS = dynamic_cast<TH1F*>(fListOfHistos->At(9));
519   hEEMCAL = dynamic_cast<TH1F*>(fListOfHistos->At(10));
520   hPtMUON = dynamic_cast<TH1F*>(fListOfHistos->At(11));
521   hMassK0 = dynamic_cast<TH1F*>(fListOfHistos->At(12));
522   hMassLambda = dynamic_cast<TH1F*>(fListOfHistos->At(13));
523   hMassLambdaBar = dynamic_cast<TH1F*>(fListOfHistos->At(14));
524   hMassXi = dynamic_cast<TH1F*>(fListOfHistos->At(15));
525   hMassOmega = dynamic_cast<TH1F*>(fListOfHistos->At(16));
526   hScalars = dynamic_cast<TH1F*>(fListOfHistos->At(17));
527   hArrayHist = dynamic_cast<TH1F*>(fListOfHistos->At(18));
528   
529   Int_t nGen = Int_t(hScalars->GetBinContent(1));
530   Int_t nRec = Int_t(hScalars->GetBinContent(2));
531   Int_t nFake = Int_t(hScalars->GetBinContent(3));
532   Int_t nGenV0s = Int_t(hScalars->GetBinContent(4));
533   Int_t nRecV0s = Int_t(hScalars->GetBinContent(5));
534   Int_t nIdentified = Int_t(hScalars->GetBinContent(6));
535   Int_t nGenCascades = Int_t(hScalars->GetBinContent(7));
536   Int_t nRecCascades = Int_t(hScalars->GetBinContent(8));
537   
538   Int_t k = 1;
539   
540   Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES];
541   for(Int_t i = 0; i < (AliPID::kSPECIES+1); i++)
542     for(Int_t j = 0; j < AliPID::kSPECIES; j++) {
543       identified[i][j] = Int_t(hArrayHist->GetBinContent(k));
544       k++;
545     }
546   
547
548   // perform checks
549   if (nGen < checkNGenLow) {
550     Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen));
551   }
552         
553   TH1F* hEff = CreateEffHisto(hGen, hRec);
554   
555   Info("CheckESD", "%d out of %d tracks reconstructed including %d "
556        "fake tracks", nRec, nGen, nFake);
557   if (nGen > 0) {
558     // efficiency
559     Double_t eff = nRec*1./nGen;
560     Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGen);
561     Double_t fake = nFake*1./nGen;
562     Double_t fakeError = TMath::Sqrt(fake*(1.-fake) / nGen);
563     Info("CheckESD", "eff = (%.1f +- %.1f) %%  fake = (%.1f +- %.1f) %%",
564          100.*eff, 100.*effError, 100.*fake, 100.*fakeError);
565     if (eff < checkEffLow - checkEffSigma*effError) {
566       Warning("CheckESD", "low efficiency: (%.1f +- %.1f) %%",100.*eff, 100.*effError);
567     }
568     if (fake > checkFakeHigh + checkFakeSigma*fakeError) {
569       Warning("CheckESD", "high fake: (%.1f +- %.1f) %%",100.*fake, 100.*fakeError);
570     }
571     // resolutions
572     Double_t res, resError;
573     if (FitHisto(hResPtInv, res, resError)) {
574       Info("CheckESD", "relative inverse pt resolution = (%.1f +- %.1f) %%",res, resError);
575       if (res > checkResPtInvHigh + checkResPtInvSigma*resError) {
576         Warning("CheckESD", "bad pt resolution: (%.1f +- %.1f) %%",res, resError);
577       }
578     }
579     if (FitHisto(hResPhi, res, resError)) {
580       Info("CheckESD", "phi resolution = (%.1f +- %.1f) mrad", res, resError);
581       if (res > checkResPhiHigh + checkResPhiSigma*resError) {
582         Warning("CheckESD", "bad phi resolution: (%.1f +- %.1f) mrad", 
583                 res, resError);
584       }
585     }
586
587     if (FitHisto(hResTheta, res, resError)) {
588       Info("CheckESD", "theta resolution = (%.1f +- %.1f) mrad", 
589            res, resError);
590       if (res > checkResThetaHigh + checkResThetaSigma*resError) {
591         Warning("CheckESD", "bad theta resolution: (%.1f +- %.1f) mrad", 
592                 res, resError);
593       }
594     }
595     // PID
596     if (nRec > 0) {
597       Double_t eff = nIdentified*1./nRec;
598       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nRec);
599       Info("CheckESD", "PID eff = (%.1f +- %.1f) %%", 
600            100.*eff, 100.*effError);
601       if (eff < checkPIDEffLow - checkPIDEffSigma*effError) {
602         Warning("CheckESD", "low PID efficiency: (%.1f +- %.1f) %%", 
603                 100.*eff, 100.*effError);
604       }
605     }
606
607     printf("%9s:", "gen\\rec");
608     for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {
609       printf("%9s", partName[iRec]);
610     }
611     printf("\n");
612     for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) {
613       printf("%9s:", partName[iGen]);
614       for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {
615         printf("%9d", identified[iGen][iRec]);
616       }
617       printf("\n");
618     }
619     
620     if (FitHisto(hResTOFRight, res, resError)) {
621       Info("CheckESD", "TOF resolution = (%.1f +- %.1f) ps", res, resError);
622       if (res > checkResTOFHigh + checkResTOFSigma*resError) {
623         Warning("CheckESD", "bad TOF resolution: (%.1f +- %.1f) ps", 
624                 res, resError);
625       }
626     }
627     
628     // calorimeters
629     if (hEPHOS->Integral() < checkPHOSNLow) {
630       Warning("CheckESD", "low number of PHOS particles: %d", 
631               Int_t(hEPHOS->Integral()));
632     }   
633     else {
634       Double_t mean = hEPHOS->GetMean();
635       if (mean < checkPHOSEnergyLow) {
636         Warning("CheckESD", "low mean PHOS energy: %.1f GeV", mean);
637       } else if (mean > checkPHOSEnergyHigh) {
638         Warning("CheckESD", "high mean PHOS energy: %.1f GeV", mean);
639       }
640     }
641     
642     if (hEEMCAL->Integral() < checkEMCALNLow) {
643       Warning("CheckESD", "low number of EMCAL particles: %d", 
644               Int_t(hEEMCAL->Integral()));
645     } 
646     else {
647       Double_t mean = hEEMCAL->GetMean();
648       if (mean < checkEMCALEnergyLow) {
649         Warning("CheckESD", "low mean EMCAL energy: %.1f GeV", mean);
650       }
651       else if (mean > checkEMCALEnergyHigh) {
652         Warning("CheckESD", "high mean EMCAL energy: %.1f GeV", mean);
653       }
654     }
655
656     // muons
657     if (hPtMUON->Integral() < checkMUONNLow) {
658       Warning("CheckESD", "low number of MUON particles: %d", 
659               Int_t(hPtMUON->Integral()));
660     }
661     else {
662       Double_t mean = hPtMUON->GetMean();
663       if (mean < checkMUONPtLow) {
664         Warning("CheckESD", "low mean MUON pt: %.1f GeV/c", mean);
665       }
666       else if (mean > checkMUONPtHigh) {
667         Warning("CheckESD", "high mean MUON pt: %.1f GeV/c", mean);
668       }
669     }
670
671     // V0s
672     if (nGenV0s > 0) {
673       Double_t eff = nRecV0s*1./nGenV0s;
674       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenV0s);
675       if (effError == 0) effError = checkV0EffLow / TMath::Sqrt(1.*nGenV0s);
676       Info("CheckESD", "V0 eff = (%.1f +- %.1f) %%", 
677            100.*eff, 100.*effError);
678       if (eff < checkV0EffLow - checkV0EffSigma*effError) {
679         Warning("CheckESD", "low V0 efficiency: (%.1f +- %.1f) %%", 
680                 100.*eff, 100.*effError);
681       }
682     }
683
684     // Cascades
685     if (nGenCascades > 0) {
686       Double_t eff = nRecCascades*1./nGenCascades;
687       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenCascades);
688       if (effError == 0) effError = checkV0EffLow / 
689                            TMath::Sqrt(1.*nGenCascades);
690       Info("CheckESD", "Cascade eff = (%.1f +- %.1f) %%", 
691            100.*eff, 100.*effError);
692       if (eff < checkCascadeEffLow - checkCascadeEffSigma*effError) {
693         Warning("CheckESD", "low Cascade efficiency: (%.1f +- %.1f) %%", 
694                 100.*eff, 100.*effError);
695       }
696     }
697     
698   }
699
700   // draw the histograms if not in batch mode
701   if (!gROOT->IsBatch()) {
702     new TCanvas;
703     hEff->DrawCopy();
704     new TCanvas;
705     hResPtInv->DrawCopy("E");
706     new TCanvas;
707     hResPhi->DrawCopy("E");
708     new TCanvas;
709     hResTheta->DrawCopy("E");
710     new TCanvas;
711     hDEdxRight->DrawCopy();
712     hDEdxWrong->DrawCopy("SAME");
713     new TCanvas;
714     hResTOFRight->DrawCopy("E");
715     hResTOFWrong->DrawCopy("SAME");
716     new TCanvas;
717     hEPHOS->DrawCopy("E");
718     new TCanvas;
719     hEEMCAL->DrawCopy("E");
720     new TCanvas;
721     hPtMUON->DrawCopy("E");
722     new TCanvas;
723     hMassK0->DrawCopy("E");
724     new TCanvas;
725     hMassLambda->DrawCopy("E");
726     new TCanvas;
727     hMassLambdaBar->DrawCopy("E");
728     new TCanvas;
729     hMassXi->DrawCopy("E");
730     new TCanvas;
731     hMassOmega->DrawCopy("E");
732   }
733
734   // write the output histograms to a file
735   TFile* outputFile = TFile::Open("check.root", "recreate");
736   if (!outputFile || !outputFile->IsOpen())
737     {
738       Error("CheckESD", "opening output file check.root failed");
739       return;
740     }
741   hEff->Write();
742   hResPtInv->Write();
743   hResPhi->Write();
744   hResTheta->Write();
745   hDEdxRight->Write();
746   hDEdxWrong->Write();
747   hResTOFRight->Write();
748   hResTOFWrong->Write();
749   hEPHOS->Write();
750   hEEMCAL->Write();
751   hPtMUON->Write();
752   hMassK0->Write();
753   hMassLambda->Write();
754   hMassLambdaBar->Write();
755   hMassXi->Write();
756   hMassOmega->Write();
757   outputFile->Close();
758   delete outputFile;
759
760   // clean up
761   /*    delete hGen;
762         delete hRec;
763         delete hEff;
764         delete hResPtInv;
765         delete hResPhi;
766         delete hResTheta;
767         delete hDEdxRight;
768         delete hDEdxWrong;
769         delete hResTOFRight;
770         delete hResTOFWrong;
771         delete hEPHOS;
772         delete hEEMCAL;
773         delete hPtMUON;
774         delete hMassK0;
775         delete hMassLambda;
776         delete hMassLambdaBar;
777         delete hMassXi;
778         delete hMassOmega;
779         delete hScalars;
780         delete hArrayHist; */
781
782   //delete esd;
783   // result of check
784   Info("CheckESD", "check of ESD was successfull");
785 }
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847