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