PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioHelper.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: ALICE Offline.                                                 *
   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 //=========================================================================//
  17 //             AliEbyE Analysis for Particle Ratio Fluctuation             //
  18 //                   Deepika Rathee  | Satyajit Jena                       //
  19 //                   drathee@cern.ch | sjena@cern.ch                       //
  20 //                  Date: Wed Jul  9 18:38:30 CEST 2014                    //
  21 //          New approch to find particle ratio to reduce memory            //
  22 //                             (Test Only)                                 //
  23 //=========================================================================//
  24
  25 #include "TMath.h"
  26 #include "TAxis.h"
  27 #include "TSystem.h"
  28 #include "TFile.h"
  29 #include "TPRegexp.h"
  30
  31 #include "AliStack.h"
  32 #include "AliMCEvent.h"
  33 #include "AliMCParticle.h"
  34 #include "AliESDtrackCuts.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliESDpid.h"
  37 #include "AliAODInputHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODMCParticle.h"
  40 #include "AliCentrality.h"
  41 #include "AliTracker.h"
  42
  43 #include "AliEbyEPidRatioHelper.h"
  44
  45 using namespace std;
  46
  47
  48 ClassImp(AliEbyEPidRatioHelper)
  49 //________________________________________________________________________
  50 AliEbyEPidRatioHelper::AliEbyEPidRatioHelper() :
  51   fModeDistCreation(0),
  52
  53   fInputEventHandler(NULL),
  54   fPIDResponse(NULL),
  55   fESD(NULL),
  56   fESDTrackCuts(NULL),
  57   fAOD(NULL),
  58   fAODtrackCutBit(1024),
  59   fIsMC(kFALSE),
  60   fMCEvent(NULL),
  61   fStack(NULL),
  62
  63   fCentralityBin(-1),
  64   fCentralityPercentile(-1.),
  65
  66   fCentralityBinMax(11),
  67   fVertexZMax(10.),
  68   fRapidityMax(0.9),
  69   fPhiMin(0.),
  70   fPhiMax(TMath::TwoPi()),
  71   fMinTrackLengthMC(70.),
  72   fNSigmaMaxCdd(3.),
  73   fNSigmaMaxCzz(3.),
  74
  75   fPIDStrategy(0),
  76   fNSigmaMaxITS(4.),
  77   fNSigmaMaxTPC(4.),
  78   fNSigmaMaxTPClow(3.),
  79   fNSigmaMaxTOF(4.),
  80   fMinPtForTOFRequired(0.69),
  81   fMaxPtForTPClow(0.69),
  82
  83   fHEventStat0(NULL),
  84   fHEventStat1(NULL),
  85   fHEventStatMax(6),
  86
  87   fHTriggerStat(NULL),
  88   fNTriggers(5),
  89
  90   fHCentralityStat(NULL),
  91   fNCentralityBins(11),
  92
  93   fRandom(NULL),
  94   fIsRatio(kFALSE),
  95   fIsPtBin(kFALSE) {
  96   // Constructor
  97
  98   AliLog::SetClassDebugLevel("AliEbyEPidRatioHelper",10);
  99 }
 100
 101 const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthRap = 0.1;
 102 const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthPt  = 0.3; // 0.08 // 300 MeV  // was 80 MeV
 103
 104 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeCent[]  = {-0.5, 10.5};
 105 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsCent    = 11 ;
 106
 107 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeEta[]   = {-0.9, 0.9};
 108 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsEta     = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeEta[1] -
 109                                                                    AliEbyEPidRatioHelper::fgkfHistRangeEta[0]) /
 110                                                                   AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
 111
 112 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeRap[]   = {-0.8, 0.8};
 113 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsRap     = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeRap[1] - AliEbyEPidRatioHelper::fgkfHistRangeRap[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
 114
 115 const Float_t AliEbyEPidRatioHelper::fgkfHistRangePhi[]   = {0.0, TMath::TwoPi()};
 116 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsPhi     = 21 ;
 117
 118 const Float_t AliEbyEPidRatioHelper::fgkfHistRangePt[]    = {0.2, 2.9}; // {0.2, 5.}; // was {0.3, 2.22}
 119 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsPt      = Int_t((AliEbyEPidRatioHelper::fgkfHistRangePt[1] - AliEbyEPidRatioHelper::fgkfHistRangePt[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthPt);
 120
 121 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeSign[]  = {-1.5, 1.5};
 122 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsSign    =  3;
 123
 124 const Char_t* AliEbyEPidRatioHelper::fgkEventNames[]         = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,100]%"};
 125 const Char_t* AliEbyEPidRatioHelper::fgkCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
 126 const Char_t* AliEbyEPidRatioHelper::fgkTriggerNames[]       = {"kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
 127 const Char_t* AliEbyEPidRatioHelper::fgkCentralityNames[]    = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%","50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
 128
 129 const Char_t* AliEbyEPidRatioHelper::fgkPidName[4]      = {"Nch","Npi","Nka","Npr"};
 130 const Char_t* AliEbyEPidRatioHelper::fgkPidLatex[4][2]  = {{"N_{-}","N_{+}"}, {"N_{#pi^{-}}","N_{#pi^{+}}"},{"N_{K^{-}}","N_{K^{+}}"}, {"N_{#bar{p}}","N_{p}"}};
 131 const Char_t* AliEbyEPidRatioHelper::fgkPidTitles[4][2] = {{"Negative","Positive"},{"Anti-Pions","Pions"},{"Anti-Kaons","Kaons"}, {"Anti-Protons","Protons"}};
 132
 133
 134
 135 //________________________________________________________________________
 136 AliEbyEPidRatioHelper::~AliEbyEPidRatioHelper() {
 137   // Destructor
 138
 139   if (fRandom)
 140     delete fRandom;
 141
 142   return;
 143 }
 144
 145 //________________________________________________________________________
 146 void AliEbyEPidRatioHelper::SetPhiRange(Float_t f1, Float_t f2) {
 147   // -- Set phi range and adopt to phi-histogram
 148
 149   fPhiMin = f1;
 150   fPhiMax = (f1 < f2) ? f2 : f2+TMath::TwoPi();
 151
 152   Float_t phiMin = fPhiMin;
 153   Float_t phiMax = fPhiMax;
 154
 155   // -- Update Ranges
 156   Float_t binWidth = (AliEbyEPidRatioHelper::fgkfHistRangePhi[1] - AliEbyEPidRatioHelper::fgkfHistRangePhi[0]) /
 157     Float_t(AliEbyEPidRatioHelper::fgkfHistNBinsPhi);
 158
 159   Float_t lowEdge  = AliEbyEPidRatioHelper::fgkfHistRangePhi[0] - binWidth;
 160   Float_t highEdge = AliEbyEPidRatioHelper::fgkfHistRangePhi[0];
 161
 162   for (Int_t ii = 1; ii <= AliEbyEPidRatioHelper::fgkfHistNBinsPhi; ++ii) {
 163     lowEdge += binWidth;
 164     highEdge += binWidth;
 165
 166     if (phiMin >= lowEdge && phiMin < highEdge )
 167       phiMin = lowEdge;
 168     if (phiMax > lowEdge && phiMax <= highEdge )
 169       phiMax = highEdge;
 170   }
 171
 172   printf(">>>> Update Phi Range : [%f,%f] -> [%f,%f]\n", fPhiMin, fPhiMax, phiMin, phiMax);
 173   fPhiMin = phiMin;
 174   fPhiMax = phiMax;
 175 }
 176
 177
 178
 179 //________________________________________________________________________
 180 Int_t AliEbyEPidRatioHelper::Initialize(AliESDtrackCuts *cuts, Bool_t isMC, Bool_t isRatio, Bool_t isPtBin, Int_t trackCutBit, Int_t modeDistCreation) {
 181   // -- Initialize helper
 182
 183   Int_t iResult = 0;
 184
 185   // -- ESD track cuts
 186   fESDTrackCuts     = cuts;
 187
 188
 189   // -- Is MC
 190   fIsMC             = isMC;
 191   fIsRatio          = isRatio;
 192   fIsPtBin          = isPtBin;
 193
 194
 195
 196
 197   // -- AOD track filter bit
 198   fAODtrackCutBit   = trackCutBit;
 199
 200   // -- mode Distribution creation
 201   fModeDistCreation = modeDistCreation;
 202
 203   // -- Setup event cut statistics
 204   InitializeEventStats();
 205
 206   // -- Setup trigger statistics
 207   InitializeTriggerStats();
 208
 209   // -- Setup centrality statistics
 210   InitializeCentralityStats();
 211
 212   // -- PRINT PID Strategy
 213   //    0 :   TPC(TPClow+TPCHigh)
 214   //    1 :   ITS
 215   //    2 :   TOF
 216   //    3 :   ITS+TPC(TPClow+TPCHigh)
 217   //    4 :   TPC(TPClow+TPCHigh)+TOF
 218   //    5 :   TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 219   //    6 :   TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
 220   //    7 :   TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 221   //    8 :   TPC(TPClow+TPCHigh)+ITS+TOF
 222   printf(">>>>  PID STRATEGY: %d || sigmaMax: ITS %.2f TPC %.2f TOF %.2f \n", fPIDStrategy, fNSigmaMaxITS, fNSigmaMaxTPC, fNSigmaMaxTOF);
 223
 224   // -- Initialize random number generator
 225   fRandom = new TRandom3();
 226   fRandom->SetSeed();
 227
 228   return iResult;
 229 }
 230
 231 //________________________________________________________________________
 232 Int_t AliEbyEPidRatioHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
 233   // -- Setup Event
 234
 235   // -- Get ESD objects
 236   if(esdHandler){
 237     fInputEventHandler = static_cast<AliInputEventHandler*>(esdHandler);
 238     fESD               = dynamic_cast<AliESDEvent*>(fInputEventHandler->GetEvent());
 239     if (!fESD) {
 240       AliError("ESD event handler not available");
 241       return -1;
 242     }
 243   }
 244
 245   // -- Get AOD objects
 246   else if(aodHandler){
 247     fInputEventHandler = static_cast<AliInputEventHandler*>(aodHandler);
 248     fAOD               = dynamic_cast<AliAODEvent*>(fInputEventHandler->GetEvent());
 249     if (!fAOD) {
 250       AliError("AOD event handler not available");
 251       return -1;
 252     }
 253   }
 254
 255   // -- Get Common objects
 256   fPIDResponse = fInputEventHandler->GetPIDResponse();
 257
 258   // -- Get MC objects
 259   fMCEvent     = mcEvent;
 260   if (fMCEvent)
 261     fStack     = fMCEvent->Stack();
 262
 263   // -- Get event centrality
 264   // >  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90|90-100 --> 11 bins
 265   // >   0   1     2     3     4     5     6     7     8     9     10
 266
 267   AliCentrality *centrality = NULL;
 268
 269   if(esdHandler)
 270     centrality = fESD->GetCentrality();
 271   else if(aodHandler)
 272     centrality = fAOD->GetHeader()->GetCentralityP();
 273
 274   if (!centrality) {
 275     AliError("Centrality not available");
 276     return -1;
 277   }
 278
 279   Int_t centBin = centrality->GetCentralityClass10("V0M");
 280   if (centBin == 0) { fCentralityBin = centrality->GetCentralityClass5("V0M"); }
 281   else if (centBin == 11 || centBin == -1.)           { fCentralityBin = -1; }
 282   else if (centBin > 0 && centBin < fNCentralityBins) { fCentralityBin = centBin + 1; }
 283   else {  fCentralityBin = -2; }
 284
 285   if (fCentralityBin >= fCentralityBinMax)
 286     fCentralityBin = -2;
 287
 288   fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
 289
 290   return 0;
 291 }
 292 //________________________________________________________________________
 293 Bool_t AliEbyEPidRatioHelper::IsEventTriggered() {
 294   // -- Check if Event is triggered and fill Trigger Histogram
 295
 296   Bool_t *aTriggerFired = new Bool_t[fNTriggers];
 297   for (Int_t ii = 0; ii < fNTriggers; ++ii)
 298     aTriggerFired[ii] = kFALSE;
 299
 300   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kMB))          aTriggerFired[0] = kTRUE;
 301   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kCentral))     aTriggerFired[1] = kTRUE;
 302   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
 303   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      aTriggerFired[3] = kTRUE;
 304   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      aTriggerFired[4] = kTRUE;
 305
 306   Bool_t isTriggered = kFALSE;
 307
 308   for (Int_t ii=0; ii<fNTriggers; ++ii) {
 309     if(aTriggerFired[ii]) {
 310       isTriggered = kTRUE;
 311       fHTriggerStat->Fill(ii);
 312     }
 313   }
 314
 315   delete[] aTriggerFired;
 316
 317   return isTriggered;
 318 }
 319
 320 //________________________________________________________________________
 321 Bool_t AliEbyEPidRatioHelper::IsEventRejected() {
 322   // -- Evaluate event statistics histograms
 323
 324   Int_t *aEventCuts = new Int_t[fHEventStatMax];
 325   // set aEventCuts[ii] to 1 in case of reject
 326
 327   for (Int_t ii=0;ii<fHEventStatMax; ++ii)
 328     aEventCuts[ii] = 0;
 329
 330   Int_t iCut = 0;
 331
 332   // -- 0 - Before Physics Selection
 333   aEventCuts[iCut] = 0;
 334
 335   // -- 1 - No Trigger fired
 336   ++iCut;
 337   if (!IsEventTriggered())
 338     aEventCuts[iCut] = 1;
 339
 340   // -- 2 - No Vertex
 341   ++iCut;
 342   const AliESDVertex* vtxESD = NULL;
 343   const AliAODVertex* vtxAOD = NULL;
 344   if (fESD){
 345     vtxESD = fESD->GetPrimaryVertexTracks();
 346     if (!vtxESD)
 347       aEventCuts[iCut] = 1;
 348   }
 349   else if (fAOD){
 350     vtxAOD = fAOD->GetPrimaryVertex();
 351     if (!vtxAOD)
 352       aEventCuts[iCut] = 1;
 353   }
 354
 355   // -- 3 - Vertex z outside cut window
 356   ++iCut;
 357   if (vtxESD){
 358     if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax)
 359       aEventCuts[iCut] = 1;
 360   }
 361   else if(vtxAOD){
 362     if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax)
 363       aEventCuts[iCut] = 1;
 364   }
 365   else
 366     aEventCuts[iCut] = 1;
 367
 368   // -- 4 - Centrality = -1  (no centrality or not hadronic)
 369   ++iCut;
 370   if(fCentralityBin == -1.)
 371     aEventCuts[iCut] = 1;
 372
 373   // -- 5 - Centrality < fCentralityMax
 374   ++iCut;
 375   if(fCentralityBin == -2.)
 376     aEventCuts[iCut] = 1;
 377
 378   // -- Fill statistics / reject event
 379   Bool_t isRejected = FillEventStats(aEventCuts);
 380
 381   // -- Cleanup
 382   delete[] aEventCuts;
 383
 384   //cout << isRejected << endl;
 385
 386   return isRejected;
 387 }
 388
 389 //________________________________________________________________________
 390 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC) {
 391   // -- Check if MC particle is accepted for basic parameters
 392
 393   if (!particle)
 394     return kFALSE;
 395
 396   // -- check if charged
 397   if (particle->Charge() == 0.0)
 398     return kFALSE;
 399
 400   // -- check if physical primary - ESD
 401   if (fESD) {
 402     if(!fStack->IsPhysicalPrimary(idxMC))
 403       return kFALSE;
 404   }
 405   // -- check if physical primary - AOD
 406   else {
 407     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
 408       return kFALSE;
 409   }
 410
 411   return kTRUE;
 412 }
 413
 414 //________________________________________________________________________
 415 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC) {
 416   // -- Check if MC particle is accepted for basic parameters
 417
 418   if (!particle)
 419     return kFALSE;
 420
 421   // -- check if charged
 422   if (particle->Charge() != 0.0)
 423     return kFALSE;
 424
 425   // -- check if physical primary - ESD
 426   if (fESD) {
 427     if(!fStack->IsPhysicalPrimary(idxMC))
 428       return kFALSE;
 429   }
 430   // -- check if physical primary - AOD
 431   else {
 432     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
 433       return kFALSE;
 434   }
 435
 436   return kTRUE;
 437 }
 438
 439 //________________________________________________________________________
 440 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP, Int_t gCurPid) {
 441
 442  if (gCurPid == 0) {
 443    yP = particle->Eta();
 444    return kTRUE;
 445  }
 446
 447  Double_t mP = AliPID::ParticleMass(AliPID::kPion);
 448  if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
 449  else if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
 450  else if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
 451
 452  // -- Calculate rapidities and kinematics
 453   Double_t p  = particle->P();
 454   Double_t pz = particle->Pz();
 455
 456   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 457   yP          = 0.5 * TMath::Log((eP + pz) / (eP - pz));
 458
 459   // -- Check Rapidity window
 460   if (TMath::Abs(yP) > fRapidityMax)
 461     return kFALSE;
 462
 463   return kTRUE;
 464 }
 465
 466 //________________________________________________________________________
 467 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedPhi(AliVParticle *particle) {
 468   // -- Check if particle is accepted
 469   // > in phi
 470   // > return 0 if not accepted
 471
 472   if (particle->Phi() > fPhiMin && particle->Phi() <= fPhiMax)
 473     return kTRUE;
 474   else if (particle->Phi() < fPhiMin && (particle->Phi() + TMath::TwoPi()) <= fPhiMax)
 475     return kTRUE;
 476   else
 477     return kFALSE;
 478 }
 479
 480 //_____________________________________________________________________________
 481 Bool_t AliEbyEPidRatioHelper::IsParticleFindable(Int_t label) {
 482   // -- Check if MC particle is findable tracks
 483
 484   AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
 485   if(!mcParticle)
 486     return kFALSE;
 487
 488   Int_t counter;
 489   Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
 490
 491   return (tpcTrackLength > fMinTrackLengthMC);
 492 }
 493 //________________________________________________________________________
 494 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedBasicCharged(AliVTrack* track) {
 495   // -- Check if track is accepted
 496   // > for basic parameters
 497
 498   if (!track)
 499     return kFALSE;
 500
 501   if (track->Charge() == 0)
 502     return kFALSE;
 503
 504   return kTRUE;
 505 }
 506
 507 //________________________________________________________________________
 508 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP, Int_t gCurPid) {
 509    if (gCurPid == 0) {
 510     yP = track->Eta();
 511     return kTRUE;
 512   }
 513
 514   Double_t mP = AliPID::ParticleMass(AliPID::kPion);
 515   if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
 516   else if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
 517   else if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
 518
 519   // -- Calculate rapidities and kinematics
 520   Double_t pvec[3];
 521   track->GetPxPyPz(pvec);
 522
 523   Double_t p  = track->P();
 524   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 525            yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
 526
 527   // -- Check Rapidity window
 528   if (TMath::Abs(yP) > fRapidityMax)
 529     return kFALSE;
 530
 531   return kTRUE;
 532 }
 533
 534 //________________________________________________________________________
 535 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedDCA(AliVTrack *vTrack) {
 536    Bool_t isAccepted = kTRUE;
 537
 538   if (!fESD)
 539     return isAccepted;
 540
 541   AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
 542   if (!track)
 543     return kFALSE;
 544
 545   // -- Get nHits SPD
 546   if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
 547
 548     // -- Get DCA nSigmas
 549     Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 550     track->GetImpactParameters(dca,cov);
 551
 552     Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
 553     Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
 554
 555     if (fNSigmaMaxCdd != 0.) {
 556       if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
 557         isAccepted = kFALSE;
 558     }
 559
 560     if (fNSigmaMaxCzz != 0.) {
 561       if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
 562         isAccepted = kFALSE;
 563     }
 564   }
 565
 566   return isAccepted;
 567 }
 568
 569 //________________________________________________________________________
 570 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid, AliPID::EParticleType gCurPid) {
 571
 572   Bool_t isAcceptedITS    = kFALSE;
 573   Bool_t isAcceptedTPC    = kFALSE;
 574   Bool_t isAcceptedTPClow = kFALSE;
 575   Bool_t isAcceptedTOF    = kFALSE;
 576   Bool_t isAccepted       = kFALSE;
 577
 578   // -- In case not PID is used
 579   if (gCurPid == 0) {
 580     pid[0] = 10.;
 581     pid[1] = 10.;
 582     pid[2] = 10.;
 583     return kTRUE;
 584   }
 585
 586   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kITS, track, gCurPid, pid[0]) == AliPIDResponse::kDetPidOk) {
 587     if (TMath::Abs(pid[0]) < fNSigmaMaxITS)
 588       isAcceptedITS = kTRUE;
 589   }
 590
 591   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, gCurPid, pid[1]) == AliPIDResponse::kDetPidOk) {
 592     if (TMath::Abs(pid[1]) < fNSigmaMaxTPC)
 593       isAcceptedTPC = kTRUE;
 594     if (TMath::Abs(pid[1]) < fNSigmaMaxTPClow)
 595       isAcceptedTPClow = kTRUE;
 596     if (track->Pt() < fMaxPtForTPClow)
 597       isAcceptedTPC = isAcceptedTPClow;
 598   }
 599
 600   Bool_t hasPIDTOF = kFALSE;
 601   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, gCurPid, pid[2]) == AliPIDResponse::kDetPidOk) {
 602     hasPIDTOF = kTRUE;
 603     if (TMath::Abs(pid[2]) < fNSigmaMaxTOF)
 604       isAcceptedTOF = kTRUE;
 605   }
 606
 607   // -- Check TOF missmatch for MC
 608
 609   //if (ESD)
 610   if (fIsMC && isAcceptedTOF) {
 611     Int_t tofLabel[3];
 612     //    AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
 613     // TODO add code for AOD
 614
 615     (dynamic_cast<AliESDtrack*>(track))->GetTOFLabel(tofLabel);
 616
 617     Bool_t hasMatchTOF = kTRUE;
 618     if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0)
 619       hasMatchTOF = kFALSE;
 620
 621     TParticle *matchedTrack = fStack->Particle(TMath::Abs(tofLabel[0]));
 622     if (TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
 623       hasMatchTOF = kTRUE;
 624
 625     isAcceptedTOF = hasMatchTOF;
 626   }
 627
 628   //    0 :   TPC(TPClow+TPCHigh)
 629   //    1 :   ITS
 630   //    2 :   TOF
 631   //    3 :   ITS+TPC(TPClow+TPCHigh)
 632   //    4 :   TPC(TPClow+TPCHigh)+TOF
 633   //    5 :   TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 634   //    6 :   TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
 635   //    7 :   TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 636   //    8 :   TPC(TPClow+TPCHigh)+ITS+TOF
 637   if (fPIDStrategy == 0) {             //  TPC PID
 638     isAccepted = isAcceptedTPC;
 639   }
 640   else if (fPIDStrategy == 1) {        //  ITS PID
 641     isAccepted = isAcceptedITS;
 642   }
 643   else if (fPIDStrategy == 2) {        //  TOF PID
 644     isAccepted = isAcceptedTOF;
 645   }
 646   else if (fPIDStrategy == 3) {        //  TPC+ITS PID
 647     isAccepted = isAcceptedTPC && isAcceptedITS;
 648   }
 649   else if (fPIDStrategy == 4) {        //  TPC+TOF PID
 650     isAccepted = isAcceptedTPC && isAcceptedTOF;
 651   }
 652   else if (fPIDStrategy == 5) {        //  TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
 653     if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
 654       isAcceptedTOF = kTRUE;
 655     isAccepted = isAcceptedTPC && isAcceptedTOF;
 656   }
 657   else if (fPIDStrategy == 6) {        //  ITS+TPC+TOF PID -- TOF only for those tracks which have TOF information
 658     isAccepted = isAcceptedTPC && isAcceptedITS;
 659     if (hasPIDTOF)
 660       isAccepted = isAccepted && isAcceptedTOF;
 661   }
 662   else if (fPIDStrategy == 7) {        //  ITS+TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
 663     if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
 664       isAcceptedTOF = kTRUE;
 665     isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
 666   }
 667   else if (fPIDStrategy == 8) {        //  ITS+TPC+TOF PID
 668     isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
 669   }
 670
 671   return isAccepted;
 672 }
 673
 674 //________________________________________________________________________
 675 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedPhi(AliVTrack *track) {
 676   // -- Check if track is accepted
 677   // > in phi
 678   // > return 0 if not accepted
 679
 680   if (track->Phi() > fPhiMin && track->Phi() <= fPhiMax)
 681     return kTRUE;
 682   else if (track->Phi() < fPhiMin && (track->Phi() + TMath::TwoPi()) <= fPhiMax)
 683     return kTRUE;
 684   else
 685     return kFALSE;
 686 }
 687
 688 //________________________________________________________________________
 689 void AliEbyEPidRatioHelper::BinLogAxis(const THnBase *hn, Int_t axisNumber, AliESDtrackCuts* cuts) {
 690     AliESDtrackCuts* esdTrackCuts = (cuts) ? cuts : fESDTrackCuts;
 691
 692   // -- Make logarithmic binning
 693   TAxis *axis = hn->GetAxis(axisNumber);
 694   Int_t  nBins = axis->GetNbins();
 695
 696   Double_t from  = axis->GetXmin();
 697   Double_t to    = axis->GetXmax();
 698   Double_t *newBins = new Double_t[nBins + 1];
 699
 700   newBins[0] = from;
 701   Double_t factor = TMath::Power(to/from, 1./nBins);
 702
 703   for (int ii = 1; ii <= nBins; ii++)
 704    newBins[ii] = factor * newBins[ii-1];
 705
 706   axis->Set(nBins, newBins);
 707
 708   delete [] newBins;
 709
 710   // -- Update Ranges
 711   // ------------------
 712   Float_t ptRange[2];
 713   Float_t oldPtRange[2];
 714   esdTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 715   esdTrackCuts->GetPtRange(oldPtRange[0],oldPtRange[1]);
 716
 717   Float_t minPtForTOFRequired = fMinPtForTOFRequired;
 718   Float_t maxPtForTPClow      = fMaxPtForTPClow;
 719
 720   // -- Update minPt Cut
 721   Int_t bin = axis->FindBin(ptRange[0]+10e-6);
 722   ptRange[0] = axis->GetBinLowEdge(bin);
 723
 724   // -- Update maxPt Cut
 725   bin = axis->FindBin(ptRange[1]-10e-6);
 726   ptRange[1] = axis->GetBinUpEdge(bin);
 727
 728   if (ptRange[0] != oldPtRange[0] || ptRange[1] != oldPtRange[1]) {
 729     printf(">>>> Update Pt Range : [%f,%f] -> [%f,%f]\n", oldPtRange[0], oldPtRange[1], ptRange[0], ptRange[1]);
 730     esdTrackCuts->SetPtRange(ptRange[0],ptRange[1]);
 731   }
 732
 733   // -- Update MinPtForTOFRequired
 734   bin = axis->FindBin(minPtForTOFRequired-10e-6);
 735   minPtForTOFRequired = axis->GetBinUpEdge(bin);
 736
 737   if (minPtForTOFRequired != fMinPtForTOFRequired) {
 738     printf(">>>> Update Min Pt for TOF : %f -> %f\n", fMinPtForTOFRequired, minPtForTOFRequired);
 739     fMinPtForTOFRequired = minPtForTOFRequired;
 740   }
 741
 742   // -- Update MaxPtForTPClow
 743   bin = axis->FindBin(maxPtForTPClow-10e-6);
 744   maxPtForTPClow = axis->GetBinUpEdge(bin);
 745
 746   if (maxPtForTPClow != fMaxPtForTPClow) {
 747     printf(">>>> Update Max Pt for TPC Low : %f -> %f\n", fMaxPtForTPClow, maxPtForTPClow);
 748     fMaxPtForTPClow = maxPtForTPClow;
 749   }
 750 }
 751
 752 //________________________________________________________________________
 753 void AliEbyEPidRatioHelper::InitializeEventStats() {
 754   // -- Initialize event statistics histograms
 755
 756   fHEventStat0 = new TH1F("hEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 757   fHEventStat1 = new TH1F("hEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 758
 759   for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
 760     fHEventStat0->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkEventNames[ii]);
 761     fHEventStat1->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkEventNames[ii]);
 762   }
 763
 764   fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
 765   fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
 766 }
 767
 768 //________________________________________________________________________
 769 void AliEbyEPidRatioHelper::InitializeTriggerStats() {
 770   // -- Initialize trigger statistics histograms
 771
 772   fHTriggerStat = new TH1F("hTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
 773
 774   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 775     fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkTriggerNames[ii]);
 776 }
 777
 778 //________________________________________________________________________
 779 void AliEbyEPidRatioHelper::InitializeCentralityStats() {
 780   // -- Initialize trigger statistics histograms
 781
 782   fHCentralityStat = new TH1F("hCentralityStat","Centrality statistics;Centrality Bins;Events",
 783                               fNCentralityBins,-0.5,fNCentralityBins-0.5);
 784
 785   for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
 786     fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkCentralityNames[ii]);
 787 }
 788 //________________________________________________________________________
 789 Bool_t AliEbyEPidRatioHelper::FillEventStats(Int_t *aEventCuts) {
 790   // -- Fill event / centrality statistics
 791
 792   Bool_t isRejected = kFALSE;
 793
 794   // -- Fill event statistics
 795   for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
 796
 797     if (aEventCuts[idx])
 798       isRejected = kTRUE;
 799     else
 800       fHEventStat0->Fill(idx);
 801   }
 802
 803   for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
 804     if (aEventCuts[idx])
 805       break;
 806     fHEventStat1->Fill(idx);
 807   }
 808
 809   // -- Fill centrality statistics of accepted events
 810   if (!isRejected)
 811     fHCentralityStat->Fill(fCentralityBin);
 812
 813   return isRejected;
 814 }