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