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