PWGCF/EBYE/NetParticle/AliAnalysisNetParticleHelper.cxx

   1 //-*- Mode: C++ -*-
   2
   3 #include "TMath.h"
   4 #include "TAxis.h"
   5 #include "TSystem.h"
   6 #include "TFile.h"
   7 #include "TPRegexp.h"
   8
   9 #include "AliStack.h"
  10 #include "AliMCEvent.h"
  11 #include "AliMCParticle.h"
  12 #include "AliESDtrackCuts.h"
  13 #include "AliESDInputHandler.h"
  14 #include "AliESDpid.h"
  15 #include "AliAODInputHandler.h"
  16 #include "AliAODEvent.h"
  17 #include "AliAODMCParticle.h"
  18 #include "AliCentrality.h"
  19 #include "AliTracker.h"
  20
  21 #include "AliAnalysisNetParticleHelper.h"
  22
  23 using namespace std;
  24
  25 /**
  26  * Class for NetParticle Distributions
  27  * -- Helper class for net particle istributions
  28  * Authors: Jochen Thaeder <jochen@thaeder.de>
  29  *          Michael Weber <m.weber@cern.ch>
  30  */
  31
  32 ClassImp(AliAnalysisNetParticleHelper)
  33
  34 /*
  35  * ---------------------------------------------------------------------------------
  36  *                            particle names
  37  * ---------------------------------------------------------------------------------
  38  */
  39
  40   // MW make fgk ... static const
  41   const Char_t* aPartNames[AliPID::kSPECIES][2]             = {
  42     {"ele",     "posi"},
  43     {"mubar",   "mu"},
  44     {"pibar",   "pi"},
  45     {"kbar",    "k"},
  46     {"pbar",    "p"}
  47   };
  48
  49   const Char_t* aPartTitles[AliPID::kSPECIES][2]            = {
  50     {"Electron",    "Positron"},
  51     {"Anti-Muon",   "Muon"},
  52     {"Anti-Pion",   "Proton"},
  53     {"Anti-Kaon",   "Kaon"},
  54     {"Anti-Proton", "Proton"}
  55   };
  56
  57   const Char_t* aPartTitlesLatex[AliPID::kSPECIES][2]        = {
  58     {"e^{-}",   "e^{+}" },
  59     {"#mu^{-}", "#mu^{+}" },
  60     {"#pi^{-}", "#pi^{+}" },
  61     {"K^{-}",   "K^{+}" },
  62     {"#bar{p}", "p"}
  63   };
  64
  65 /*
  66  * ---------------------------------------------------------------------------------
  67  *                            Constructor / Destructor
  68  * ---------------------------------------------------------------------------------
  69  */
  70
  71 //________________________________________________________________________
  72 AliAnalysisNetParticleHelper::AliAnalysisNetParticleHelper() :
  73   fModeDistCreation(0),
  74
  75   fInputEventHandler(NULL),
  76   fPIDResponse(NULL),
  77   fESD(NULL),
  78   fESDTrackCuts(NULL),
  79   fAOD(NULL),
  80   fAODtrackCutBit(1024),
  81   fIsMC(kFALSE),
  82   fMCEvent(NULL),
  83   fStack(NULL),
  84
  85   fCentralityBin(-1),
  86   fCentralityPercentile(-1.),
  87
  88   fCentralityBinMax(7),
  89   fVertexZMax(10.),
  90   fRapidityMax(0.5),
  91   fPhiMin(0.),
  92   fPhiMax(TMath::TwoPi()),
  93   fMinTrackLengthMC(70.),
  94   fNSigmaMaxCdd(3.),
  95   fNSigmaMaxCzz(3.),
  96
  97   fParticleSpecies(AliPID::kProton),
  98
  99   fUsePID(kTRUE),
 100   fPIDStrategy(0),
 101   fNSigmaMaxITS(4.),
 102   fNSigmaMaxTPC(4.),
 103   fNSigmaMaxTPClow(3.),
 104   fNSigmaMaxTOF(4.),
 105   fMinPtForTOFRequired(0.69),
 106   fMaxPtForTPClow(0.69),
 107
 108   fHEventStat0(NULL),
 109   fHEventStat1(NULL),
 110   fHEventStatMax(6),
 111
 112   fHTriggerStat(NULL),
 113   fNTriggers(5),
 114
 115   fHCentralityStat(NULL),
 116   fNCentralityBins(10),
 117
 118   fNSubSamples(20),
 119   fSubSampleIdx(0),
 120   fRandom(NULL) {
 121   // Constructor
 122
 123   AliLog::SetClassDebugLevel("AliAnalysisNetParticleHelper",10);
 124 }
 125
 126 const Float_t AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap = 0.075;
 127 const Float_t AliAnalysisNetParticleHelper::fgkfHistBinWitdthPt  = 0.3; // 0.08 // 300 MeV  // was 80 MeV
 128
 129 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeCent[]  = {-0.5, 8.5};
 130 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsCent    = 9 ;
 131
 132 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeEta[]   = {-0.9, 0.9};
 133 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsEta     = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangeEta[1] -
 134                                                                           AliAnalysisNetParticleHelper::fgkfHistRangeEta[0]) /
 135                                                                          AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap) +1;
 136
 137 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeRap[]   = {-0.5, 0.5};
 138 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsRap     = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangeRap[1] -
 139                                                                           AliAnalysisNetParticleHelper::fgkfHistRangeRap[0]) /
 140                                                                          AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap) +1;
 141
 142 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangePhi[]   = {0.0, static_cast<Float_t>(TMath::TwoPi())};
 143 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsPhi     = 42 ;
 144
 145 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangePt[]    = {0.2, 2.9}; // {0.2, 5.}; // was {0.3, 2.22}
 146 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsPt      = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangePt[1] -
 147                                                                           AliAnalysisNetParticleHelper::fgkfHistRangePt[0]) /
 148                                                                          AliAnalysisNetParticleHelper::fgkfHistBinWitdthPt);
 149
 150 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeSign[]  = {-1.5, 1.5};
 151 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsSign    =  3;
 152
 153 const Char_t* AliAnalysisNetParticleHelper::fgkEventNames[]         = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,90]%"};
 154 const Char_t* AliAnalysisNetParticleHelper::fgkCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
 155 const Char_t* AliAnalysisNetParticleHelper::fgkTriggerNames[]       = {"kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
 156 const Char_t* AliAnalysisNetParticleHelper::fgkCentralityNames[]    = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%",
 157                                                                        "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
 158
 159 //________________________________________________________________________
 160 AliAnalysisNetParticleHelper::~AliAnalysisNetParticleHelper() {
 161   // Destructor
 162
 163   if (fRandom)
 164     delete fRandom;
 165
 166   return;
 167 }
 168
 169 /*
 170  * ---------------------------------------------------------------------------------
 171  *                                    Setter
 172  * ---------------------------------------------------------------------------------
 173  */
 174
 175 //________________________________________________________________________
 176 void AliAnalysisNetParticleHelper::SetPhiRange(Float_t f1, Float_t f2) {
 177   // -- Set phi range and adopt to phi-histogram
 178
 179   fPhiMin = f1;
 180   fPhiMax = (f1 < f2) ? f2 : f2+TMath::TwoPi();
 181
 182   Float_t phiMin = fPhiMin;
 183   Float_t phiMax = fPhiMax;
 184
 185   // -- Update Ranges
 186   Float_t binWidth = (AliAnalysisNetParticleHelper::fgkfHistRangePhi[1] - AliAnalysisNetParticleHelper::fgkfHistRangePhi[0]) /
 187     Float_t(AliAnalysisNetParticleHelper::fgkfHistNBinsPhi);
 188
 189   Float_t lowEdge  = AliAnalysisNetParticleHelper::fgkfHistRangePhi[0] - binWidth;
 190   Float_t highEdge = AliAnalysisNetParticleHelper::fgkfHistRangePhi[0];
 191
 192   for (Int_t ii = 1; ii <= AliAnalysisNetParticleHelper::fgkfHistNBinsPhi; ++ii) {
 193     lowEdge += binWidth;
 194     highEdge += binWidth;
 195
 196     if (phiMin >= lowEdge && phiMin < highEdge )
 197       phiMin = lowEdge;
 198     if (phiMax > lowEdge && phiMax <= highEdge )
 199       phiMax = highEdge;
 200   }
 201
 202   printf(">>>> Update Phi Range : [%f,%f] -> [%f,%f]\n", fPhiMin, fPhiMax, phiMin, phiMax);
 203   fPhiMin = phiMin;
 204   fPhiMax = phiMax;
 205 }
 206
 207
 208 //________________________________________________________________________
 209 void AliAnalysisNetParticleHelper::SetParticleSpecies(AliPID::EParticleType pid) {
 210   // -- Set particle species (ID, Name, Title, Title LATEX)
 211
 212   if ( Int_t(pid) < 0 || Int_t(pid) >= AliPID::kSPECIES) {
 213     AliWarning("Particle ID not in AliPID::kSPECIES --> Set to protons");
 214     pid = AliPID::kProton;
 215   }
 216
 217   fParticleSpecies = pid;
 218
 219   for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
 220     fPartName[idxPart]       = aPartNames[fParticleSpecies][idxPart];
 221     fPartTitle[idxPart]      = aPartTitles[fParticleSpecies][idxPart];
 222     fPartTitleLatex[idxPart] = aPartTitlesLatex[fParticleSpecies][idxPart];
 223   }
 224 }
 225
 226 //________________________________________________________________________
 227 void AliAnalysisNetParticleHelper::SetUsePID(Bool_t usePID) {
 228   // -- Set usage of PID
 229   //    > if turn off, set charge types (ID, Name, Title, Title LATEX)
 230
 231   fUsePID = usePID;
 232
 233   if (!usePID) {
 234     fParticleSpecies   = AliPID::kUnknown;
 235
 236     fPartName[0]       = "neg";
 237     fPartName[1]       = "pos";
 238     fPartTitle[0]      = "Negative";
 239     fPartTitle[1]      = "Positive";
 240     fPartTitleLatex[0] = "Negative";
 241     fPartTitleLatex[1] = "Positive";
 242   }
 243 }
 244
 245 /*
 246  * ---------------------------------------------------------------------------------
 247  *                                    Getter
 248  * ---------------------------------------------------------------------------------
 249  */
 250
 251 //________________________________________________________________________
 252 TString AliAnalysisNetParticleHelper::GetParticleName(Int_t idxPart) {
 253   // -- Get particle Name
 254
 255   if( idxPart != 0 && idxPart != 1){
 256     AliWarning("Particle type not known --> Set to antiparticles");
 257     idxPart = 0;
 258   }
 259
 260   return fPartName[idxPart];
 261 }
 262
 263 //________________________________________________________________________
 264 TString AliAnalysisNetParticleHelper::GetParticleTitle(Int_t idxPart) {
 265   // -- Get particle Title
 266
 267   if( idxPart != 0 && idxPart != 1){
 268     AliWarning("Particle type not known --> Set to antiparticles");
 269     idxPart = 0;
 270   }
 271
 272   return fPartTitle[idxPart];
 273 }
 274
 275 //________________________________________________________________________
 276 TString AliAnalysisNetParticleHelper::GetParticleTitleLatex(Int_t idxPart) {
 277   // -- Get particle Title LATEX
 278
 279   if( idxPart != 0 && idxPart != 1){
 280     AliWarning("Particle type not known --> Set to antiparticles");
 281     idxPart = 0;
 282   }
 283
 284   return fPartTitleLatex[idxPart];
 285 }
 286
 287 /*
 288  * ---------------------------------------------------------------------------------
 289  *                                 Public Methods
 290  * ---------------------------------------------------------------------------------
 291  */
 292
 293 //________________________________________________________________________
 294 Int_t AliAnalysisNetParticleHelper::Initialize(AliESDtrackCuts *cuts, Bool_t isMC, Int_t trackCutBit, Int_t modeDistCreation) {
 295   // -- Initialize helper
 296
 297   Int_t iResult = 0;
 298
 299   // -- ESD track cuts
 300   fESDTrackCuts     = cuts;
 301
 302   // -- Is MC
 303   fIsMC             = isMC;
 304
 305   // -- AOD track filter bit
 306   fAODtrackCutBit   = trackCutBit;
 307
 308   // -- mode Distribution creation
 309   fModeDistCreation = modeDistCreation;
 310
 311   // -- Setup event cut statistics
 312   InitializeEventStats();
 313
 314   // -- Setup trigger statistics
 315   InitializeTriggerStats();
 316
 317   // -- Setup centrality statistics
 318   InitializeCentralityStats();
 319
 320   // -- PRINT PID Strategy
 321   //    0 :   TPC(TPClow+TPCHigh)
 322   //    1 :   ITS
 323   //    2 :   TOF
 324   //    3 :   ITS+TPC(TPClow+TPCHigh)
 325   //    4 :   TPC(TPClow+TPCHigh)+TOF
 326   //    5 :   TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 327   //    6 :   TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
 328   //    7 :   TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 329   //    8 :   TPC(TPClow+TPCHigh)+ITS+TOF
 330   printf(">>>> USE PID %d || PID STRATEGY: %d || sigmaMax: ITS %.2f TPC %.2f TOF %.2f \n", fUsePID, fPIDStrategy, fNSigmaMaxITS, fNSigmaMaxTPC, fNSigmaMaxTOF);
 331
 332   // -- Initialize random number generator
 333   fRandom = new TRandom3();
 334   fRandom->SetSeed();
 335
 336   return iResult;
 337 }
 338
 339 //________________________________________________________________________
 340 Int_t AliAnalysisNetParticleHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
 341   // -- Setup Event
 342
 343   // -- Get ESD objects
 344   if(esdHandler){
 345     fInputEventHandler = static_cast<AliInputEventHandler*>(esdHandler);
 346     fESD               = dynamic_cast<AliESDEvent*>(fInputEventHandler->GetEvent());
 347     if (!fESD) {
 348       AliError("ESD event handler not available");
 349       return -1;
 350     }
 351   }
 352
 353   // -- Get AOD objects
 354   else if(aodHandler){
 355     fInputEventHandler = static_cast<AliInputEventHandler*>(aodHandler);
 356     fAOD               = dynamic_cast<AliAODEvent*>(fInputEventHandler->GetEvent());
 357     if (!fAOD) {
 358       AliError("AOD event handler not available");
 359       return -1;
 360     }
 361   }
 362
 363   // -- Get Common objects
 364   fPIDResponse = fInputEventHandler->GetPIDResponse();
 365
 366   // -- Get MC objects
 367   fMCEvent     = mcEvent;
 368   if (fMCEvent)
 369     fStack     = fMCEvent->Stack();
 370
 371   // -- Get event centrality
 372   // >  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 373   // >   0   1     2     3     4     5     6     7     8     9
 374
 375   AliCentrality *centrality = NULL;
 376
 377   if(esdHandler)
 378     centrality = fESD->GetCentrality();
 379   else if(aodHandler)
 380     centrality = ((AliVAODHeader*)fAOD->GetHeader())->GetCentralityP();
 381
 382   if (!centrality) {
 383     AliError("Centrality not available");
 384     return -1;
 385   }
 386
 387   Int_t centBin = centrality->GetCentralityClass10("V0M");
 388   if (centBin == 0)
 389     fCentralityBin = centrality->GetCentralityClass5("V0M");
 390   else if (centBin == 10 || centBin == -1.)
 391     fCentralityBin = -1;
 392   else if (centBin > 0 && centBin < fNCentralityBins)
 393     fCentralityBin = centBin + 1;
 394   else
 395     fCentralityBin = -2;
 396
 397   // -- Stay within the max centrality bin
 398   if (fCentralityBin >= fCentralityBinMax)
 399     fCentralityBin = -2;
 400
 401   fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
 402
 403   // -- Get current subsample idx
 404   fSubSampleIdx = fRandom->Integer(fNSubSamples);
 405
 406   return 0;
 407 }
 408
 409 /*
 410  * ---------------------------------------------------------------------------------
 411  *                         Event / Trigger Statistics
 412  * ---------------------------------------------------------------------------------
 413  */
 414
 415 //________________________________________________________________________
 416 Bool_t AliAnalysisNetParticleHelper::IsEventTriggered() {
 417   // -- Check if Event is triggered and fill Trigger Histogram
 418
 419   Bool_t *aTriggerFired = new Bool_t[fNTriggers];
 420   for (Int_t ii = 0; ii < fNTriggers; ++ii)
 421     aTriggerFired[ii] = kFALSE;
 422
 423   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kMB))          aTriggerFired[0] = kTRUE;
 424   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kCentral))     aTriggerFired[1] = kTRUE;
 425   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
 426   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      aTriggerFired[3] = kTRUE;
 427   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      aTriggerFired[4] = kTRUE;
 428
 429   Bool_t isTriggered = kFALSE;
 430
 431   for (Int_t ii=0; ii<fNTriggers; ++ii) {
 432     if(aTriggerFired[ii]) {
 433       isTriggered = kTRUE;
 434       fHTriggerStat->Fill(ii);
 435     }
 436   }
 437
 438   delete[] aTriggerFired;
 439
 440   return isTriggered;
 441 }
 442
 443 //________________________________________________________________________
 444 Bool_t AliAnalysisNetParticleHelper::IsEventRejected() {
 445   // -- Evaluate event statistics histograms
 446
 447   Int_t *aEventCuts = new Int_t[fHEventStatMax];
 448   // set aEventCuts[ii] to 1 in case of reject
 449
 450   for (Int_t ii=0;ii<fHEventStatMax; ++ii)
 451     aEventCuts[ii] = 0;
 452
 453   Int_t iCut = 0;
 454
 455   // -- 0 - Before Physics Selection
 456   aEventCuts[iCut] = 0;
 457
 458   // -- 1 - No Trigger fired
 459   ++iCut;
 460   if (!IsEventTriggered())
 461     aEventCuts[iCut] = 1;
 462
 463   // -- 2 - No Vertex
 464   ++iCut;
 465   const AliESDVertex* vtxESD = NULL;
 466   const AliAODVertex* vtxAOD = NULL;
 467   if (fESD){
 468     vtxESD = fESD->GetPrimaryVertexTracks();
 469     if (!vtxESD)
 470       aEventCuts[iCut] = 1;
 471   }
 472   else if (fAOD){
 473     vtxAOD = fAOD->GetPrimaryVertex();
 474     if (!vtxAOD)
 475       aEventCuts[iCut] = 1;
 476   }
 477
 478   // -- 3 - Vertex z outside cut window
 479   ++iCut;
 480   if (vtxESD){
 481     if(TMath::Abs(vtxESD->GetZ()) > fVertexZMax)
 482       aEventCuts[iCut] = 1;
 483   }
 484   else if(vtxAOD){
 485     if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax)
 486       aEventCuts[iCut] = 1;
 487   }
 488   else
 489     aEventCuts[iCut] = 1;
 490
 491   // -- 4 - Centrality = -1  (no centrality or not hadronic)
 492   ++iCut;
 493   if(fCentralityBin == -1.)
 494     aEventCuts[iCut] = 1;
 495
 496   // -- 5 - Centrality < fCentralityMax
 497   ++iCut;
 498   if(fCentralityBin == -2.)
 499     aEventCuts[iCut] = 1;
 500
 501   // -- Fill statistics / reject event
 502   Bool_t isRejected = FillEventStats(aEventCuts);
 503
 504   // -- Cleanup
 505   delete[] aEventCuts;
 506
 507   return isRejected;
 508 }
 509
 510 /*
 511  * ---------------------------------------------------------------------------------
 512  *                         Accept Particle Methods - private
 513  * ---------------------------------------------------------------------------------
 514  */
 515
 516 //________________________________________________________________________
 517 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC) {
 518   // -- Check if MC particle is accepted for basic parameters
 519
 520   if (!particle)
 521     return kFALSE;
 522
 523   // -- check if charged
 524   if (particle->Charge() == 0.0)
 525     return kFALSE;
 526
 527   // -- check if physical primary - ESD
 528   if (fESD) {
 529     if(!fStack->IsPhysicalPrimary(idxMC))
 530       return kFALSE;
 531   }
 532   // -- check if physical primary - AOD
 533   else {
 534     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
 535       return kFALSE;
 536   }
 537
 538   return kTRUE;
 539 }
 540
 541 //________________________________________________________________________
 542 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC) {
 543   // -- Check if MC particle is accepted for basic parameters
 544
 545   if (!particle)
 546     return kFALSE;
 547
 548   // -- check if charged
 549   if (particle->Charge() != 0.0)
 550     return kFALSE;
 551
 552   // -- check if physical primary - ESD
 553   if (fESD) {
 554     if(!fStack->IsPhysicalPrimary(idxMC))
 555       return kFALSE;
 556   }
 557   // -- check if physical primary - AOD
 558   else {
 559     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
 560       return kFALSE;
 561   }
 562
 563   return kTRUE;
 564 }
 565
 566 //________________________________________________________________________
 567 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP) {
 568   // -- Check if particle is accepted
 569   // > in rapidity
 570   // > if no pid : return kTRUE, yP = eta
 571   // > return 0 if not accepted
 572
 573   if (!fUsePID) {
 574     yP = particle->Eta();
 575     return kTRUE;
 576   }
 577
 578   Double_t mP = AliPID::ParticleMass(fParticleSpecies);
 579
 580   // -- Calculate rapidities and kinematics
 581   Double_t p  = particle->P();
 582   Double_t pz = particle->Pz();
 583
 584   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 585   yP          = 0.5 * TMath::Log((eP + pz) / (eP - pz));
 586
 587   // -- Check Rapidity window
 588   if (TMath::Abs(yP) > fRapidityMax)
 589     return kFALSE;
 590
 591   return kTRUE;
 592 }
 593
 594 //________________________________________________________________________
 595 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedPhi(AliVParticle *particle) {
 596   // -- Check if particle is accepted
 597   // > in phi
 598   // > return 0 if not accepted
 599
 600   if (particle->Phi() > fPhiMin && particle->Phi() <= fPhiMax)
 601     return kTRUE;
 602   else if (particle->Phi() < fPhiMin && (particle->Phi() + TMath::TwoPi()) <= fPhiMax)
 603     return kTRUE;
 604   else
 605     return kFALSE;
 606 }
 607
 608 //_____________________________________________________________________________
 609 Bool_t AliAnalysisNetParticleHelper::IsParticleFindable(Int_t label) {
 610   // -- Check if MC particle is findable tracks
 611
 612   AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
 613   if(!mcParticle)
 614     return kFALSE;
 615
 616   Int_t counter;
 617   Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
 618
 619   return (tpcTrackLength > fMinTrackLengthMC);
 620 }
 621
 622 /*
 623  * ---------------------------------------------------------------------------------
 624  *                            Accept Track Methods - public
 625  * ---------------------------------------------------------------------------------
 626  */
 627
 628 //________________________________________________________________________
 629 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliVTrack* track) {
 630   // -- Check if track is accepted
 631   // > for basic parameters
 632
 633   if (!track)
 634     return kFALSE;
 635
 636   if (track->Charge() == 0)
 637     return kFALSE;
 638
 639   return kTRUE;
 640 }
 641
 642 //________________________________________________________________________
 643 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP) {
 644   // -- Check if track is accepted
 645   // > in rapidity
 646   // > if no pid : return kTRUE
 647   // > return 0 if not accepted
 648
 649   if (!fUsePID) {
 650     yP = track->Eta();
 651     return kTRUE;
 652   }
 653
 654   Double_t mP = AliPID::ParticleMass(fParticleSpecies);
 655
 656   // -- Calculate rapidities and kinematics
 657   Double_t pvec[3];
 658   track->GetPxPyPz(pvec);
 659
 660   Double_t p  = track->P();
 661   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 662            yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
 663
 664   // -- Check Rapidity window
 665   if (TMath::Abs(yP) > fRapidityMax)
 666     return kFALSE;
 667
 668   return kTRUE;
 669 }
 670
 671 //________________________________________________________________________
 672 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedDCA(AliVTrack *vTrack) {
 673   // -- Check if track is accepted - ONLY FOR ESDs so far
 674   // > for DCA, if both SPD layers have hits
 675   // > For now only Implemented for ESDs
 676
 677   Bool_t isAccepted = kTRUE;
 678
 679   if (!fESD)
 680     return isAccepted;
 681
 682   AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
 683   if (!track)
 684     return kFALSE;
 685
 686   // -- Get nHits SPD
 687   if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
 688
 689     // -- Get DCA nSigmas
 690     Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 691     track->GetImpactParameters(dca,cov);
 692
 693     Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
 694     Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
 695
 696     if (fNSigmaMaxCdd != 0.) {
 697       if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
 698         isAccepted = kFALSE;
 699     }
 700
 701     if (fNSigmaMaxCzz != 0.) {
 702       if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
 703         isAccepted = kFALSE;
 704     }
 705   }
 706
 707   return isAccepted;
 708 }
 709
 710 //________________________________________________________________________
 711 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid) {
 712   // -- Check if track is accepted
 713   // > provides ITS, TPC and TOF nSigmas to the argument
 714
 715   Bool_t isAcceptedITS    = kFALSE;
 716   Bool_t isAcceptedTPC    = kFALSE;
 717   Bool_t isAcceptedTPClow = kFALSE;
 718   Bool_t isAcceptedTOF    = kFALSE;
 719   Bool_t isAccepted       = kFALSE;
 720
 721   // -- In case not PID is used
 722   if (!fUsePID) {
 723     pid[0] = 10.;
 724     pid[1] = 10.;
 725     pid[2] = 10.;
 726     return kTRUE;
 727   }
 728
 729   // -- Get PID with ITS and check
 730   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kITS, track, fParticleSpecies, pid[0]) == AliPIDResponse::kDetPidOk) {
 731     if (TMath::Abs(pid[0]) < fNSigmaMaxITS)
 732       isAcceptedITS = kTRUE;
 733   }
 734
 735   // -- Get PID with TPC and check
 736   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, fParticleSpecies, pid[1]) == AliPIDResponse::kDetPidOk) {
 737     if (TMath::Abs(pid[1]) < fNSigmaMaxTPC)
 738       isAcceptedTPC = kTRUE;
 739     if (TMath::Abs(pid[1]) < fNSigmaMaxTPClow)
 740       isAcceptedTPClow = kTRUE;
 741     if (track->Pt() < fMaxPtForTPClow)
 742       isAcceptedTPC = isAcceptedTPClow;
 743   }
 744
 745   // -- Get PID with TOF and check
 746   Bool_t hasPIDTOF = kFALSE;
 747   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, fParticleSpecies, pid[2]) == AliPIDResponse::kDetPidOk) {
 748     hasPIDTOF = kTRUE;
 749     if (TMath::Abs(pid[2]) < fNSigmaMaxTOF)
 750       isAcceptedTOF = kTRUE;
 751   }
 752
 753   // -- Check TOF missmatch for MC
 754
 755   //if (ESD)
 756   if (fIsMC && isAcceptedTOF) {
 757     Int_t tofLabel[3];
 758     //    AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
 759     // TODO add code for AOD
 760
 761     (dynamic_cast<AliESDtrack*>(track))->GetTOFLabel(tofLabel);
 762
 763     Bool_t hasMatchTOF = kTRUE;
 764     if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0)
 765       hasMatchTOF = kFALSE;
 766
 767     TParticle *matchedTrack = fStack->Particle(TMath::Abs(tofLabel[0]));
 768     if (TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel()))
 769       hasMatchTOF = kTRUE;
 770
 771     isAcceptedTOF = hasMatchTOF;
 772   }
 773
 774   //    0 :   TPC(TPClow+TPCHigh)
 775   //    1 :   ITS
 776   //    2 :   TOF
 777   //    3 :   ITS+TPC(TPClow+TPCHigh)
 778   //    4 :   TPC(TPClow+TPCHigh)+TOF
 779   //    5 :   TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 780   //    6 :   TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
 781   //    7 :   TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
 782   //    8 :   TPC(TPClow+TPCHigh)+ITS+TOF
 783   if (fPIDStrategy == 0) {             //  TPC PID
 784     isAccepted = isAcceptedTPC;
 785   }
 786   else if (fPIDStrategy == 1) {        //  ITS PID
 787     isAccepted = isAcceptedITS;
 788   }
 789   else if (fPIDStrategy == 2) {        //  TOF PID
 790     isAccepted = isAcceptedTOF;
 791   }
 792   else if (fPIDStrategy == 3) {        //  TPC+ITS PID
 793     isAccepted = isAcceptedTPC && isAcceptedITS;
 794   }
 795   else if (fPIDStrategy == 4) {        //  TPC+TOF PID
 796     isAccepted = isAcceptedTPC && isAcceptedTOF;
 797   }
 798   else if (fPIDStrategy == 5) {        //  TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
 799     if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
 800       isAcceptedTOF = kTRUE;
 801     isAccepted = isAcceptedTPC && isAcceptedTOF;
 802   }
 803   else if (fPIDStrategy == 6) {        //  ITS+TPC+TOF PID -- TOF only for those tracks which have TOF information
 804     isAccepted = isAcceptedTPC && isAcceptedITS;
 805     if (hasPIDTOF)
 806       isAccepted = isAccepted && isAcceptedTOF;
 807   }
 808   else if (fPIDStrategy == 7) {        //  ITS+TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
 809     if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
 810       isAcceptedTOF = kTRUE;
 811     isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
 812   }
 813   else if (fPIDStrategy == 8) {        //  ITS+TPC+TOF PID
 814     isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
 815   }
 816
 817   return isAccepted;
 818 }
 819
 820 //________________________________________________________________________
 821 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPhi(AliVTrack *track) {
 822   // -- Check if track is accepted
 823   // > in phi
 824   // > return 0 if not accepted
 825
 826   if (track->Phi() > fPhiMin && track->Phi() <= fPhiMax)
 827     return kTRUE;
 828   else if (track->Phi() < fPhiMin && (track->Phi() + TMath::TwoPi()) <= fPhiMax)
 829     return kTRUE;
 830   else
 831     return kFALSE;
 832 }
 833
 834 /*
 835  * ---------------------------------------------------------------------------------
 836  *                         Helper Methods
 837  * ---------------------------------------------------------------------------------
 838  */
 839
 840 //________________________________________________________________________
 841 void AliAnalysisNetParticleHelper::BinLogAxis(const THnBase *hn, Int_t axisNumber, AliESDtrackCuts* cuts) {
 842   // -- Method for the correct logarithmic binning of histograms
 843   // -- and update fMinPtForTOFRequired using the logarithmic scale
 844
 845   AliESDtrackCuts* esdTrackCuts = (cuts) ? cuts : fESDTrackCuts;
 846
 847   // -- Make logarithmic binning
 848   TAxis *axis = hn->GetAxis(axisNumber);
 849   Int_t  nBins = axis->GetNbins();
 850
 851   Double_t from  = axis->GetXmin();
 852   Double_t to    = axis->GetXmax();
 853   Double_t *newBins = new Double_t[nBins + 1];
 854
 855   newBins[0] = from;
 856   Double_t factor = TMath::Power(to/from, 1./nBins);
 857
 858   for (int ii = 1; ii <= nBins; ii++)
 859    newBins[ii] = factor * newBins[ii-1];
 860
 861   axis->Set(nBins, newBins);
 862
 863   delete [] newBins;
 864
 865   // -- Update Ranges
 866   // ------------------
 867   Float_t ptRange[2];
 868   Float_t oldPtRange[2];
 869   esdTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 870   esdTrackCuts->GetPtRange(oldPtRange[0],oldPtRange[1]);
 871
 872   Float_t minPtForTOFRequired = fMinPtForTOFRequired;
 873   Float_t maxPtForTPClow      = fMaxPtForTPClow;
 874
 875   // -- Update minPt Cut
 876   Int_t bin = axis->FindBin(ptRange[0]+10e-6);
 877   ptRange[0] = axis->GetBinLowEdge(bin);
 878
 879   // -- Update maxPt Cut
 880   bin = axis->FindBin(ptRange[1]-10e-6);
 881   ptRange[1] = axis->GetBinUpEdge(bin);
 882
 883   if (ptRange[0] != oldPtRange[0] || ptRange[1] != oldPtRange[1]) {
 884     printf(">>>> Update Pt Range : [%f,%f] -> [%f,%f]\n", oldPtRange[0], oldPtRange[1], ptRange[0], ptRange[1]);
 885     esdTrackCuts->SetPtRange(ptRange[0],ptRange[1]);
 886   }
 887
 888   // -- Update MinPtForTOFRequired
 889   bin = axis->FindBin(minPtForTOFRequired-10e-6);
 890   minPtForTOFRequired = axis->GetBinUpEdge(bin);
 891
 892   if (minPtForTOFRequired != fMinPtForTOFRequired) {
 893     printf(">>>> Update Min Pt for TOF : %f -> %f\n", fMinPtForTOFRequired, minPtForTOFRequired);
 894     fMinPtForTOFRequired = minPtForTOFRequired;
 895   }
 896
 897   // -- Update MaxPtForTPClow
 898   bin = axis->FindBin(maxPtForTPClow-10e-6);
 899   maxPtForTPClow = axis->GetBinUpEdge(bin);
 900
 901   if (maxPtForTPClow != fMaxPtForTPClow) {
 902     printf(">>>> Update Max Pt for TPC Low : %f -> %f\n", fMaxPtForTPClow, maxPtForTPClow);
 903     fMaxPtForTPClow = maxPtForTPClow;
 904   }
 905 }
 906
 907 ///////////////////////////////////////////////////////////////////////////////////
 908
 909 /*
 910  * ---------------------------------------------------------------------------------
 911  *                           Initialize - Private
 912  * ---------------------------------------------------------------------------------
 913  */
 914
 915 //________________________________________________________________________
 916 void AliAnalysisNetParticleHelper::InitializeEventStats() {
 917   // -- Initialize event statistics histograms
 918
 919   fHEventStat0 = new TH1F("hEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 920   fHEventStat1 = new TH1F("hEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 921
 922   for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
 923     fHEventStat0->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkEventNames[ii]);
 924     fHEventStat1->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkEventNames[ii]);
 925   }
 926
 927   fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliAnalysisNetParticleHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
 928   fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliAnalysisNetParticleHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
 929 }
 930
 931 //________________________________________________________________________
 932 void AliAnalysisNetParticleHelper::InitializeTriggerStats() {
 933   // -- Initialize trigger statistics histograms
 934
 935   fHTriggerStat = new TH1F("hTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
 936
 937   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 938     fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkTriggerNames[ii]);
 939 }
 940
 941 //________________________________________________________________________
 942 void AliAnalysisNetParticleHelper::InitializeCentralityStats() {
 943   // -- Initialize trigger statistics histograms
 944
 945   fHCentralityStat = new TH1F("hCentralityStat","Centrality statistics;Centrality Bins;Events",
 946                               fNCentralityBins,-0.5,fNCentralityBins-0.5);
 947
 948   for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
 949     fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkCentralityNames[ii]);
 950 }
 951
 952 /*
 953  * ---------------------------------------------------------------------------------
 954  *                         Event / Trigger Statistics - private
 955  * ---------------------------------------------------------------------------------
 956  */
 957
 958 //________________________________________________________________________
 959 Bool_t AliAnalysisNetParticleHelper::FillEventStats(Int_t *aEventCuts) {
 960   // -- Fill event / centrality statistics
 961
 962   Bool_t isRejected = kFALSE;
 963
 964   // -- Fill event statistics
 965   for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
 966
 967     if (aEventCuts[idx])
 968       isRejected = kTRUE;
 969     else
 970       fHEventStat0->Fill(idx);
 971   }
 972
 973   for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
 974     if (aEventCuts[idx])
 975       break;
 976     fHEventStat1->Fill(idx);
 977   }
 978
 979   // -- Fill centrality statistics of accepted events
 980   if (!isRejected)
 981     fHCentralityStat->Fill(fCentralityBin);
 982
 983   return isRejected;
 984 }