AliPhysics  86c65ee (86c65ee)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Check.C
Go to the documentation of this file.
1 
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 //____________________________________________________________________
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 //____________________________________________________________________
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 //____________________________________________________________________
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 //====================================================================
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()) -
380  1./particle->Pt()) *
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  {
492  for (Int_t iCluster=0;iCluster<esd->GetNumberOfCaloClusters();iCluster++){
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 }
746 //
747 // EOF
748 //
Check
Bool_t Check(const char *gAliceFileName="galice.root", const char *esdFileName="AliESDs.root")
Definition: Check.C:147
Double_t
double Double_t
Definition: External.C:58
TH2F
Definition: External.C:236
title
const char * title
Definition: MakeQAPdf.C:26
Int_t
int Int_t
Definition: External.C:63
Float_t
float Float_t
Definition: External.C:68
CreateEffHisto
TH1 * CreateEffHisto(TH1 *hGen, TH1 *hRec)
Definition: Check.C:77
kProton
Definition: FastCentEstimators.C:23
CreateHisto
TH1 * CreateHisto(const char *name, const char *title, Int_t nBins, Double_t xMin, Double_t xMax, const char *xLabel=NULL, const char *yLabel=NULL)
Definition: Check.C:51
FitHisto
Bool_t FitHisto(TH1 *histo, Double_t &res, Double_t &resError)
Definition: Check.C:118
Bool_t
bool Bool_t
Definition: External.C:53
TH1
Definition: External.C:196