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