]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/sim/Check.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Various small fixes:
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / sim / Check.C
CommitLineData
421f877c 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 */
f8b7a926 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
ee275c29 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 */
50TH1*
51CreateHisto(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)
f8b7a926 58{
ee275c29 59 // create a histogram
f8b7a926 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
ee275c29 68//____________________________________________________________________
69/**
70 * Create efficiency histogram
71 *
72 * @param hGen Generated
73 * @param hRec Reconstructed
74 *
75 * @return Newly allocated histogram
76 */
77TH1* CreateEffHisto(TH1* hGen, TH1* hRec)
f8b7a926 78{
ee275c29 79 // create an efficiency histogram
80
f8b7a926 81 Int_t nBins = hGen->GetNbinsX();
ee275c29 82 TH1* hEff = static_cast<TH1*>(hGen->Clone("hEff"));
f8b7a926 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
ee275c29 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 */
117Bool_t
118FitHisto(TH1* histo, Double_t& res, Double_t& resError)
f8b7a926 119{
f8b7a926 120 static TF1* fitFunc = new TF1("fitFunc", "gaus");
121 fitFunc->SetLineWidth(2);
122 fitFunc->SetFillStyle(0);
123 Double_t maxFitRange = 2;
124
ee275c29 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;
f8b7a926 136}
137
ee275c29 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 */
f8b7a926 147Bool_t Check(const char* gAliceFileName = "galice.root",
148 const char* esdFileName = "AliESDs.root")
149{
ee275c29 150 // check the content of the ESD
f8b7a926 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
15cb4a2b 379 hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) -
380 1./particle->Pt()) *
f8b7a926 381 particle->Pt());
382 hResPhi->Fill(1000. * (track->Phi() - particle->Phi()));
383 hResTheta->Fill(1000. * (track->Theta() - particle->Theta()));
384
385 // PID
386 if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
387 Int_t iGen = 5;
388 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
389 if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i;
390 }
391 Double_t probability[AliPID::kSPECIES];
392 track->GetESDpid(probability);
393 Double_t pMax = 0;
394 Int_t iRec = 0;
395 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
396 probability[i] *= partFrac[i];
397 if (probability[i] > pMax) {
398 pMax = probability[i];
399 iRec = i;
400 }
401 }
402 identified[iGen][iRec]++;
403 if (iGen == iRec) nIdentified++;
404
405 // dE/dx and TOF
406 Double_t time[AliPID::kSPECIES];
407 track->GetIntegratedTimes(time);
408 if (iGen == iRec) {
409 hDEdxRight->Fill(particle->P(), track->GetTPCsignal());
410 if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {
411 hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]);
412 }
413 } else {
414 hDEdxWrong->Fill(particle->P(), track->GetTPCsignal());
415 if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {
416 hResTOFWrong->Fill(track->GetTOFsignal() - time[iRec]);
417 }
418 }
419 }
420
421 // loop over muon tracks
422 {
423 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
424 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
425 Double_t ptInv = TMath::Abs(muonTrack->GetInverseBendingMomentum());
426 if (ptInv > 0.001) {
427 hPtMUON->Fill(1./ptInv);
428 }
429 }
430 }
431
432 // loop over V0s
433 for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
434 AliESDv0* v0 = esd->GetV0(iV0);
435 if (v0->GetOnFlyStatus()) continue;
436 v0->ChangeMassHypothesis(kK0Short);
437 hMassK0->Fill(v0->GetEffMass());
438 v0->ChangeMassHypothesis(kLambda0);
439 hMassLambda->Fill(v0->GetEffMass());
440 v0->ChangeMassHypothesis(kLambda0Bar);
441 hMassLambdaBar->Fill(v0->GetEffMass());
442
443 Int_t negLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel());
444 if (negLabel > stack->GetNtrack()) continue; // background
445 Int_t negMother = stack->Particle(negLabel)->GetMother(0);
446 if (negMother < 0) continue;
447 Int_t posLabel = TMath::Abs(esd->GetTrack(v0->GetPindex())->GetLabel());
448 if (posLabel > stack->GetNtrack()) continue; // background
449 Int_t posMother = stack->Particle(posLabel)->GetMother(0);
450 if (negMother != posMother) continue;
451 TParticle* particle = stack->Particle(negMother);
452 if (!selV0s.Contains(particle)) continue;
453 selV0s.Remove(particle);
454 nRecV0s++;
455 }
456
457 // loop over Cascades
458 for (Int_t iCascade = 0; iCascade < esd->GetNumberOfCascades();
459 iCascade++) {
460 AliESDcascade* cascade = esd->GetCascade(iCascade);
461 Double_t v0q;
462 cascade->ChangeMassHypothesis(v0q,kXiMinus);
463 hMassXi->Fill(cascade->GetEffMassXi());
464 cascade->ChangeMassHypothesis(v0q,kOmegaMinus);
465 hMassOmega->Fill(cascade->GetEffMassXi());
466
467 Int_t negLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex())
468 ->GetLabel());
469 if (negLabel > stack->GetNtrack()) continue; // background
470 Int_t negMother = stack->Particle(negLabel)->GetMother(0);
471 if (negMother < 0) continue;
472 Int_t posLabel = TMath::Abs(esd->GetTrack(cascade->GetPindex())
473 ->GetLabel());
474 if (posLabel > stack->GetNtrack()) continue; // background
475 Int_t posMother = stack->Particle(posLabel)->GetMother(0);
476 if (negMother != posMother) continue;
477 Int_t v0Mother = stack->Particle(negMother)->GetMother(0);
478 if (v0Mother < 0) continue;
479 Int_t bacLabel = TMath::Abs(esd->GetTrack(cascade->GetBindex())
480 ->GetLabel());
481 if (bacLabel > stack->GetNtrack()) continue; // background
482 Int_t bacMother = stack->Particle(bacLabel)->GetMother(0);
483 if (v0Mother != bacMother) continue;
484 TParticle* particle = stack->Particle(v0Mother);
485 if (!selCascades.Contains(particle)) continue;
486 selCascades.Remove(particle);
487 nRecCascades++;
488 }
489
490 // loop over the clusters
491 {
15cb4a2b 492 for (Int_t iCluster=0;iCluster<esd->GetNumberOfCaloClusters();iCluster++){
f8b7a926 493 AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster);
494 if (clust->IsPHOS()) hEPHOS->Fill(clust->E());
495 if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E());
496 }
497 }
498
499 }
500
501 // perform checks
502 if (nGen < checkNGenLow) {
503 Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen));
504 }
505
506 TH1F* hEff = CreateEffHisto(hGen, hRec);
507
508 Info("CheckESD", "%d out of %d tracks reconstructed including %d "
509 "fake tracks", nRec, nGen, nFake);
510 if (nGen > 0) {
511 // efficiency
512 Double_t eff = nRec*1./nGen;
513 Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGen);
514 Double_t fake = nFake*1./nGen;
515 Double_t fakeError = TMath::Sqrt(fake*(1.-fake) / nGen);
516 Info("CheckESD", "eff = (%.1f +- %.1f) %% fake = (%.1f +- %.1f) %%",
517 100.*eff, 100.*effError, 100.*fake, 100.*fakeError);
518
519 if (eff < checkEffLow - checkEffSigma*effError) {
520 Warning("CheckESD", "low efficiency: (%.1f +- %.1f) %%",
521 100.*eff, 100.*effError);
522 }
523 if (fake > checkFakeHigh + checkFakeSigma*fakeError) {
524 Warning("CheckESD", "high fake: (%.1f +- %.1f) %%",
525 100.*fake, 100.*fakeError);
526 }
527
528 // resolutions
529 Double_t res, resError;
530 if (FitHisto(hResPtInv, res, resError)) {
531 Info("CheckESD", "relative inverse pt resolution = (%.1f +- %.1f) %%",
532 res, resError);
533 if (res > checkResPtInvHigh + checkResPtInvSigma*resError) {
534 Warning("CheckESD", "bad pt resolution: (%.1f +- %.1f) %%",
535 res, resError);
536 }
537 }
538
539 if (FitHisto(hResPhi, res, resError)) {
540 Info("CheckESD", "phi resolution = (%.1f +- %.1f) mrad", res, resError);
541 if (res > checkResPhiHigh + checkResPhiSigma*resError) {
542 Warning("CheckESD", "bad phi resolution: (%.1f +- %.1f) mrad",
543 res, resError);
544 }
545 }
546
547 if (FitHisto(hResTheta, res, resError)) {
548 Info("CheckESD", "theta resolution = (%.1f +- %.1f) mrad",
549 res, resError);
550 if (res > checkResThetaHigh + checkResThetaSigma*resError) {
551 Warning("CheckESD", "bad theta resolution: (%.1f +- %.1f) mrad",
552 res, resError);
553 }
554 }
555
556 // PID
557 if (nRec > 0) {
558 Double_t eff = nIdentified*1./nRec;
559 Double_t effError = TMath::Sqrt(eff*(1.-eff) / nRec);
560 Info("CheckESD", "PID eff = (%.1f +- %.1f) %%",
561 100.*eff, 100.*effError);
562 if (eff < checkPIDEffLow - checkPIDEffSigma*effError) {
563 Warning("CheckESD", "low PID efficiency: (%.1f +- %.1f) %%",
564 100.*eff, 100.*effError);
565 }
566 }
567
568 printf("%9s:", "gen\\rec");
569 for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {
570 printf("%9s", partName[iRec]);
571 }
572 printf("\n");
573 for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) {
574 printf("%9s:", partName[iGen]);
575 for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {
576 printf("%9d", identified[iGen][iRec]);
577 }
578 printf("\n");
579 }
580
581 if (FitHisto(hResTOFRight, res, resError)) {
582 Info("CheckESD", "TOF resolution = (%.1f +- %.1f) ps", res, resError);
583 if (res > checkResTOFHigh + checkResTOFSigma*resError) {
584 Warning("CheckESD", "bad TOF resolution: (%.1f +- %.1f) ps",
585 res, resError);
586 }
587 }
588
589 // calorimeters
590 if (hEPHOS->Integral() < checkPHOSNLow) {
591 Warning("CheckESD", "low number of PHOS particles: %d",
592 Int_t(hEPHOS->Integral()));
593 } else {
594 Double_t mean = hEPHOS->GetMean();
595 if (mean < checkPHOSEnergyLow) {
596 Warning("CheckESD", "low mean PHOS energy: %.1f GeV", mean);
597 } else if (mean > checkPHOSEnergyHigh) {
598 Warning("CheckESD", "high mean PHOS energy: %.1f GeV", mean);
599 }
600 }
601
602 if (hEEMCAL->Integral() < checkEMCALNLow) {
603 Warning("CheckESD", "low number of EMCAL particles: %d",
604 Int_t(hEEMCAL->Integral()));
605 } else {
606 Double_t mean = hEEMCAL->GetMean();
607 if (mean < checkEMCALEnergyLow) {
608 Warning("CheckESD", "low mean EMCAL energy: %.1f GeV", mean);
609 } else if (mean > checkEMCALEnergyHigh) {
610 Warning("CheckESD", "high mean EMCAL energy: %.1f GeV", mean);
611 }
612 }
613
614 // muons
615 if (hPtMUON->Integral() < checkMUONNLow) {
616 Warning("CheckESD", "low number of MUON particles: %d",
617 Int_t(hPtMUON->Integral()));
618 } else {
619 Double_t mean = hPtMUON->GetMean();
620 if (mean < checkMUONPtLow) {
621 Warning("CheckESD", "low mean MUON pt: %.1f GeV/c", mean);
622 } else if (mean > checkMUONPtHigh) {
623 Warning("CheckESD", "high mean MUON pt: %.1f GeV/c", mean);
624 }
625 }
626
627 // V0s
628 if (nGenV0s > 0) {
629 Double_t eff = nRecV0s*1./nGenV0s;
630 Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenV0s);
631 if (effError == 0) effError = checkV0EffLow / TMath::Sqrt(1.*nGenV0s);
632 Info("CheckESD", "V0 eff = (%.1f +- %.1f) %%",
633 100.*eff, 100.*effError);
634 if (eff < checkV0EffLow - checkV0EffSigma*effError) {
635 Warning("CheckESD", "low V0 efficiency: (%.1f +- %.1f) %%",
636 100.*eff, 100.*effError);
637 }
638 }
639
640 // Cascades
641 if (nGenCascades > 0) {
642 Double_t eff = nRecCascades*1./nGenCascades;
643 Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenCascades);
644 if (effError == 0) effError = checkV0EffLow /
645 TMath::Sqrt(1.*nGenCascades);
646 Info("CheckESD", "Cascade eff = (%.1f +- %.1f) %%",
647 100.*eff, 100.*effError);
648 if (eff < checkCascadeEffLow - checkCascadeEffSigma*effError) {
649 Warning("CheckESD", "low Cascade efficiency: (%.1f +- %.1f) %%",
650 100.*eff, 100.*effError);
651 }
652 }
653 }
654
655 // draw the histograms if not in batch mode
656 if (!gROOT->IsBatch()) {
657 new TCanvas;
658 hEff->DrawCopy();
659 new TCanvas;
660 hResPtInv->DrawCopy("E");
661 new TCanvas;
662 hResPhi->DrawCopy("E");
663 new TCanvas;
664 hResTheta->DrawCopy("E");
665 new TCanvas;
666 hDEdxRight->DrawCopy();
667 hDEdxWrong->DrawCopy("SAME");
668 new TCanvas;
669 hResTOFRight->DrawCopy("E");
670 hResTOFWrong->DrawCopy("SAME");
671 new TCanvas;
672 hEPHOS->DrawCopy("E");
673 new TCanvas;
674 hEEMCAL->DrawCopy("E");
675 new TCanvas;
676 hPtMUON->DrawCopy("E");
677 new TCanvas;
678 hMassK0->DrawCopy("E");
679 new TCanvas;
680 hMassLambda->DrawCopy("E");
681 new TCanvas;
682 hMassLambdaBar->DrawCopy("E");
683 new TCanvas;
684 hMassXi->DrawCopy("E");
685 new TCanvas;
686 hMassOmega->DrawCopy("E");
687 }
688
689 // write the output histograms to a file
690 TFile* outputFile = TFile::Open("check.root", "recreate");
691 if (!outputFile || !outputFile->IsOpen()) {
692 Error("CheckESD", "opening output file check.root failed");
693 return kFALSE;
694 }
695 hEff->Write();
696 hResPtInv->Write();
697 hResPhi->Write();
698 hResTheta->Write();
699 hDEdxRight->Write();
700 hDEdxWrong->Write();
701 hResTOFRight->Write();
702 hResTOFWrong->Write();
703 hEPHOS->Write();
704 hEEMCAL->Write();
705 hPtMUON->Write();
706 hMassK0->Write();
707 hMassLambda->Write();
708 hMassLambdaBar->Write();
709 hMassXi->Write();
710 hMassOmega->Write();
711 outputFile->Close();
712 delete outputFile;
713
714 // clean up
715 delete hGen;
716 delete hRec;
717 delete hEff;
718 delete hResPtInv;
719 delete hResPhi;
720 delete hResTheta;
721 delete hDEdxRight;
722 delete hDEdxWrong;
723 delete hResTOFRight;
724 delete hResTOFWrong;
725 delete hEPHOS;
726 delete hEEMCAL;
727 delete hPtMUON;
728 delete hMassK0;
729 delete hMassLambda;
730 delete hMassLambdaBar;
731 delete hMassXi;
732 delete hMassOmega;
733
734 delete esd;
735 esdFile->Close();
736 delete esdFile;
737
738 runLoader->UnloadHeader();
739 runLoader->UnloadKinematics();
740 delete runLoader;
741
742 // result of check
743 Info("CheckESD", "check of ESD was successfull");
744 return kTRUE;
745}
ee275c29 746//
747// EOF
748//