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