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