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 #include "TDatabasePDG.h"
   9
  10 #include "AliStack.h"
  11 #include "AliMCEvent.h"
  12 #include "AliMCParticle.h"
  13 #include "AliESDtrackCuts.h"
  14 #include "AliESDInputHandler.h"
  15 #include "AliESDpid.h"
  16 #include "AliAODInputHandler.h"
  17 #include "AliAODEvent.h"
  18 #include "AliAODMCParticle.h"
  19 #include "AliCentrality.h"
  20 #include "AliTracker.h"
  21
  22 #include "AliAnalysisNetParticleHelper.h"
  23
  24 using namespace std;
  25
  26 // Task for NetParticle checks
  27 // Author: Jochen Thaeder <jochen@thaeder.de>
  28
  29 ClassImp(AliAnalysisNetParticleHelper)
  30
  31 /*
  32  * ---------------------------------------------------------------------------------
  33  *                            particle names
  34  * ---------------------------------------------------------------------------------
  35  */
  36
  37   const Char_t* aPartNames[AliPID::kSPECIES][2]             = {
  38     {"ele",     "posi"},
  39     {"mubar",   "mu"},
  40     {"pibar",   "pi"},
  41     {"kbar",    "k"},
  42     {"pbar",    "p"}
  43   };
  44
  45   const Char_t* aPartTitles[AliPID::kSPECIES][2]            = {
  46     {"Electron",    "Positron"},
  47     {"Anti-Muon",   "Muon"},
  48     {"Anti-Pion",   "Proton"},
  49     {"Anti-Kaon",   "Kaon"},
  50     {"Anti-Proton", "Proton"}
  51   };
  52
  53   const Char_t* aPartTitlesLatex[AliPID::kSPECIES][2]        = {
  54     {"e^{-}",   "e^{+}" },
  55     {"#mu^{-}", "#mu^{+}" },
  56     {"#pi^{-}", "#pi^{+}" },
  57     {"K^{-}",   "K^{+}" },
  58     {"#bar{p}", "p"}
  59   };
  60
  61   const Char_t* aControlPartNames[2]      = {"lambdabar",     "lambda"};
  62   const Char_t* aControlPartTitles[2]     = {"Anti-Lambda",   "Lambda"};
  63
  64 /*
  65  * ---------------------------------------------------------------------------------
  66  *                            Constructor / Destructor
  67  * ---------------------------------------------------------------------------------
  68  */
  69
  70 //________________________________________________________________________
  71 AliAnalysisNetParticleHelper::AliAnalysisNetParticleHelper() :
  72   fESDHandler(NULL),
  73   fPIDResponse(NULL),
  74   fESD(NULL),
  75   fAODHandler(NULL),
  76   fAOD(NULL),
  77   fMCEvent(NULL),
  78   fStack(NULL),
  79
  80   fCentralityBin(-1),
  81   fCentralityPercentile(-1.),
  82
  83   fCentralityBinMax(7),
  84   fVertexZMax(10.),
  85   fRapidityMax(0.5),
  86   fMinTrackLengthMC(70.),
  87   fNSigmaMaxCdd(3.),
  88   fNSigmaMaxCzz(3.),
  89
  90   fParticleSpecies(AliPID::kProton),
  91   fControlParticleCode(3122),
  92   fControlParticleIsNeutral(kTRUE),
  93   fControlParticleName("Lambda"),
  94
  95   fNSigmaMaxTPC(2.5),
  96   fNSigmaMaxTOF(2.5),
  97   fMinPtForTOFRequired(0.8),
  98
  99   fHEventStat0(NULL),
 100   fHEventStat1(NULL),
 101   fHEventStatMax(6),
 102
 103   fHTriggerStat(NULL),
 104   fNTriggers(5),
 105
 106   fHCentralityStat(NULL),
 107   fNCentralityBins(10),
 108
 109   fEtaCorrFunc(NULL),
 110   fCorr0(NULL),
 111   fCorr1(NULL) {
 112   // Constructor
 113
 114   AliLog::SetClassDebugLevel("AliAnalysisNetParticleHelper",10);
 115 }
 116
 117 //________________________________________________________________________
 118 AliAnalysisNetParticleHelper::~AliAnalysisNetParticleHelper() {
 119   // Destructor
 120
 121   for (Int_t jj = 0; jj < 2; ++jj) {
 122     if (fCorr0[jj]) delete[] fCorr0[jj];
 123     if (fCorr1[jj]) delete[] fCorr1[jj];
 124   }
 125   if (fCorr0) delete[] fCorr0;
 126   if (fCorr1) delete[] fCorr1;
 127
 128   return;
 129 }
 130
 131 /*
 132  * ---------------------------------------------------------------------------------
 133  *                                    Setter
 134  * ---------------------------------------------------------------------------------
 135  */
 136
 137 //________________________________________________________________________
 138 void AliAnalysisNetParticleHelper::SetParticleSpecies(AliPID::EParticleType pid) {
 139   // -- Set particle species (ID, Name, Title, Title LATEX)
 140
 141   if( (Int_t)pid < 0 || (Int_t)pid > AliPID::kSPECIES){
 142     AliWarning("Particle ID not in AliPID::kSPECIES --> Set to protons");
 143     pid = AliPID::kProton;
 144   }
 145
 146   fParticleSpecies     = pid;
 147
 148   for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
 149     fPartName[idxPart]       = aPartNames[fParticleSpecies][idxPart];
 150     fPartTitle[idxPart]      = aPartTitles[fParticleSpecies][idxPart];
 151     fPartTitleLatex[idxPart] = aPartTitlesLatex[fParticleSpecies][idxPart];
 152   }
 153 }
 154
 155
 156 /*
 157  * ---------------------------------------------------------------------------------
 158  *                                    Getter
 159  * ---------------------------------------------------------------------------------
 160  */
 161
 162 //________________________________________________________________________
 163 TString AliAnalysisNetParticleHelper::GetParticleName(Int_t idxPart) {
 164   // -- Get particle Name
 165
 166   if( idxPart != 0 && idxPart != 1){
 167     AliWarning("Particle type not known --> Set to antiparticles");
 168     idxPart = 0;
 169   }
 170
 171   return fPartName[idxPart];
 172 }
 173
 174 //________________________________________________________________________
 175 TString AliAnalysisNetParticleHelper::GetParticleTitle(Int_t idxPart) {
 176   // -- Get particle Title
 177
 178   if( idxPart != 0 && idxPart != 1){
 179     AliWarning("Particle type not known --> Set to antiparticles");
 180     idxPart = 0;
 181   }
 182
 183   return fPartTitle[idxPart];
 184 }
 185
 186 //________________________________________________________________________
 187 TString AliAnalysisNetParticleHelper::GetParticleTitleLatex(Int_t idxPart) {
 188   // -- Get particle Title LATEX
 189
 190   if( idxPart != 0 && idxPart != 1){
 191     AliWarning("Particle type not known --> Set to antiparticles");
 192     idxPart = 0;
 193   }
 194
 195   return fPartTitleLatex[idxPart];
 196 }
 197
 198 //________________________________________________________________________
 199 TString AliAnalysisNetParticleHelper::GetControlParticleName(Int_t idxPart) {
 200   // -- Get particle Title LATEX
 201
 202   if( idxPart != 0 && idxPart != 1){
 203     AliWarning("Particle type not known --> Set to antiparticles");
 204     idxPart = 0;
 205   }
 206
 207   return aControlPartNames[idxPart];
 208 }
 209
 210 //________________________________________________________________________
 211 TString AliAnalysisNetParticleHelper::GetControlParticleTitle(Int_t idxPart) {
 212   // -- Get particle Title LATEX
 213
 214   if( idxPart != 0 && idxPart != 1){
 215     AliWarning("Particle type not known --> Set to antiparticles");
 216     idxPart = 0;
 217   }
 218
 219   return aControlPartTitles[idxPart];
 220 }
 221
 222 /*
 223  * ---------------------------------------------------------------------------------
 224  *                                 Public Methods
 225  * ---------------------------------------------------------------------------------
 226  */
 227
 228 //________________________________________________________________________
 229 Int_t AliAnalysisNetParticleHelper::Initialize(Bool_t isMC) {
 230   // -- Initialize helper
 231
 232   Int_t iResult = 0;
 233
 234   // -- Setup event cut statistics
 235   InitializeEventStats();
 236
 237   // -- Setup trigger statistics
 238   InitializeTriggerStats();
 239
 240   // -- Setup centrality statistics
 241   InitializeCentralityStats();
 242
 243   // -- Load Eta correction function
 244   iResult = InitializeEtaCorrection(isMC);
 245
 246   // -- Load track by track correction function
 247   iResult = InitializeTrackbyTrackCorrection();
 248
 249   return iResult;
 250 }
 251
 252 //________________________________________________________________________
 253 Int_t AliAnalysisNetParticleHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
 254   // -- Setup Event
 255
 256   // -- Get ESD objects
 257   if(esdHandler){
 258     fESDHandler  = esdHandler;
 259     fPIDResponse = esdHandler->GetPIDResponse();
 260     fESD         = fESDHandler->GetEvent();
 261   }
 262
 263   // -- Get AOD objects
 264   else if(aodHandler){
 265     fAODHandler  = aodHandler;
 266     fPIDResponse = aodHandler->GetPIDResponse();
 267     fAOD         = fAODHandler->GetEvent();
 268   }
 269
 270   // -- Get MC objects
 271   fMCEvent     = mcEvent;
 272   if (fMCEvent)
 273     fStack     = fMCEvent->Stack();
 274
 275   // -- Get event centrality
 276   // >  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 277   // >   0   1     2     3     4     5     6     7     8     9
 278
 279   AliCentrality *centrality = NULL;
 280
 281   if(esdHandler)
 282     centrality = fESD->GetCentrality();
 283   else if(aodHandler)
 284     centrality = fAOD->GetHeader()->GetCentralityP();
 285
 286   Int_t centBin = centrality->GetCentralityClass10("V0M");
 287   if (centBin == 0)
 288     fCentralityBin = centrality->GetCentralityClass5("V0M");
 289   else if (centBin == 10 || centBin == -1.)
 290     fCentralityBin = -1;
 291   else if (centBin > 0 && centBin < fNCentralityBins)
 292     fCentralityBin = centBin + 1;
 293   else
 294     fCentralityBin = -2;
 295
 296   // -- Stay within the max centrality bin
 297   if (fCentralityBin >= fCentralityBinMax)
 298     fCentralityBin = -2;
 299
 300   fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
 301
 302   if(esdHandler){
 303     // -- Update TPC pid with eta correction (only for ESDs?) XXX
 304     UpdateEtaCorrectedTPCPid();
 305   }
 306
 307   return 0;
 308 }
 309
 310 /*
 311  * ---------------------------------------------------------------------------------
 312  *                         Event / Trigger Statistics
 313  * ---------------------------------------------------------------------------------
 314  */
 315
 316 //________________________________________________________________________
 317 Bool_t AliAnalysisNetParticleHelper::IsEventTriggered() {
 318   // -- Check if Event is triggered and fill Trigger Histogram
 319
 320   Bool_t *aTriggerFired = new Bool_t[fNTriggers];
 321   for (Int_t ii = 0; ii < fNTriggers; ++ii)
 322     aTriggerFired[ii] = kFALSE;
 323
 324   if(fESDHandler){
 325     if ((fESDHandler->IsEventSelected() & AliVEvent::kMB))          aTriggerFired[0] = kTRUE;
 326     if ((fESDHandler->IsEventSelected() & AliVEvent::kCentral))     aTriggerFired[1] = kTRUE;
 327     if ((fESDHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
 328     if ((fESDHandler->IsEventSelected() & AliVEvent::kEMCEJE))      aTriggerFired[3] = kTRUE;
 329     if ((fESDHandler->IsEventSelected() & AliVEvent::kEMCEGA))      aTriggerFired[4] = kTRUE;
 330   }
 331   else if(fAODHandler){
 332     if ((fAODHandler->IsEventSelected() & AliVEvent::kMB))          aTriggerFired[0] = kTRUE;
 333     if ((fAODHandler->IsEventSelected() & AliVEvent::kCentral))     aTriggerFired[1] = kTRUE;
 334     if ((fAODHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
 335     if ((fAODHandler->IsEventSelected() & AliVEvent::kEMCEJE))      aTriggerFired[3] = kTRUE;
 336     if ((fAODHandler->IsEventSelected() & AliVEvent::kEMCEGA))      aTriggerFired[4] = kTRUE;
 337   }
 338
 339   Bool_t isTriggered = kFALSE;
 340
 341   for (Int_t ii=0; ii<fNTriggers; ++ii) {
 342     if(aTriggerFired[ii]) {
 343       isTriggered = kTRUE;
 344       fHTriggerStat->Fill(ii);
 345     }
 346   }
 347
 348   delete[] aTriggerFired;
 349
 350   return isTriggered;
 351 }
 352
 353 //________________________________________________________________________
 354 Bool_t AliAnalysisNetParticleHelper::IsEventRejected() {
 355   // -- Evaluate event statistics histograms
 356
 357   Int_t *aEventCuts = new Int_t[fHEventStatMax];
 358   // set aEventCuts[ii] to 1 in case of reject
 359
 360   for (Int_t ii=0;ii<fHEventStatMax; ++ii)
 361     aEventCuts[ii] = 0;
 362
 363   Int_t iCut = 0;
 364
 365   // -- 0 - Before Physics Selection
 366   aEventCuts[iCut] = 0;
 367
 368   // -- 1 - No Trigger fired
 369   ++iCut;
 370   if (!IsEventTriggered())
 371     aEventCuts[iCut] = 1;
 372
 373   // -- 2 - No Vertex
 374   ++iCut;
 375   const AliESDVertex* vtxESD = NULL;
 376   const AliAODVertex* vtxAOD = NULL;
 377   if (fESD){
 378     vtxESD = fESD->GetPrimaryVertexTracks();
 379     if (!vtxESD)
 380       aEventCuts[iCut] = 1;
 381   }
 382   else if (fAOD){
 383     vtxAOD = fAOD->GetPrimaryVertex();
 384     if (!vtxAOD)
 385       aEventCuts[iCut] = 1;
 386   }
 387
 388   // -- 3 - Vertex z outside cut window
 389   ++iCut;
 390   if (vtxESD){
 391     if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax)
 392       aEventCuts[iCut] = 1;
 393   }
 394   else if(vtxAOD){
 395     if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax)
 396       aEventCuts[iCut] = 1;
 397   }
 398   else
 399     aEventCuts[iCut] = 1;
 400
 401   // -- 4 - Centrality = -1  (no centrality or not hadronic)
 402   ++iCut;
 403   if(fCentralityBin == -1.)
 404     aEventCuts[iCut] = 1;
 405
 406   // -- 5 - Centrality < fCentralityMax
 407   ++iCut;
 408   if(fCentralityBin == -2.)
 409     aEventCuts[iCut] = 1;
 410
 411   // -- Fill statistics / reject event
 412   Bool_t isRejected = FillEventStats(aEventCuts);
 413
 414   // -- Cleanup
 415   delete[] aEventCuts;
 416
 417   return isRejected;
 418 }
 419
 420 /*
 421  * ---------------------------------------------------------------------------------
 422  *                         Accept Particle Methods - private
 423  * ---------------------------------------------------------------------------------
 424  */
 425
 426 //________________________________________________________________________
 427 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(TParticle *particle, Int_t idxMC) {
 428   // -- Check if MC particle is accepted for basic parameters
 429
 430   if (!particle)
 431     return kFALSE;
 432
 433   // -- check if PDF code exists
 434   if (!particle->GetPDG())
 435     return kFALSE;
 436
 437   // -- check if charged
 438   if (particle->GetPDG()->Charge() == 0.0)
 439     return kFALSE;
 440
 441   // -- check if physical primary
 442   if(!fStack->IsPhysicalPrimary(idxMC))
 443     return kFALSE;
 444
 445   return kTRUE;
 446 }
 447 //________________________________________________________________________
 448 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicCharged(AliAODMCParticle *particle) {
 449   // -- Check if MC particle is accepted for basic parameters
 450
 451   if (!particle)
 452     return kFALSE;
 453
 454   // -- check if PDF code exists
 455   if (!(TDatabasePDG::Instance()->GetParticle(particle->PdgCode())))
 456     return kFALSE;
 457
 458   // -- check if charged
 459   if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() == 0.0)
 460     return kFALSE;
 461
 462   // -- check if physical primary
 463   if(!particle->IsPhysicalPrimary())
 464     return kFALSE;
 465
 466   return kTRUE;
 467 }
 468 //________________________________________________________________________
 469 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(TParticle *particle, Int_t idxMC) {
 470   // -- Check if MC particle is accepted for basic parameters
 471
 472   if (!particle)
 473     return kFALSE;
 474
 475   // -- check if PDF code exists
 476   if (!particle->GetPDG())
 477     return kFALSE;
 478
 479   // -- check if neutral
 480   if (particle->GetPDG()->Charge() != 0.0)
 481     return kFALSE;
 482
 483   // -- check if physical primary
 484   if(!fStack->IsPhysicalPrimary(idxMC))
 485     return kFALSE;
 486
 487   return kTRUE;
 488 }
 489 //________________________________________________________________________
 490 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedBasicNeutral(AliAODMCParticle *particle) {
 491   // -- Check if MC particle is accepted for basic parameters
 492
 493   if (!particle)
 494     return kFALSE;
 495
 496   // -- check if PDF code exists
 497   if (!(TDatabasePDG::Instance()->GetParticle(particle->PdgCode())))
 498     return kFALSE;
 499
 500   // -- check if charged
 501   if ((TDatabasePDG::Instance()->GetParticle(particle->PdgCode()))->Charge() != 0.0)
 502     return kFALSE;
 503
 504   // -- check if physical primary
 505   if(!particle->IsPhysicalPrimary())
 506     return kFALSE;
 507
 508   return kTRUE;
 509 }
 510 //________________________________________________________________________
 511 Bool_t AliAnalysisNetParticleHelper::IsParticleAcceptedRapidity(TParticle *particle, Double_t &yP) {
 512   // -- Check if particle is accepted
 513   // > in rapidity
 514   // > return 0 if not accepted
 515
 516   Double_t mP = AliPID::ParticleMass(fParticleSpecies);
 517
 518   // -- Calculate rapidities and kinematics
 519   Double_t p  = particle->P();
 520   Double_t pz = particle->Pz();
 521
 522   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 523   yP          = 0.5 * TMath::Log((eP + pz) / (eP - pz));
 524
 525   // -- Check Rapidity window
 526   if (TMath::Abs(yP) > fRapidityMax)
 527     return kFALSE;
 528
 529   return kTRUE;
 530 }
 531
 532 //_____________________________________________________________________________
 533 Bool_t AliAnalysisNetParticleHelper::IsParticleFindable(Int_t label) {
 534   // -- Check if MC particle is findable tracks
 535
 536   AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
 537   if(!mcParticle)
 538     return kFALSE;
 539
 540   Int_t counter;
 541   Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0);
 542
 543   return (tpcTrackLength > fMinTrackLengthMC);
 544 }
 545
 546 /*
 547  * ---------------------------------------------------------------------------------
 548  *                            Accept Track Methods - public
 549  * ---------------------------------------------------------------------------------
 550  */
 551
 552 //________________________________________________________________________
 553 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliESDtrack *track) {
 554   // -- Check if track is accepted
 555   // > for basic parameters
 556
 557   if (!track)
 558     return kFALSE;
 559
 560   if (track->Charge() == 0)
 561     return kFALSE;
 562
 563   // -- Get momentum for dEdx
 564   if (!track->GetInnerParam())
 565     return kFALSE;
 566
 567   return kTRUE;
 568 }
 569
 570 //________________________________________________________________________
 571 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedBasicCharged(AliAODTrack *track) {
 572   // -- Check if track is accepted
 573   // > for basic parameters
 574
 575   if (!track)
 576     return kFALSE;
 577
 578   if (track->Charge() == 0)
 579     return kFALSE;
 580
 581   //// -- Get momentum for dEdx --> returns always ZERO XXX
 582   //if (!track->GetInnerParam())
 583   //  return kFALSE;
 584
 585   return kTRUE;
 586 }
 587
 588 //________________________________________________________________________
 589 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP) {
 590   // -- Check if track is accepted
 591   // > in rapidity
 592   // > return 0 if not accepted
 593
 594   Double_t mP = AliPID::ParticleMass(fParticleSpecies);
 595
 596   // -- Calculate rapidities and kinematics
 597   Double_t pvec[3];
 598   track->GetPxPyPz(pvec);
 599
 600   Double_t p  = track->P();
 601   Double_t eP = TMath::Sqrt(p*p + mP*mP);
 602            yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
 603
 604   // -- Check Rapidity window
 605   if (TMath::Abs(yP) > fRapidityMax)
 606     return kFALSE;
 607
 608   return kTRUE;
 609 }
 610
 611 //________________________________________________________________________
 612 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedDCA(AliESDtrack *track) {  //ONLY FOR ESDs so far XXX
 613   // -- Check if track is accepted
 614   // > for DCA, if both SPD layers have hits
 615
 616   Bool_t isAccepted = kTRUE;
 617
 618   // -- Get nHits SPD
 619   if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
 620
 621     // -- Get DCA nSigmas
 622     Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 623     track->GetImpactParameters(dca,cov);
 624
 625     Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
 626     Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
 627
 628     if (fNSigmaMaxCdd != 0.) {
 629       if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
 630         isAccepted = kFALSE;
 631     }
 632
 633     if (fNSigmaMaxCzz != 0.) {
 634       if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
 635         isAccepted = kFALSE;
 636     }
 637   }
 638
 639   return isAccepted;
 640 }
 641
 642 //________________________________________________________________________
 643 Bool_t AliAnalysisNetParticleHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid) {
 644   // -- Check if track is accepted
 645   // > provides TPC and TOF nSigmas to the argument
 646
 647   Bool_t isAcceptedTPC = kFALSE;
 648   Bool_t isAcceptedTOF = kTRUE;
 649
 650   // -- Get PID
 651   pid[0] = fPIDResponse->NumberOfSigmasTPC(track, fParticleSpecies);
 652   pid[1] = fPIDResponse->NumberOfSigmasTOF(track, fParticleSpecies);
 653
 654   // -- Check PID with TPC
 655   if (TMath::Abs(pid[0]) < fNSigmaMaxTPC)
 656     isAcceptedTPC = kTRUE;
 657
 658   // -- Check PID with TOF
 659   if (isAcceptedTPC) {
 660
 661     // Has TOF PID availible
 662     if (track->GetStatus() & AliVTrack::kTOFpid) {
 663       if (TMath::Abs(pid[1]) < fNSigmaMaxTOF)
 664         isAcceptedTOF = kTRUE;
 665       else
 666         isAcceptedTOF = kFALSE;
 667     }
 668
 669     // Has no TOF PID availible
 670     else {
 671       if (track->Pt() > fMinPtForTOFRequired)
 672         isAcceptedTOF = kFALSE;
 673       else
 674         isAcceptedTOF = kTRUE;
 675     }
 676   } // if (isAcceptedTPC) {
 677
 678   return (isAcceptedTPC && isAcceptedTOF);
 679 }
 680
 681 /*
 682  * ---------------------------------------------------------------------------------
 683  *                         Helper Methods
 684  * ---------------------------------------------------------------------------------
 685  */
 686
 687 //________________________________________________________________________
 688 void AliAnalysisNetParticleHelper::UpdateEtaCorrectedTPCPid() {
 689   // -- Update eta corrected TPC pid
 690
 691   if (!fEtaCorrFunc)
 692     return;
 693
 694   for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 695     AliESDtrack *track = fESD->GetTrack(idxTrack);
 696
 697     // -- Check if charged track is accepted for basic parameters
 698     if (!IsTrackAcceptedBasicCharged(track))
 699       continue;
 700
 701     Double_t newTPCSignal = track->GetTPCsignal() / fEtaCorrFunc->Eval(track->Eta());
 702     track->SetTPCsignal(newTPCSignal, track->GetTPCsignalSigma(), track->GetTPCsignalN());
 703   }
 704 }
 705
 706 //________________________________________________________________________
 707 Double_t AliAnalysisNetParticleHelper::GetTrackbyTrackCorrectionFactor(Double_t *aTrack,  Int_t flag) {
 708   // -- Get efficiency correctionf of particle dependent on (eta, phi, pt, centrality)
 709
 710   Int_t idxPart = (aTrack[2] < 0) ? 0 : 1;
 711   THnSparseF* corrMatrix = (flag == 0) ? fCorr0[idxPart][fCentralityBin] : fCorr1[idxPart][fCentralityBin];
 712
 713   Double_t dimBin[3] = {aTrack[3], aTrack[4], aTrack[1]}; // eta, phi, pt
 714
 715   Double_t corr = corrMatrix->GetBinContent(corrMatrix->GetBin(dimBin));
 716   if (corr == 0.) {
 717     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",
 718                   aTrack[3], aTrack[4], aTrack[1], fCentralityBin));
 719     corr = 1.;
 720   }
 721
 722   return corr;
 723 }
 724
 725 //________________________________________________________________________
 726 void AliAnalysisNetParticleHelper::BinLogAxis(const THnSparseF *h, Int_t axisNumber) {
 727   // -- Method for the correct logarithmic binning of histograms
 728
 729   TAxis *axis = h->GetAxis(axisNumber);
 730   Int_t  bins = axis->GetNbins();
 731
 732   Double_t from  = axis->GetXmin();
 733   Double_t to    = axis->GetXmax();
 734   Double_t *newBins = new Double_t[bins + 1];
 735
 736   newBins[0] = from;
 737   Double_t factor = pow(to/from, 1./bins);
 738
 739   for (int i = 1; i <= bins; i++) {
 740    newBins[i] = factor * newBins[i-1];
 741   }
 742   axis->Set(bins, newBins);
 743   delete [] newBins;
 744 }
 745
 746 ///////////////////////////////////////////////////////////////////////////////////
 747
 748 /*
 749  * ---------------------------------------------------------------------------------
 750  *                           Initialize - Private
 751  * ---------------------------------------------------------------------------------
 752  */
 753
 754 //________________________________________________________________________
 755 void AliAnalysisNetParticleHelper::InitializeEventStats() {
 756   // -- Initialize event statistics histograms
 757
 758   const Char_t* aEventNames[]      = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,90]%"};
 759   const Char_t* aCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
 760
 761   fHEventStat0 = new TH1F("fHEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 762   fHEventStat1 = new TH1F("fHEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
 763
 764   for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
 765     fHEventStat0->GetXaxis()->SetBinLabel(ii+1, aEventNames[ii]);
 766     fHEventStat1->GetXaxis()->SetBinLabel(ii+1, aEventNames[ii]);
 767   }
 768   fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", aCentralityMaxNames[fCentralityBinMax-1]));
 769   fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", aCentralityMaxNames[fCentralityBinMax-1]));
 770 }
 771
 772 //________________________________________________________________________
 773 void AliAnalysisNetParticleHelper::InitializeTriggerStats() {
 774   // -- Initialize trigger statistics histograms
 775
 776   const Char_t* aTriggerNames[] = { "kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" };
 777
 778   fHTriggerStat = new TH1F("fHTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
 779
 780   for ( Int_t ii=0; ii < fNTriggers; ii++ )
 781     fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, aTriggerNames[ii]);
 782 }
 783
 784 //________________________________________________________________________
 785 void AliAnalysisNetParticleHelper::InitializeCentralityStats() {
 786   // -- Initialize trigger statistics histograms
 787
 788   const Char_t* aCentralityNames[] = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%",
 789                                       "50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
 790
 791   fHCentralityStat = new TH1F("fHCentralityStat","Centrality statistics;Centrality Bins;Events",
 792                               fNCentralityBins,-0.5,fNCentralityBins-0.5);
 793
 794   for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
 795     fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, aCentralityNames[ii]);
 796 }
 797
 798 //________________________________________________________________________
 799 Int_t AliAnalysisNetParticleHelper::InitializeEtaCorrection(Bool_t isMC) {
 800   // -- Initialize eta correction maps for TPC pid
 801
 802   TFile fileEtaCorrMaps("${ALICE_ROOT}/PWGDQ/dielectron/files/EtaCorrMaps.root");
 803   if (!fileEtaCorrMaps.IsOpen())
 804     return -1;
 805
 806   TList *keyList = fileEtaCorrMaps.GetListOfKeys();
 807
 808   TString sList("");
 809   sList = (isMC) ? "LHC11a10" :  "LHC10h.pass2";
 810
 811   for (Int_t idx = 0; idx < keyList->GetEntries(); ++idx) {
 812     TString keyName = keyList->At(idx)->GetName();
 813     TPRegexp reg(keyName);
 814     if (reg.MatchB(sList)){
 815       AliInfo(Form("Using Eta Correction Function: %s",keyName.Data()));
 816       fEtaCorrFunc = static_cast<TF1*>(fileEtaCorrMaps.Get(keyName.Data()));
 817       return 0;
 818     }
 819   }
 820
 821   return -2;
 822 }
 823
 824 //________________________________________________________________________
 825 Int_t AliAnalysisNetParticleHelper::InitializeTrackbyTrackCorrection() {
 826   // -- Initialize track by track correction matrices
 827
 828   if( fParticleSpecies != AliPID::kProton){
 829     AliWarning("Only proton corrections (preliminary available) --> Will correct with these");
 830   }
 831
 832   // In future --> !!!!!!!!!!!!!!!!!!!!!!!
 833   //TFile* corrFile  = TFile::Open("${ALICE_ROOT}/PWGCF/EBYE/NetParticle/eff/effectiveCorrection.root");
 834   TFile* corrFile  = TFile::Open("${ALICE_ROOT}/PWGCF/EBYE/NetParticle/eff/effectiveCorrectionProton.root");
 835   if (!corrFile) {
 836     AliError("Track-by-track correction file can not be opened!");
 837     return -1;
 838   }
 839
 840   // -- correction - not cross section corrected
 841   fCorr0    = new THnSparseF**[2];
 842   fCorr0[0] = new THnSparseF*[fCentralityBinMax];
 843   fCorr0[1] = new THnSparseF*[fCentralityBinMax];
 844
 845   for (Int_t idxCent = 0; idxCent < fCentralityBinMax; ++idxCent) {
 846
 847     // In future --> !!!!!!!!!!!!!!!!!!!!!!!
 848     // THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr0_Cent_%d", fPartName[0].Data(), idxCent)));
 849     // THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr0_Cent_%d", fPartName[1].Data(), idxCent)));
 850     THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("pbar_Corr0_Cent_%d", idxCent)));
 851     THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("p_Corr0_Cent_%d", idxCent)));
 852
 853     if (!sp0 || !sp1) {
 854       AliError(Form("Effective correction objects 0 - idxCent %d can not be retrieved!",idxCent));
 855       return -1;
 856     }
 857
 858     fCorr0[0][idxCent] = static_cast<THnSparseF*>(sp0->Clone());
 859     fCorr0[1][idxCent] = static_cast<THnSparseF*>(sp1->Clone());
 860   }
 861
 862   // -- correction - ross section corrected
 863   fCorr1    = new THnSparseF**[2];
 864   fCorr1[0] = new THnSparseF*[fCentralityBinMax];
 865   fCorr1[1] = new THnSparseF*[fCentralityBinMax];
 866
 867   for (Int_t idxCent = 0; idxCent < fCentralityBinMax; ++idxCent) {
 868
 869     // In future --> !!!!!!!!!!!!!!!!!!!!!!!
 870     //THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr1_Cent_%d", fPartName[0].Data(), idxCent)));
 871     //THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("%s_Corr1_Cent_%d", fPartName[1].Data(), idxCent)));
 872     THnSparseF *sp0 = static_cast<THnSparseF*>(corrFile->Get(Form("pbar_Corr1_Cent_%d", idxCent)));
 873     THnSparseF *sp1 = static_cast<THnSparseF*>(corrFile->Get(Form("p_Corr1_Cent_%d", idxCent)));
 874
 875     if (!sp0 || !sp1) {
 876       AliError(Form("Effective correction objects 1 - idxCent %d can not be retrieved!",idxCent));
 877       return -1;
 878     }
 879
 880     fCorr1[0][idxCent] = static_cast<THnSparseF*>(sp0->Clone());
 881     fCorr1[1][idxCent] = static_cast<THnSparseF*>(sp1->Clone());
 882   }
 883
 884   return 0;
 885 }
 886
 887 /*
 888  * ---------------------------------------------------------------------------------
 889  *                         Event / Trigger Statistics - private
 890  * ---------------------------------------------------------------------------------
 891  */
 892
 893 //________________________________________________________________________
 894 Bool_t AliAnalysisNetParticleHelper::FillEventStats(Int_t *aEventCuts) {
 895   // -- Fill event / centrality statistics
 896
 897   Bool_t isRejected = kFALSE;
 898
 899   // -- Fill event statistics
 900   for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
 901
 902     if (aEventCuts[idx])
 903       isRejected = kTRUE;
 904     else
 905       fHEventStat0->Fill(idx);
 906   }
 907
 908   for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
 909     if (aEventCuts[idx])
 910       break;
 911     fHEventStat1->Fill(idx);
 912   }
 913
 914   // -- Fill centrality statistics of accepted events
 915   if (!isRejected)
 916     fHCentralityStat->Fill(fCentralityBin);
 917
 918   return isRejected;
 919 }