test/QA/CheckESD.C

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