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