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 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   fAOD(NULL),
  79   fMCEvent(NULL),
  80   fStack(NULL),
  81
  82   fCentralityBin(-1),
  83   fCentralityPercentile(-1.),
  84
  85   fCentralityBinMax(7),
  86   fVertexZMax(10.),
  87   fRapidityMax(0.5),
  88   fPhiMin(0.),
  89   fPhiMax(TMath::TwoPi()),
  90   fMinTrackLengthMC(70.),
  91   fNSigmaMaxCdd(3.),
  92   fNSigmaMaxCzz(3.),
  93
  94   fParticleSpecies(AliPID::kProton),
  95
  96   fUsePID(kTRUE),
  97   fNSigmaMaxTPC(2.5),
  98   fNSigmaMaxTOF(2.5),
  99   fMinPtForTOFRequired(0.8),
 100
 101   fHEventStat0(NULL),
 102   fHEventStat1(NULL),
 103   fHEventStatMax(6),
 104
 105   fHTriggerStat(NULL),
 106   fNTriggers(5),
 107
 108   fHCentralityStat(NULL),
 109   fNCentralityBins(10),
 110
 111   fEtaCorrFunc(NULL) {
 112   // Constructor
 113
 114   AliLog::SetClassDebugLevel("AliAnalysisNetParticleHelper",10);
 115 }
 116
 117 const Float_t AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap = 0.075;
 118 const Float_t AliAnalysisNetParticleHelper::fgkfHistBinWitdthPt  = 0.15;   // 150 MeV
 119
 120 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeCent[]  = {-0.5, 8.5};
 121 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsCent    = 9 ;
 122
 123 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeEta[]   = {-0.9, 0.9};
 124 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsEta     = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangeEta[1] -
 125                                                                           AliAnalysisNetParticleHelper::fgkfHistRangeEta[0]) /
 126                                                                          AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap) +1;
 127
 128 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeRap[]   = {-0.5, 0.5};
 129 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsRap     = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangeRap[1] -
 130                                                                           AliAnalysisNetParticleHelper::fgkfHistRangeRap[0]) /
 131                                                                          AliAnalysisNetParticleHelper::fgkfHistBinWitdthRap) +1;
 132
 133 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangePhi[]   = {0.0, static_cast<Float_t>(TMath::TwoPi())};
 134 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsPhi     = 42 ;
 135
 136 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangePt[]    = {0.1, 3.0};
 137 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsPt      = Int_t((AliAnalysisNetParticleHelper::fgkfHistRangePt[1] -
 138                                                                           AliAnalysisNetParticleHelper::fgkfHistRangePt[0]) /
 139                                                                          AliAnalysisNetParticleHelper::fgkfHistBinWitdthPt);
 140
 141 const Float_t AliAnalysisNetParticleHelper::fgkfHistRangeSign[]  = {-1.5, 1.5};
 142 const Int_t   AliAnalysisNetParticleHelper::fgkfHistNBinsSign    =  3;
 143
 144 const Char_t* AliAnalysisNetParticleHelper::fgkEventNames[]         = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,90]%"};
 145 const Char_t* AliAnalysisNetParticleHelper::fgkCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
 146 const Char_t* AliAnalysisNetParticleHelper::fgkTriggerNames[]       = {"kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
 147 const Char_t* AliAnalysisNetParticleHelper::fgkCentralityNames[]    = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%",
 148                                                                        "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
 149
 150 //________________________________________________________________________
 151 AliAnalysisNetParticleHelper::~AliAnalysisNetParticleHelper() {
 152   // Destructor
 153
 154   return;
 155 }
 156
 157 /*
 158  * ---------------------------------------------------------------------------------
 159  *                                    Setter
 160  * ---------------------------------------------------------------------------------
 161  */
 162
 163 //________________________________________________________________________
 164 void AliAnalysisNetParticleHelper::SetPhiRange(Float_t f1, Float_t f2) {
 165   // -- Set phi range and adopt to phi-histogram
 166
 167   fPhiMin = f1;
 168   fPhiMax = (f1 < f2) ? f2 : f2+TMath::TwoPi();
 169
 170   // -- Update Ranges
 171   Float_t binWidth = (AliAnalysisNetParticleHelper::fgkfHistRangePhi[1] - AliAnalysisNetParticleHelper::fgkfHistRangePhi[0]) /
 172     Float_t(AliAnalysisNetParticleHelper::fgkfHistNBinsPhi);
 173
 174   Float_t lowEdge  = AliAnalysisNetParticleHelper::fgkfHistRangePhi[0] - binWidth;
 175   Float_t highEdge = AliAnalysisNetParticleHelper::fgkfHistRangePhi[0];
 176
 177   for (Int_t ii = 1; ii <= AliAnalysisNetParticleHelper::fgkfHistNBinsPhi; ++ii) {
 178     lowEdge += binWidth;
 179     highEdge += binWidth;
 180
 181     if (fPhiMin >= lowEdge && fPhiMin < highEdge )
 182       fPhiMin = lowEdge;
 183     if (fPhiMax > lowEdge && fPhiMax <= highEdge )
 184       fPhiMax = highEdge;
 185   }
 186
 187 }
 188
 189
 190 //________________________________________________________________________
 191 void AliAnalysisNetParticleHelper::SetParticleSpecies(AliPID::EParticleType pid) {
 192   // -- Set particle species (ID, Name, Title, Title LATEX)
 193
 194   if ( Int_t(pid) < 0 || Int_t(pid) >= AliPID::kSPECIES) {
 195     AliWarning("Particle ID not in AliPID::kSPECIES --> Set to protons");
 196     pid = AliPID::kProton;
 197   }
 198
 199   fParticleSpecies = pid;
 200
 201   for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
 202     fPartName[idxPart]       = aPartNames[fParticleSpecies][idxPart];
 203     fPartTitle[idxPart]      = aPartTitles[fParticleSpecies][idxPart];
 204     fPartTitleLatex[idxPart] = aPartTitlesLatex[fParticleSpecies][idxPart];
 205   }
 206 }
 207
 208 //________________________________________________________________________
 209 void AliAnalysisNetParticleHelper::SetUsePID(Bool_t usePID) {
 210   // -- Set usage of PID
 211   //    > if turn off, set charge types (ID, Name, Title, Title LATEX)
 212
 213   fUsePID = usePID;
 214
 215   if (!usePID) {
 216     fParticleSpecies   = AliPID::kUnknown;
 217
 218     fPartName[0]       = "neg";
 219     fPartName[1]       = "pos";
 220     fPartTitle[0]      = "Negative";
 221     fPartTitle[1]      = "Positive";
 222     fPartTitleLatex[0] = "Negative";
 223     fPartTitleLatex[1] = "Positive";
 224   }
 225 }
 226
 227 /*
 228  * ---------------------------------------------------------------------------------
 229  *                                    Getter
 230  * ---------------------------------------------------------------------------------
 231  */
 232
 233 //________________________________________________________________________
 234 TString AliAnalysisNetParticleHelper::GetParticleName(Int_t idxPart) {
 235   // -- Get particle Name
 236
 237   if( idxPart != 0 && idxPart != 1){
 238     AliWarning("Particle type not known --> Set to antiparticles");
 239     idxPart = 0;
 240   }
 241
 242   return fPartName[idxPart];
 243 }
 244
 245 //________________________________________________________________________
 246 TString AliAnalysisNetParticleHelper::GetParticleTitle(Int_t idxPart) {
 247   // -- Get particle Title
 248
 249   if( idxPart != 0 && idxPart != 1){
 250     AliWarning("Particle type not known --> Set to antiparticles");
 251     idxPart = 0;
 252   }
 253
 254   return fPartTitle[idxPart];
 255 }
 256
 257 //________________________________________________________________________
 258 TString AliAnalysisNetParticleHelper::GetParticleTitleLatex(Int_t idxPart) {
 259   // -- Get particle Title LATEX
 260
 261   if( idxPart != 0 && idxPart != 1){
 262     AliWarning("Particle type not known --> Set to antiparticles");
 263     idxPart = 0;
 264   }
 265
 266   return fPartTitleLatex[idxPart];
 267 }
 268
 269 /*
 270  * ---------------------------------------------------------------------------------
 271  *                                 Public Methods
 272  * ---------------------------------------------------------------------------------
 273  */
 274
 275 //________________________________________________________________________
 276 Int_t AliAnalysisNetParticleHelper::Initialize(Bool_t isMC, Int_t modeDistCreation) {
 277   // -- Initialize helper
 278
 279   Int_t iResult = 0;
 280   fModeDistCreation = modeDistCreation;
 281
 282   // -- Setup event cut statistics
 283   InitializeEventStats();
 284
 285   // -- Setup trigger statistics
 286   InitializeTriggerStats();
 287
 288   // -- Setup centrality statistics
 289   InitializeCentralityStats();
 290
 291   // -- Load Eta correction function
 292   iResult = InitializeEtaCorrection(isMC);
 293
 294   return iResult;
 295 }
 296
 297 //________________________________________________________________________
 298 Int_t AliAnalysisNetParticleHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
 299   // -- Setup Event
 300
 301   // -- Get ESD objects
 302   if(esdHandler){
 303     fInputEventHandler = static_cast<AliInputEventHandler*>(esdHandler);
 304     fESD               = dynamic_cast<AliESDEvent*>(fInputEventHandler->GetEvent());
 305     if (!fESD) {
 306       AliError("ESD event handler not available");
 307       return -1;
 308     }
 309   }
 310
 311   // -- Get AOD objects
 312   else if(aodHandler){
 313     fInputEventHandler = static_cast<AliInputEventHandler*>(aodHandler);
 314     fAOD               = dynamic_cast<AliAODEvent*>(fInputEventHandler->GetEvent());
 315     if (!fAOD) {
 316       AliError("AOD event handler not available");
 317       return -1;
 318     }
 319   }
 320
 321   // -- Get Common objects
 322   fPIDResponse = fInputEventHandler->GetPIDResponse();
 323
 324   // -- Get MC objects
 325   fMCEvent     = mcEvent;
 326   if (fMCEvent)
 327     fStack     = fMCEvent->Stack();
 328
 329   // -- Get event centrality
 330   // >  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 331   // >   0   1     2     3     4     5     6     7     8     9
 332
 333   AliCentrality *centrality = NULL;
 334
 335   if(esdHandler)
 336     centrality = fESD->GetCentrality();
 337   else if(aodHandler)
 338     centrality = fAOD->GetHeader()->GetCentralityP();
 339
 340   if (!centrality) {
 341     AliError("Centrality not available");
 342     return -1;
 343   }
 344
 345   Int_t centBin = centrality->GetCentralityClass10("V0M");
 346   if (centBin == 0)
 347     fCentralityBin = centrality->GetCentralityClass5("V0M");
 348   else if (centBin == 10 || centBin == -1.)
 349     fCentralityBin = -1;
 350   else if (centBin > 0 && centBin < fNCentralityBins)
 351     fCentralityBin = centBin + 1;
 352   else
 353     fCentralityBin = -2;
 354
 355   // -- Stay within the max centrality bin
 356   if (fCentralityBin >= fCentralityBinMax)
 357     fCentralityBin = -2;
 358
 359   fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
 360
 361   // -- Update TPC pid with eta correction (only for ESDs!!)
 362   if(esdHandler)
 363     UpdateEtaCorrectedTPCPid();
 364
 365   return 0;
 366 }
 367
 368 /*
 369  * ---------------------------------------------------------------------------------
 370  *                         Event / Trigger Statistics
 371  * ---------------------------------------------------------------------------------
 372  */
 373
 374 //________________________________________________________________________
 375 Bool_t AliAnalysisNetParticleHelper::IsEventTriggered() {
 376   // -- Check if Event is triggered and fill Trigger Histogram
 377
 378   Bool_t *aTriggerFired = new Bool_t[fNTriggers];
 379   for (Int_t ii = 0; ii < fNTriggers; ++ii)
 380     aTriggerFired[ii] = kFALSE;
 381
 382   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kMB))          aTriggerFired[0] = kTRUE;
 383   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kCentral))     aTriggerFired[1] = kTRUE;
 384   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
 385   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      aTriggerFired[3] = kTRUE;
 386   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      aTriggerFired[4] = kTRUE;
 387
 388   Bool_t isTriggered = kFALSE;
 389
 390   for (Int_t ii=0; ii<fNTriggers; ++ii) {
 391     if(aTriggerFired[ii]) {
 392       isTriggered = kTRUE;
 393       fHTriggerStat->Fill(ii);
 394     }
 395   }
 396
 397   delete[] aTriggerFired;
 398
 399   return isTriggered;
 400 }
 401
 402 //________________________________________________________________________
 403 Bool_t AliAnalysisNetParticleHelper::IsEventRejected() {
 404   // -- Evaluate event statistics histograms
 405
 406   Int_t *aEventCuts = new Int_t[fHEventStatMax];
 407   // set aEventCuts[ii] to 1 in case of reject
 408
 409   for (Int_t ii=0;ii<fHEventStatMax; ++ii)
 410     aEventCuts[ii] = 0;
 411
 412   Int_t iCut = 0;
 413
 414   // -- 0 - Before Physics Selection
 415   aEventCuts[iCut] = 0;
 416
 417   // -- 1 - No Trigger fired
 418   ++iCut;
 419   if (!IsEventTriggered())
 420     aEventCuts[iCut] = 1;
 421
 422   // -- 2 - No Vertex
 423   ++iCut;
 424   const AliESDVertex* vtxESD = NULL;
 425   const AliAODVertex* vtxAOD = NULL;
 426   if (fESD){
 427     vtxESD = fESD->GetPrimaryVertexTracks();
 428     if (!vtxESD)
 429       aEventCuts[iCut] = 1;
 430   }
 431   else if (fAOD){
 432     vtxAOD = fAOD->GetPrimaryVertex();
 433     if (!vtxAOD)
 434       aEventCuts[iCut] = 1;
 435   }
 436
 437   // -- 3 - Vertex z outside cut window
 438   ++iCut;
 439   if (vtxESD){
 440     if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax)
 441       aEventCuts[iCut] = 1;
 442   }
 443   else if(vtxAOD){
 444     if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax)
 445       aEventCuts[iCut] = 1;
 446   }
 447   else
 448     aEventCuts[iCut] = 1;
 449
 450   // -- 4 - Centrality = -1  (no centrality or not hadronic)
 451   ++iCut;
 452   if(fCentralityBin == -1.)
 453     aEventCuts[iCut] = 1;
 454
 455   // -- 5 - Centrality < fCentralityMax
 456   ++iCut;
 457   if(fCentralityBin == -2.)
 458     aEventCuts[iCut] = 1;
 459
 460   // -- Fill statistics / reject event
 461   Bool_t isRejected = FillEventStats(aEventCuts);
 462
 463   // -- Cleanup
 464   delete[] aEventCuts;
 465
 466   return isRejected;
 467 }
 468
 469 /*
 470  * ---------------------------------------------------------------------------------
 471  *                         Accept Particle Methods - private
 472  * ---------------------------------------------------------------------------------
 473  */
 474
 475 //________________________________________________________________________
 476 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC) {
 477   // -- Check if MC particle is accepted for basic parameters
 478
 479   if (!particle)
 480     return kFALSE;
 481
 482   // -- check if charged
 483   if (particle->Charge() == 0.0)
 484     return kFALSE;
 485
 486   // -- check if physical primary - ESD
 487   if (fESD) {
 488     if(!fStack->IsPhysicalPrimary(idxMC))
 489       return kFALSE;
 490   }
 491   // -- check if physical primary - AOD
 492   else {
 493     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
 494       return kFALSE;
 495   }
 496
 497   return kTRUE;
 498 }
 499
 500 //________________________________________________________________________
 501 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC) {
 502   // -- Check if MC particle is accepted for basic parameters
 503
 504   if (!particle)
 505     return kFALSE;
 506
 507   // -- check if charged
 508   if (particle->Charge() != 0.0)
 509     return kFALSE;
 510
 511   // -- check if physical primary - ESD
 512   if (fESD) {
 513     if(!fStack->IsPhysicalPrimary(idxMC))
 514       return kFALSE;
 515   }
 516   // -- check if physical primary - AOD
 517   else {
 518     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary())
 519       return kFALSE;
 520   }
 521
 522   return kTRUE;
 523 }
 524
 525 //________________________________________________________________________
 526 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP) {
 527   // -- Check if particle is accepted
 528   // > in rapidity
 529   // > if no pid : return kTRUE, yP = eta
 530   // > return 0 if not accepted
 531
 532   if (!fUsePID) {
 533     yP = particle->Eta();
 534     return kTRUE;
 535   }
 536
 537   Double_t mP = AliPID::ParticleMass(fParticleSpecies);
 538
 539   // -- Calculate rapidities and kinematics
 540   Double_t p  = particle->P();
 541   Double_t pz = particle->Pz();
 542
 543   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 544   yP          = 0.5 * TMath::Log((eP + pz) / (eP - pz));
 545
 546   // -- Check Rapidity window
 547   if (TMath::Abs(yP) > fRapidityMax)
 548     return kFALSE;
 549
 550   return kTRUE;
 551 }
 552
 553 //________________________________________________________________________
 554 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedPhi(AliVParticle *particle) {
 555   // -- Check if particle is accepted
 556   // > in phi
 557   // > return 0 if not accepted
 558
 559   if (particle->Phi() > fPhiMin && particle->Phi() <= fPhiMax)
 560     return kTRUE;
 561   else if (particle->Phi() < fPhiMin && (particle->Phi() + TMath::TwoPi()) <= fPhiMax)
 562     return kTRUE;
 563   else
 564     return kFALSE;
 565 }
 566
 567 //_____________________________________________________________________________
 568 Bool_t AliAnalysisNetParticleHelper::IsParticleFindable(Int_t label) {
 569   // -- Check if MC particle is findable tracks
 570
 571   AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
 572   if(!mcParticle)
 573     return kFALSE;
 574
 575   Int_t counter;
 576   Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
 577
 578   return (tpcTrackLength > fMinTrackLengthMC);
 579 }
 580
 581 /*
 582  * ---------------------------------------------------------------------------------
 583  *                            Accept Track Methods - public
 584  * ---------------------------------------------------------------------------------
 585  */
 586
 587 //________________________________________________________________________
 588 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliVTrack* track) {
 589   // -- Check if track is accepted
 590   // > for basic parameters
 591
 592   if (!track)
 593     return kFALSE;
 594
 595   if (track->Charge() == 0)
 596     return kFALSE;
 597
 598   return kTRUE;
 599 }
 600
 601 //________________________________________________________________________
 602 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP) {
 603   // -- Check if track is accepted
 604   // > in rapidity
 605   // > if no pid : return kTRUE
 606   // > return 0 if not accepted
 607
 608   if (!fUsePID) {
 609     yP = track->Eta();
 610     return kTRUE;
 611   }
 612
 613   Double_t mP = AliPID::ParticleMass(fParticleSpecies);
 614
 615   // -- Calculate rapidities and kinematics
 616   Double_t pvec[3];
 617   track->GetPxPyPz(pvec);
 618
 619   Double_t p  = track->P();
 620   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 621            yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
 622
 623   // -- Check Rapidity window
 624   if (TMath::Abs(yP) > fRapidityMax)
 625     return kFALSE;
 626
 627   return kTRUE;
 628 }
 629
 630 //________________________________________________________________________
 631 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedDCA(AliVTrack *vTrack) {
 632   // -- Check if track is accepted - ONLY FOR ESDs so far
 633   // > for DCA, if both SPD layers have hits
 634   // > For now only Implemented for ESDs
 635
 636   Bool_t isAccepted = kTRUE;
 637
 638   if (!fESD)
 639     return isAccepted;
 640
 641   AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
 642   if (!track)
 643     return kFALSE;
 644
 645   // -- Get nHits SPD
 646   if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
 647
 648     // -- Get DCA nSigmas
 649     Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 650     track->GetImpactParameters(dca,cov);
 651
 652     Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
 653     Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
 654
 655     if (fNSigmaMaxCdd != 0.) {
 656       if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
 657         isAccepted = kFALSE;
 658     }
 659
 660     if (fNSigmaMaxCzz != 0.) {
 661       if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
 662         isAccepted = kFALSE;
 663     }
 664   }
 665
 666   return isAccepted;
 667 }
 668
 669 //________________________________________________________________________
 670 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid) {
 671   // -- Check if track is accepted
 672   // > provides TPC and TOF nSigmas to the argument
 673
 674   Bool_t isAcceptedTPC = kFALSE;
 675   Bool_t isAcceptedTOF = kTRUE;
 676
 677   // -- In case not PID is used
 678   if (!fUsePID) {
 679     pid[0] = 10.;
 680     pid[1] = 10.;
 681     return kTRUE;
 682   }
 683
 684   // -- Get PID
 685   pid[0] = fPIDResponse->NumberOfSigmasTPC(track, fParticleSpecies);
 686   pid[1] = fPIDResponse->NumberOfSigmasTOF(track, fParticleSpecies);
 687
 688   // -- Check PID with TPC
 689   if (TMath::Abs(pid[0]) < fNSigmaMaxTPC)
 690     isAcceptedTPC = kTRUE;
 691
 692   // -- Check PID with TOF
 693   if (isAcceptedTPC) {
 694
 695     // Has TOF PID availible
 696     if (fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track) == AliPIDResponse::kDetPidOk) {
 697       if (TMath::Abs(pid[1]) < fNSigmaMaxTOF)
 698         isAcceptedTOF = kTRUE;
 699       else
 700         isAcceptedTOF = kFALSE;
 701     }
 702
 703     // Has no TOF PID availible
 704     else {
 705       if (track->Pt() >= fMinPtForTOFRequired)
 706         isAcceptedTOF = kFALSE;
 707       else
 708         isAcceptedTOF = kTRUE;
 709     }
 710   } // if (isAcceptedTPC) {
 711
 712   return (isAcceptedTPC && isAcceptedTOF);
 713 }
 714
 715 //________________________________________________________________________
 716 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPhi(AliVTrack *track) {
 717   // -- Check if track is accepted
 718   // > in phi
 719   // > return 0 if not accepted
 720
 721   if (track->Phi() > fPhiMin && track->Phi() <= fPhiMax)
 722     return kTRUE;
 723   else if (track->Phi() < fPhiMin && (track->Phi() + TMath::TwoPi()) <= fPhiMax)
 724     return kTRUE;
 725   else
 726     return kFALSE;
 727 }
 728
 729 /*
 730  * ---------------------------------------------------------------------------------
 731  *                         Helper Methods
 732  * ---------------------------------------------------------------------------------
 733  */
 734
 735 //________________________________________________________________________
 736 void AliAnalysisNetParticleHelper::UpdateEtaCorrectedTPCPid() {
 737   // -- Update eta corrected TPC pid
 738
 739   if (!fEtaCorrFunc)
 740     return;
 741
 742   for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 743     AliESDtrack *track = fESD->GetTrack(idxTrack);
 744
 745     // -- Check if charged track is accepted for basic parameters
 746     if (!IsTrackAcceptedBasicCharged(track))
 747       continue;
 748
 749     Double_t newTPCSignal = track->GetTPCsignal() / fEtaCorrFunc->Eval(track->Eta());
 750     track->SetTPCsignal(newTPCSignal, track->GetTPCsignalSigma(), track->GetTPCsignalN());
 751   }
 752 }
 753
 754 //________________________________________________________________________
 755 void AliAnalysisNetParticleHelper::BinLogAxis(const THnBase *hn, Int_t axisNumber) {
 756   // -- Method for the correct logarithmic binning of histograms
 757   // -- and update fMinPtForTOFRequired using the logarithmic scale
 758
 759   // -- Make logarithmic binning
 760   TAxis *axis = hn->GetAxis(axisNumber);
 761   Int_t  nBins = axis->GetNbins();
 762
 763   Double_t from  = axis->GetXmin();
 764   Double_t to    = axis->GetXmax();
 765   Double_t *newBins = new Double_t[nBins + 1];
 766
 767   newBins[0] = from;
 768   Double_t factor = TMath::Power(to/from, 1./nBins);
 769
 770   for (int ii = 1; ii <= nBins; ii++)
 771    newBins[ii] = factor * newBins[ii-1];
 772
 773   axis->Set(nBins, newBins);
 774
 775   delete [] newBins;
 776
 777   // -- Update MinPtForTOFRequired
 778   for (Int_t ii = 1; ii <= nBins; ++ii)
 779     if (axis->GetBinLowEdge(ii) <= fMinPtForTOFRequired && axis->GetBinLowEdge(ii) + axis->GetBinWidth(ii) >  fMinPtForTOFRequired)
 780       fMinPtForTOFRequired = axis->GetBinLowEdge(ii) + axis->GetBinWidth(ii);
 781 }
 782
 783 ///////////////////////////////////////////////////////////////////////////////////
 784
 785 /*
 786  * ---------------------------------------------------------------------------------
 787  *                           Initialize - Private
 788  * ---------------------------------------------------------------------------------
 789  */
 790
 791 //________________________________________________________________________
 792 void AliAnalysisNetParticleHelper::InitializeEventStats() {
 793   // -- Initialize event statistics histograms
 794
 795   fHEventStat0 = new TH1F("fHEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 796   fHEventStat1 = new TH1F("fHEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 797
 798   for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
 799     fHEventStat0->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkEventNames[ii]);
 800     fHEventStat1->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkEventNames[ii]);
 801   }
 802
 803   fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliAnalysisNetParticleHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
 804   fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliAnalysisNetParticleHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
 805 }
 806
 807 //________________________________________________________________________
 808 void AliAnalysisNetParticleHelper::InitializeTriggerStats() {
 809   // -- Initialize trigger statistics histograms
 810
 811   fHTriggerStat = new TH1F("fHTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
 812
 813   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 814     fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkTriggerNames[ii]);
 815 }
 816
 817 //________________________________________________________________________
 818 void AliAnalysisNetParticleHelper::InitializeCentralityStats() {
 819   // -- Initialize trigger statistics histograms
 820
 821   fHCentralityStat = new TH1F("fHCentralityStat","Centrality statistics;Centrality Bins;Events",
 822                               fNCentralityBins,-0.5,fNCentralityBins-0.5);
 823
 824   for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
 825     fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, AliAnalysisNetParticleHelper::fgkCentralityNames[ii]);
 826 }
 827
 828 //________________________________________________________________________
 829 Int_t AliAnalysisNetParticleHelper::InitializeEtaCorrection(Bool_t isMC) {
 830   // -- Initialize eta correction maps for TPC pid
 831
 832   TFile fileEtaCorrMaps("${ALICE_ROOT}/PWGDQ/dielectron/files/EtaCorrMaps.root");
 833   if (!fileEtaCorrMaps.IsOpen())
 834     return -1;
 835
 836   TList *keyList = fileEtaCorrMaps.GetListOfKeys();
 837
 838   TString sList("");
 839   sList = (isMC) ? "LHC11a10" :  "LHC10h.pass2";
 840
 841   for (Int_t idx = 0; idx < keyList->GetEntries(); ++idx) {
 842     TString keyName = keyList->At(idx)->GetName();
 843     TPRegexp reg(keyName);
 844     if (reg.MatchB(sList)){
 845       AliInfo(Form("Using Eta Correction Function: %s",keyName.Data()));
 846       fEtaCorrFunc = static_cast<TF1*>(fileEtaCorrMaps.Get(keyName.Data()));
 847       return 0;
 848     }
 849   }
 850
 851   return -2;
 852 }
 853
 854 /*
 855  * ---------------------------------------------------------------------------------
 856  *                         Event / Trigger Statistics - private
 857  * ---------------------------------------------------------------------------------
 858  */
 859
 860 //________________________________________________________________________
 861 Bool_t AliAnalysisNetParticleHelper::FillEventStats(Int_t *aEventCuts) {
 862   // -- Fill event / centrality statistics
 863
 864   Bool_t isRejected = kFALSE;
 865
 866   // -- Fill event statistics
 867   for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
 868
 869     if (aEventCuts[idx])
 870       isRejected = kTRUE;
 871     else
 872       fHEventStat0->Fill(idx);
 873   }
 874
 875   for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
 876     if (aEventCuts[idx])
 877       break;
 878     fHEventStat1->Fill(idx);
 879   }
 880
 881   // -- Fill centrality statistics of accepted events
 882   if (!isRejected)
 883     fHCentralityStat->Fill(fCentralityBin);
 884
 885   return isRejected;
 886 }