Use PDG code incude
[u/mrichter/AliRoot.git] / PWG1 / comparison / AliAnalysisTaskCheckESD.cxx
CommitLineData
7d245e03 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
8858ee12 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
43ClassImp(AliAnalysisTaskCheckESD)
44
45AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD():
46AliAnalysisTaskSE("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
76AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD(const char* name):
77AliAnalysisTaskSE(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
108TH1F * 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
121TH1F * 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
151Bool_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
173void 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
8d071f77 211 hEPHOS = CreateHisto("hEPHOS", "PHOS", 100, 0, 50, "E [GeV]", "N");
8858ee12 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
249void AliAnalysisTaskCheckESD::UserExec(Option_t *option)
250{
8858ee12 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
480void 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