prod/LHC09a14/CheckESD.C

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