PWGPP/comparison/AliAnalysisTaskCheckESD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //------------------------------
  19 //
  20 // Proof-enabled
  21 // version
  22 // of CheckESD.C
  23 //
  24 //------------------------------
  25
  26 #include "TChain.h"
  27 #include "TROOT.h"
  28 #include "TFile.h"
  29 #include "TH1F.h"
  30 #include "TH2F.h"
  31 #include "TF1.h"
  32 #include "TCanvas.h"
  33 #include "TVector3.h"
  34 #include "TParticle.h"
  35 #include "AliESDEvent.h"
  36 #include "AliESDv0.h"
  37 #include "AliESDcascade.h"
  38 #include "AliESDMuonTrack.h"
  39 #include "AliESDCaloCluster.h"
  40 #include "AliRun.h"
  41 #include "AliMCEvent.h"
  42 #include "AliGenEventHeader.h"
  43 #include "AliPID.h"
  44
  45 #include "AliAnalysisTaskCheckESD.h"
  46
  47 ClassImp(AliAnalysisTaskCheckESD)
  48
  49 AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD():
  50 AliAnalysisTaskSE("AliAnalysisTaskCheckESD"),
  51   fListOfHistos(0),
  52   fGen(0),
  53   fRec(0),
  54   fResPtInv(0),
  55   fResPhi(0),
  56   fResTheta(0),
  57   fDEdxRight(0),
  58   fDEdxWrong(0),
  59   fResTOFRight(0),
  60   fResTOFWrong(0),
  61   fEPHOS(0),
  62   fEEMCAL(0),
  63   fPtMUON(0),
  64   fMassK0(0),
  65   fMassLambda(0),
  66   fMassLambdaBar(0),
  67   fMassXi(0),
  68   fMassOmega(0),
  69   fScalars(0),
  70   fArrayHist(0)
  71 {
  72   // Default constructor
  73   // Define input and output slots here
  74   // Input slot #0 works with a TChain
  75   DefineInput(0, TChain::Class());
  76   // Output slot #1 TList
  77   DefineOutput(1, TList::Class());
  78 }
  79
  80 AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD(const char* name):
  81 AliAnalysisTaskSE(name),
  82   fListOfHistos(0),
  83   fGen(0),
  84   fRec(0),
  85   fResPtInv(0),
  86   fResPhi(0),
  87   fResTheta(0),
  88   fDEdxRight(0),
  89   fDEdxWrong(0),
  90   fResTOFRight(0),
  91   fResTOFWrong(0),
  92   fEPHOS(0),
  93   fEEMCAL(0),
  94   fPtMUON(0),
  95   fMassK0(0),
  96   fMassLambda(0),
  97   fMassLambdaBar(0),
  98   fMassXi(0),
  99   fMassOmega(0),
 100   fScalars(0),
 101   fArrayHist(0)
 102 {
 103   // Constructor
 104   AliInfo("Constructor AliAnalysisTaskCheckESD");
 105   // Define input and output slots here
 106   // Input slot #0 works with a TChain
 107   DefineInput(0, TChain::Class());
 108   // Output slot #1 TList
 109   DefineOutput(1, TList::Class());
 110 }
 111
 112 TH1F * AliAnalysisTaskCheckESD::CreateHisto(const char* name, const char* title,Int_t nBins,
 113                                             Double_t xMin, Double_t xMax,
 114                                             const char* xLabel, const char* yLabel)
 115 {
 116   // create a histogram
 117   TH1F* result = new TH1F(name, title, nBins, xMin, xMax);
 118   result->SetOption("E");
 119   if (xLabel) result->GetXaxis()->SetTitle(xLabel);
 120   if (yLabel) result->GetYaxis()->SetTitle(yLabel);
 121   result->SetMarkerStyle(kFullCircle);
 122   return result;
 123 }
 124
 125 TH1F *AliAnalysisTaskCheckESD::CreateEffHisto(const TH1F* hGen, const TH1F* hRec)
 126 {
 127   // create an efficiency histogram
 128   Int_t nBins = hGen->GetNbinsX();
 129   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
 130   hEff->SetTitle("");
 131   hEff->SetStats(kFALSE);
 132   hEff->SetMinimum(0.);
 133   hEff->SetMaximum(110.);
 134   hEff->GetYaxis()->SetTitle("#epsilon [%]");
 135
 136   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
 137     Double_t nGenEff = hGen->GetBinContent(iBin);
 138     Double_t nRecEff = hRec->GetBinContent(iBin);
 139     if (nGenEff > 0) {
 140       Double_t eff = nRecEff/nGenEff;
 141       hEff->SetBinContent(iBin, 100. * eff);
 142       Double_t error = sqrt(eff*(1.-eff) / nGenEff);
 143       if (TMath::Abs(error) < 1.e-33) error = 0.0001;
 144       hEff->SetBinError(iBin, 100. * error);
 145     }
 146     else {
 147       hEff->SetBinContent(iBin, -100.);
 148       hEff->SetBinError(iBin, 0);
 149     }
 150   }
 151   return hEff;
 152 }
 153
 154
 155 Bool_t AliAnalysisTaskCheckESD::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
 156 {
 157   // fit a gaussian to a histogram
 158   static TF1* fitFunc = new TF1("fitFunc", "gaus");
 159   fitFunc->SetLineWidth(2);
 160   fitFunc->SetFillStyle(0);
 161   Double_t maxFitRange = 2;
 162
 163   if (histo->Integral() > 50) {
 164     Float_t mean = histo->GetMean();
 165     Float_t rms = histo->GetRMS();
 166     fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms);
 167     fitFunc->SetParameters(mean, rms);
 168     histo->Fit(fitFunc, "QRI0");
 169     histo->GetFunction("fitFunc")->ResetBit(1<<9);
 170     res = TMath::Abs(fitFunc->GetParameter(2));
 171     resError = TMath::Abs(fitFunc->GetParError(2));
 172     return kTRUE;
 173   }
 174   return kFALSE;
 175 }
 176
 177 void AliAnalysisTaskCheckESD::UserCreateOutputObjects()
 178 {
 179   // Create histograms
 180   AliInfo("AliAnalysisTaskCheckESD::UserCreateOutputObjects");
 181   // Create output container
 182   fListOfHistos = new TList();
 183
 184   // efficiency and resolution histog
 185   Int_t nBinsPt = 15;
 186   Float_t minPt = 0.1;
 187   Float_t maxPt = 3.1;
 188
 189   fGen = CreateHisto("hGen", "generated tracks", nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N");
 190   fRec = CreateHisto("hRec", "reconstructed tracks", nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N");
 191   fResPtInv = CreateHisto("hResPtInv", "", 100, -10, 10, "(p_{t,rec}^{-1}-p_{t,sim}^{-1}) / p_{t,sim}^{-1} [%]", "N");
 192   fResPhi = CreateHisto("hResPhi", "", 100, -20, 20, "#phi_{rec}-#phi_{sim} [mrad]", "N");
 193   fResTheta = CreateHisto("hResTheta", "", 100, -20, 20, "#theta_{rec}-#theta_{sim} [mrad]", "N");
 194
 195   // dE/dx and TOF
 196   fDEdxRight = new TH2F("hDEdxRight", "", 300, 0, 3, 100, 0, 400);
 197   fDEdxRight->SetStats(kFALSE);
 198   fDEdxRight->GetXaxis()->SetTitle("p [GeV/c]");
 199   fDEdxRight->GetYaxis()->SetTitle("dE/dx_{TPC}");
 200   fDEdxRight->SetMarkerStyle(kFullCircle);
 201   fDEdxRight->SetMarkerSize(0.4);
 202   fDEdxWrong = new TH2F("hDEdxWrong", "", 300, 0, 3, 100, 0, 400);
 203   fDEdxWrong->SetStats(kFALSE);
 204   fDEdxWrong->GetXaxis()->SetTitle("p [GeV/c]");
 205   fDEdxWrong->GetYaxis()->SetTitle("dE/dx_{TPC}");
 206   fDEdxWrong->SetMarkerStyle(kFullCircle);
 207   fDEdxWrong->SetMarkerSize(0.4);
 208   fDEdxWrong->SetMarkerColor(kRed);
 209
 210   fResTOFRight = CreateHisto("hResTOFRight", "", 100, -1000, 1000, "t_{TOF}-t_{track} [ps]", "N");
 211   fResTOFWrong = CreateHisto("hResTOFWrong", "", 100, -1000, 1000, "t_{TOF}-t_{track} [ps]", "N");
 212   fResTOFWrong->SetLineColor(kRed);
 213
 214   // calorimeters
 215   fEPHOS = CreateHisto("hEPHOS", "PHOS", 100, 0, 50, "E [GeV]", "N");
 216   fEEMCAL = CreateHisto("hEEMCAL", "EMCAL", 100, 0, 50, "E [GeV]", "N");
 217
 218   // muons
 219   fPtMUON = CreateHisto("hPtMUON", "MUON", 100, 0, 20, "p_{t} [GeV/c]", "N");
 220
 221   // V0s and cascades
 222   fMassK0 = CreateHisto("hMassK0", "K^{0}", 100, 0.4, 0.6, "M(#pi^{+}#pi^{-}) [GeV/c^{2}]", "N");
 223   fMassLambda = CreateHisto("hMassLambda", "#Lambda", 100, 1.0, 1.2, "M(p#pi^{-}) [GeV/c^{2}]", "N");
 224
 225   fMassLambdaBar = CreateHisto("hMassLambdaBar", "#bar{#Lambda}", 100, 1.0, 1.2, "M(#bar{p}#pi^{+}) [GeV/c^{2}]", "N");
 226   fMassXi = CreateHisto("hMassXi", "#Xi", 100, 1.2, 1.5, "M(#Lambda#pi) [GeV/c^{2}]", "N");
 227   fMassOmega = CreateHisto("hMassOmega", "#Omega", 100, 1.5, 1.8, "M(#LambdaK) [GeV/c^{2}]", "N");
 228   fScalars = new TH1F("hScalars","Container of scalars",8,0,8);
 229   fArrayHist = new TH1F("hArrayHist","Container for Array",
 230                         (AliPID::kSPECIES+1)*AliPID::kSPECIES,0,(AliPID::kSPECIES+1)*AliPID::kSPECIES);
 231
 232   fListOfHistos->Add(fGen);
 233   fListOfHistos->Add(fRec);
 234   fListOfHistos->Add(fResPtInv);
 235   fListOfHistos->Add(fResPhi);
 236   fListOfHistos->Add(fResTheta);
 237   fListOfHistos->Add(fDEdxRight);
 238   fListOfHistos->Add(fDEdxWrong);
 239   fListOfHistos->Add(fResTOFRight);
 240   fListOfHistos->Add(fResTOFWrong);
 241   fListOfHistos->Add(fEPHOS);
 242   fListOfHistos->Add(fEEMCAL);
 243   fListOfHistos->Add(fPtMUON);
 244   fListOfHistos->Add(fMassK0);
 245   fListOfHistos->Add(fMassLambda);
 246   fListOfHistos->Add(fMassLambdaBar);
 247   fListOfHistos->Add(fMassXi);
 248   fListOfHistos->Add(fMassOmega);
 249   fListOfHistos->Add(fScalars);
 250   fListOfHistos->Add(fArrayHist);
 251 }
 252
 253 void AliAnalysisTaskCheckESD::UserExec(Option_t */*option*/)
 254 {
 255   // check the content of the ESD
 256   Double_t cutPtV0 = 0.3;
 257   Double_t cutPtCascade = 0.5;
 258   Float_t minPt = 0.1;
 259
 260   // PID
 261   Int_t partCode[AliPID::kSPECIES] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
 262   Double_t partFrac[AliPID::kSPECIES] = {0.01, 0.01, 0.85, 0.10, 0.05};
 263
 264   AliMCEvent* mcEvent = MCEvent();
 265   if (!mcEvent) {
 266     Printf("ERROR: Could not retrieve MC event");
 267     return;
 268   }
 269
 270   //Primary vertex needed
 271   TArrayF vertex(3);
 272   mcEvent->GenEventHeader()->PrimaryVertex(vertex);
 273
 274   TObjArray selParticles; //Use TClonesArray?
 275   TObjArray selV0s;
 276   TObjArray selCascades;
 277   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again
 278     AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
 279     TParticle* particle = track->Particle();
 280     if (!particle) continue;
 281     if (particle->Pt() < 0.001) continue;
 282     if (TMath::Abs(particle->Eta()) > 0.9) continue;
 283     TVector3 dVertex(particle->Vx() - vertex[0], particle->Vy() - vertex[1], particle->Vz() - vertex[2]);
 284     if (dVertex.Mag() > 0.0001) continue;
 285
 286     switch (TMath::Abs(particle->GetPdgCode())) {
 287     case kElectron:
 288     case kMuonMinus:
 289     case kPiPlus:
 290     case kKPlus:
 291     case kProton:
 292       {
 293         if (particle->Pt() > minPt) {
 294           selParticles.Add(particle);
 295           fScalars->Fill(0);
 296           fGen->Fill(particle->Pt());
 297         }
 298         break;
 299       }
 300     case kK0Short:
 301     case kLambda0:
 302       {
 303         if (particle->Pt() > cutPtV0) {
 304           fScalars->Fill(3);
 305           selV0s.Add(particle);
 306         }
 307         break;
 308       }
 309     case kXiMinus:
 310     case kOmegaMinus:
 311       {
 312         if (particle->Pt() > cutPtCascade) {
 313           fScalars->Fill(6);
 314           selCascades.Add(particle);
 315         }
 316         break;
 317       }
 318     default: break;
 319     }
 320   }
 321
 322   // ESD information
 323   AliVEvent* event = InputEvent();
 324   if (!event) {
 325     Printf("ERROR: Could not retrieve event");
 326     return;
 327   }
 328
 329   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 330   // loop over tracks
 331   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 332     AliESDtrack* track = esd->GetTrack(iTrack);
 333
 334     // select tracks of selected particles
 335     Int_t label = TMath::Abs(track->GetLabel());
 336     if (label > mcEvent->GetNumberOfTracks()) continue;     // background
 337     AliMCParticle* mctrack = (AliMCParticle*)mcEvent->GetTrack(label);
 338     TParticle* particle = mctrack->Particle();
 339     if (!selParticles.Contains(particle)) continue;
 340     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
 341     if (track->GetConstrainedChi2() > 1e9) continue;
 342     selParticles.Remove(particle);   // don't count multiple tracks
 343
 344     fScalars->Fill(1);
 345     fRec->Fill(particle->Pt());
 346     if (track->GetLabel() < 0) {
 347       fScalars->Fill(2);
 348     }
 349
 350     // resolutions
 351     fResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) *particle->Pt());
 352     fResPhi->Fill(1000. * (track->Phi() - particle->Phi()));
 353     fResTheta->Fill(1000. * (track->Theta() - particle->Theta()));
 354
 355     // PID
 356     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
 357     Int_t iGen = 5;
 358
 359     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 360       if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i;
 361     }
 362     Double_t probability[AliPID::kSPECIES];
 363     track->GetESDpid(probability);
 364     Double_t pMax = 0;
 365     Int_t iRec = 0;
 366     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 367       probability[i] *= partFrac[i];
 368       if (probability[i] > pMax) {
 369         pMax = probability[i];
 370         iRec = i;
 371       }
 372
 373     }
 374     fArrayHist->Fill(AliPID::kSPECIES*iGen + iRec);
 375     if (iGen == iRec) {
 376       fScalars->Fill(5);
 377     }
 378
 379     // dE/dx and TOF
 380     Double_t time[AliPID::kSPECIESC];
 381     track->GetIntegratedTimes(time,AliPID::kSPECIESC);
 382     if (iGen == iRec) {
 383       fDEdxRight->Fill(particle->P(), track->GetTPCsignal());
 384       if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {
 385         fResTOFRight->Fill(track->GetTOFsignal() - time[iRec]);
 386       }
 387     }
 388     else {
 389       fDEdxWrong->Fill(particle->P(), track->GetTPCsignal());
 390       if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) {
 391         fResTOFWrong->Fill(track->GetTOFsignal() - time[iRec]);
 392       }
 393     }
 394   }
 395
 396   // loop over muon tracks
 397   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
 398     AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
 399     Double_t ptInv = TMath::Abs(muonTrack->GetInverseBendingMomentum());
 400     if (ptInv > 0.001) {
 401       fPtMUON->Fill(1./ptInv);
 402     }
 403   }
 404
 405   // loop over V0s
 406   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
 407     AliESDv0* v0 = esd->GetV0(iV0);
 408     if (v0->GetOnFlyStatus()) continue;
 409     v0->ChangeMassHypothesis(kK0Short);
 410     fMassK0->Fill(v0->GetEffMass());
 411     v0->ChangeMassHypothesis(kLambda0);
 412     fMassLambda->Fill(v0->GetEffMass());
 413     v0->ChangeMassHypothesis(kLambda0Bar);
 414     fMassLambdaBar->Fill(v0->GetEffMass());
 415
 416     Int_t negLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel());
 417     if (negLabel > mcEvent->GetNumberOfTracks()) continue;     // background
 418     AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(negLabel);
 419     TParticle* particle = track->Particle();
 420     Int_t negMother = particle->GetMother(0);
 421     if (negMother < 0) continue;
 422     Int_t posLabel = TMath::Abs(esd->GetTrack(v0->GetPindex())->GetLabel());
 423     if (posLabel > mcEvent->GetNumberOfTracks()) continue;     // background
 424     track = (AliMCParticle*)mcEvent->GetTrack(posLabel);
 425     particle = track->Particle();
 426     Int_t posMother = particle->GetMother(0);
 427     if (negMother != posMother) continue;
 428     track = (AliMCParticle*)mcEvent->GetTrack(negMother);
 429     particle = track->Particle();
 430     if (!selV0s.Contains(particle)) continue;
 431     selV0s.Remove(particle);
 432     fScalars->Fill(4);
 433   }
 434
 435   // loop over Cascades
 436   for (Int_t iCascade = 0; iCascade < esd->GetNumberOfCascades(); iCascade++) {
 437     AliESDcascade* cascade = esd->GetCascade(iCascade);
 438     Double_t v0q;
 439     cascade->ChangeMassHypothesis(v0q,kXiMinus);
 440     fMassXi->Fill(cascade->GetEffMassXi());
 441     cascade->ChangeMassHypothesis(v0q,kOmegaMinus);
 442     fMassOmega->Fill(cascade->GetEffMassXi());
 443
 444     Int_t negLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex())->GetLabel());
 445     if (negLabel > mcEvent->GetNumberOfTracks()) continue;     // background
 446     AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(negLabel);
 447     TParticle* particle = track->Particle();
 448     Int_t negMother = particle->GetMother(0);
 449     if (negMother < 0) continue;
 450     Int_t posLabel = TMath::Abs(esd->GetTrack(cascade->GetPindex())->GetLabel());
 451     if (posLabel > mcEvent->GetNumberOfTracks()) continue;     // background
 452     track = (AliMCParticle*)mcEvent->GetTrack(posLabel);
 453     particle = track->Particle();
 454     Int_t posMother = particle->GetMother(0);
 455     if (negMother != posMother) continue;
 456     track = (AliMCParticle*)mcEvent->GetTrack(negMother);
 457     particle = track->Particle();
 458     Int_t v0Mother = particle->GetMother(0);
 459     if (v0Mother < 0) continue;
 460     Int_t bacLabel = TMath::Abs(esd->GetTrack(cascade->GetBindex())->GetLabel());
 461     if (bacLabel > mcEvent->GetNumberOfTracks()) continue;     // background
 462     track = (AliMCParticle*)mcEvent->GetTrack(bacLabel);
 463     particle = track->Particle();
 464     Int_t bacMother = particle->GetMother(0);
 465     if (v0Mother != bacMother) continue;
 466     track = (AliMCParticle*)mcEvent->GetTrack(v0Mother);
 467     particle = track->Particle();
 468     if (!selCascades.Contains(particle)) continue;
 469     selCascades.Remove(particle);
 470     fScalars->Fill(7);
 471   }
 472
 473   // loop over the clusters
 474   for (Int_t iCluster=0; iCluster<esd->GetNumberOfCaloClusters(); iCluster++) {
 475     AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster);
 476     if (clust->IsPHOS()) fEPHOS->Fill(clust->E());
 477     if (clust->IsEMCAL()) fEEMCAL->Fill(clust->E());
 478   }
 479
 480   // Post output data.
 481   PostData(1, fListOfHistos);
 482 }
 483
 484 void AliAnalysisTaskCheckESD::Terminate(Option_t *)
 485 {
 486   // check values
 487   Int_t    checkNGenLow = 1;
 488   Double_t checkEffLow = 0.5;
 489   Double_t checkEffSigma = 3;
 490   Double_t checkFakeHigh = 0.5;
 491   Double_t checkFakeSigma = 3;
 492
 493   Double_t checkResPtInvHigh = 5;
 494   Double_t checkResPtInvSigma = 3;
 495   Double_t checkResPhiHigh = 10;
 496   Double_t checkResPhiSigma = 3;
 497   Double_t checkResThetaHigh = 10;
 498   Double_t checkResThetaSigma = 3;
 499
 500   Double_t checkPIDEffLow = 0.5;
 501   Double_t checkPIDEffSigma = 3;
 502   Double_t checkResTOFHigh = 500;
 503   Double_t checkResTOFSigma = 3;
 504
 505   Double_t checkPHOSNLow = 5;
 506   Double_t checkPHOSEnergyLow = 0.3;
 507   Double_t checkPHOSEnergyHigh = 1.0;
 508   Double_t checkEMCALNLow = 50;
 509   Double_t checkEMCALEnergyLow = 0.05;
 510   Double_t checkEMCALEnergyHigh = 1.0;
 511
 512   Double_t checkMUONNLow = 1;
 513   Double_t checkMUONPtLow = 0.5;
 514   Double_t checkMUONPtHigh = 10.;
 515
 516   Double_t checkV0EffLow = 0.02;
 517   Double_t checkV0EffSigma = 3;
 518   Double_t checkCascadeEffLow = 0.01;
 519   Double_t checkCascadeEffSigma = 3;
 520
 521   const char* partName[AliPID::kSPECIES+1] = {"electron", "muon", "pion", "kaon", "proton", "other"};
 522   fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
 523   if (!fListOfHistos) {
 524     Printf("ERROR: fListOfHistos not available");
 525     return;
 526   }
 527
 528   fGen = dynamic_cast<TH1F*>(fListOfHistos->At(0));
 529   fRec = dynamic_cast<TH1F*>(fListOfHistos->At(1));
 530   fResPtInv = dynamic_cast<TH1F*>(fListOfHistos->At(2));
 531   fResPhi = dynamic_cast<TH1F*>(fListOfHistos->At(3));
 532   fResTheta = dynamic_cast<TH1F*>(fListOfHistos->At(4));
 533   fDEdxRight = dynamic_cast<TH2F*>(fListOfHistos->At(5));
 534   fDEdxWrong = dynamic_cast<TH2F*>(fListOfHistos->At(6));
 535   fResTOFRight = dynamic_cast<TH1F*>(fListOfHistos->At(7));
 536   fResTOFWrong = dynamic_cast<TH1F*>(fListOfHistos->At(8));
 537   fEPHOS = dynamic_cast<TH1F*>(fListOfHistos->At(9));
 538   fEEMCAL = dynamic_cast<TH1F*>(fListOfHistos->At(10));
 539   fPtMUON = dynamic_cast<TH1F*>(fListOfHistos->At(11));
 540   fMassK0 = dynamic_cast<TH1F*>(fListOfHistos->At(12));
 541   fMassLambda = dynamic_cast<TH1F*>(fListOfHistos->At(13));
 542   fMassLambdaBar = dynamic_cast<TH1F*>(fListOfHistos->At(14));
 543   fMassXi = dynamic_cast<TH1F*>(fListOfHistos->At(15));
 544   fMassOmega = dynamic_cast<TH1F*>(fListOfHistos->At(16));
 545   fScalars = dynamic_cast<TH1F*>(fListOfHistos->At(17));
 546   fArrayHist = dynamic_cast<TH1F*>(fListOfHistos->At(18));
 547
 548   Int_t nGen = Int_t(fScalars->GetBinContent(1));
 549   Int_t nRec = Int_t(fScalars->GetBinContent(2));
 550   Int_t nFake = Int_t(fScalars->GetBinContent(3));
 551   Int_t nGenV0s = Int_t(fScalars->GetBinContent(4));
 552   Int_t nRecV0s = Int_t(fScalars->GetBinContent(5));
 553   Int_t nIdentified = Int_t(fScalars->GetBinContent(6));
 554   Int_t nGenCascades = Int_t(fScalars->GetBinContent(7));
 555   Int_t nRecCascades = Int_t(fScalars->GetBinContent(8));
 556
 557   Int_t k = 1;
 558
 559   Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES];
 560   for(Int_t i = 0; i < (AliPID::kSPECIES+1); i++)
 561     for(Int_t j = 0; j < AliPID::kSPECIES; j++) {
 562       identified[i][j] = Int_t(fArrayHist->GetBinContent(k));
 563       k++;
 564     }
 565
 566
 567   // perform checks
 568   if (nGen < checkNGenLow) {
 569     Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen));
 570   }
 571
 572   TH1F* hEff = CreateEffHisto(fGen, fRec);
 573
 574   Info("CheckESD", "%d out of %d tracks reconstructed including %d "
 575        "fake tracks", nRec, nGen, nFake);
 576   if (nGen > 0) {
 577     // efficiency
 578     Double_t eff = nRec*1./nGen;
 579     Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGen);
 580     Double_t fake = nFake*1./nGen;
 581     Double_t fakeError = TMath::Sqrt(fake*(1.-fake) / nGen);
 582     Info("CheckESD", "eff = (%.1f +- %.1f) %%  fake = (%.1f +- %.1f) %%",
 583          100.*eff, 100.*effError, 100.*fake, 100.*fakeError);
 584     if (eff < checkEffLow - checkEffSigma*effError) {
 585       Warning("CheckESD", "low efficiency: (%.1f +- %.1f) %%",100.*eff, 100.*effError);
 586     }
 587     if (fake > checkFakeHigh + checkFakeSigma*fakeError) {
 588       Warning("CheckESD", "high fake: (%.1f +- %.1f) %%",100.*fake, 100.*fakeError);
 589     }
 590     // resolutions
 591     Double_t res, resError;
 592     if (FitHisto(fResPtInv, res, resError)) {
 593       Info("CheckESD", "relative inverse pt resolution = (%.1f +- %.1f) %%",res, resError);
 594       if (res > checkResPtInvHigh + checkResPtInvSigma*resError) {
 595         Warning("CheckESD", "bad pt resolution: (%.1f +- %.1f) %%",res, resError);
 596       }
 597     }
 598     if (FitHisto(fResPhi, res, resError)) {
 599       Info("CheckESD", "phi resolution = (%.1f +- %.1f) mrad", res, resError);
 600       if (res > checkResPhiHigh + checkResPhiSigma*resError) {
 601         Warning("CheckESD", "bad phi resolution: (%.1f +- %.1f) mrad",
 602                 res, resError);
 603       }
 604     }
 605
 606     if (FitHisto(fResTheta, res, resError)) {
 607       Info("CheckESD", "theta resolution = (%.1f +- %.1f) mrad",
 608            res, resError);
 609       if (res > checkResThetaHigh + checkResThetaSigma*resError) {
 610         Warning("CheckESD", "bad theta resolution: (%.1f +- %.1f) mrad",
 611                 res, resError);
 612       }
 613     }
 614     // PID
 615     if (nRec > 0) {
 616       Double_t eff = nIdentified*1./nRec;
 617       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nRec);
 618       Info("CheckESD", "PID eff = (%.1f +- %.1f) %%",
 619            100.*eff, 100.*effError);
 620       if (eff < checkPIDEffLow - checkPIDEffSigma*effError) {
 621         Warning("CheckESD", "low PID efficiency: (%.1f +- %.1f) %%",
 622                 100.*eff, 100.*effError);
 623       }
 624     }
 625
 626     printf("%9s:", "gen\\rec");
 627     for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {
 628       printf("%9s", partName[iRec]);
 629     }
 630     printf("\n");
 631     for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) {
 632       printf("%9s:", partName[iGen]);
 633       for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) {
 634         printf("%9d", identified[iGen][iRec]);
 635       }
 636       printf("\n");
 637     }
 638
 639     if (FitHisto(fResTOFRight, res, resError)) {
 640       Info("CheckESD", "TOF resolution = (%.1f +- %.1f) ps", res, resError);
 641       if (res > checkResTOFHigh + checkResTOFSigma*resError) {
 642         Warning("CheckESD", "bad TOF resolution: (%.1f +- %.1f) ps",
 643                 res, resError);
 644       }
 645     }
 646
 647     // calorimeters
 648     if (fEPHOS->Integral() < checkPHOSNLow) {
 649       Warning("CheckESD", "low number of PHOS particles: %d",
 650               Int_t(fEPHOS->Integral()));
 651     }
 652     else {
 653       Double_t mean = fEPHOS->GetMean();
 654       if (mean < checkPHOSEnergyLow) {
 655         Warning("CheckESD", "low mean PHOS energy: %.1f GeV", mean);
 656       } else if (mean > checkPHOSEnergyHigh) {
 657         Warning("CheckESD", "high mean PHOS energy: %.1f GeV", mean);
 658       }
 659     }
 660
 661     if (fEEMCAL->Integral() < checkEMCALNLow) {
 662       Warning("CheckESD", "low number of EMCAL particles: %d",
 663               Int_t(fEEMCAL->Integral()));
 664     }
 665     else {
 666       Double_t mean = fEEMCAL->GetMean();
 667       if (mean < checkEMCALEnergyLow) {
 668         Warning("CheckESD", "low mean EMCAL energy: %.1f GeV", mean);
 669       }
 670       else if (mean > checkEMCALEnergyHigh) {
 671         Warning("CheckESD", "high mean EMCAL energy: %.1f GeV", mean);
 672       }
 673     }
 674
 675     // muons
 676     if (fPtMUON->Integral() < checkMUONNLow) {
 677       Warning("CheckESD", "low number of MUON particles: %d",
 678               Int_t(fPtMUON->Integral()));
 679     }
 680     else {
 681       Double_t mean = fPtMUON->GetMean();
 682       if (mean < checkMUONPtLow) {
 683         Warning("CheckESD", "low mean MUON pt: %.1f GeV/c", mean);
 684       }
 685       else if (mean > checkMUONPtHigh) {
 686         Warning("CheckESD", "high mean MUON pt: %.1f GeV/c", mean);
 687       }
 688     }
 689
 690     // V0s
 691     if (nGenV0s > 0) {
 692       Double_t eff = nRecV0s*1./nGenV0s;
 693       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenV0s);
 694       if (TMath::Abs(effError) < 1.e-33) effError = checkV0EffLow / TMath::Sqrt(1.*nGenV0s);
 695       Info("CheckESD", "V0 eff = (%.1f +- %.1f) %%",
 696            100.*eff, 100.*effError);
 697       if (eff < checkV0EffLow - checkV0EffSigma*effError) {
 698         Warning("CheckESD", "low V0 efficiency: (%.1f +- %.1f) %%",
 699                 100.*eff, 100.*effError);
 700       }
 701     }
 702
 703     // Cascades
 704     if (nGenCascades > 0) {
 705       Double_t eff = nRecCascades*1./nGenCascades;
 706       Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenCascades);
 707       if (effError == 0) effError = checkV0EffLow /
 708                            TMath::Sqrt(1.*nGenCascades);
 709       Info("CheckESD", "Cascade eff = (%.1f +- %.1f) %%",
 710            100.*eff, 100.*effError);
 711       if (eff < checkCascadeEffLow - checkCascadeEffSigma*effError) {
 712         Warning("CheckESD", "low Cascade efficiency: (%.1f +- %.1f) %%",
 713                 100.*eff, 100.*effError);
 714       }
 715     }
 716
 717   }
 718
 719   // draw the histograms if not in batch mode
 720   if (!gROOT->IsBatch()) {
 721     new TCanvas;
 722     hEff->DrawCopy();
 723     new TCanvas;
 724     fResPtInv->DrawCopy("E");
 725     new TCanvas;
 726     fResPhi->DrawCopy("E");
 727     new TCanvas;
 728     fResTheta->DrawCopy("E");
 729     new TCanvas;
 730     fDEdxRight->DrawCopy();
 731     fDEdxWrong->DrawCopy("SAME");
 732     new TCanvas;
 733     fResTOFRight->DrawCopy("E");
 734     fResTOFWrong->DrawCopy("SAME");
 735     new TCanvas;
 736     fEPHOS->DrawCopy("E");
 737     new TCanvas;
 738     fEEMCAL->DrawCopy("E");
 739     new TCanvas;
 740     fPtMUON->DrawCopy("E");
 741     new TCanvas;
 742     fMassK0->DrawCopy("E");
 743     new TCanvas;
 744     fMassLambda->DrawCopy("E");
 745     new TCanvas;
 746     fMassLambdaBar->DrawCopy("E");
 747     new TCanvas;
 748     fMassXi->DrawCopy("E");
 749     new TCanvas;
 750     fMassOmega->DrawCopy("E");
 751   }
 752
 753   // write the output histograms to a file
 754   TFile* outputFile = TFile::Open("check.root", "recreate");
 755   if (!outputFile || !outputFile->IsOpen())
 756     {
 757       Error("CheckESD", "opening output file check.root failed");
 758       return;
 759     }
 760   hEff->Write();
 761   fResPtInv->Write();
 762   fResPhi->Write();
 763   fResTheta->Write();
 764   fDEdxRight->Write();
 765   fDEdxWrong->Write();
 766   fResTOFRight->Write();
 767   fResTOFWrong->Write();
 768   fEPHOS->Write();
 769   fEEMCAL->Write();
 770   fPtMUON->Write();
 771   fMassK0->Write();
 772   fMassLambda->Write();
 773   fMassLambdaBar->Write();
 774   fMassXi->Write();
 775   fMassOmega->Write();
 776   outputFile->Close();
 777   delete outputFile;
 778
 779   // clean up
 780   /*    delete hGen;
 781         delete hRec;
 782         delete hEff;
 783         delete hResPtInv;
 784         delete hResPhi;
 785         delete hResTheta;
 786         delete hDEdxRight;
 787         delete hDEdxWrong;
 788         delete hResTOFRight;
 789         delete hResTOFWrong;
 790         delete hEPHOS;
 791         delete hEEMCAL;
 792         delete hPtMUON;
 793         delete hMassK0;
 794         delete hMassLambda;
 795         delete hMassLambdaBar;
 796         delete hMassXi;
 797         delete hMassOmega;
 798         delete hScalars;
 799         delete hArrayHist; */
 800
 801   //delete esd;
 802   // result of check
 803   Info("CheckESD", "check of ESD was successfull");
 804 }
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866