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