Update PR task: drathee
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidRatioHelper.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: ALICE Offline.                                                 *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //=========================================================================//
17 //             AliEbyE Analysis for Particle Ratio Fluctuation             //
18 //                   Deepika Rathee  | Satyajit Jena                       //
19 //                   drathee@cern.ch | sjena@cern.ch                       //
20 //                  Date: Wed Jul  9 18:38:30 CEST 2014                    // 
21 //          New approch to find particle ratio to reduce memory            //
22 //                             (Test Only)                                 //
23 //=========================================================================//
24
25 #include "TMath.h"
26 #include "TAxis.h"
27 #include "TSystem.h" 
28 #include "TFile.h" 
29 #include "TPRegexp.h"
30
31 #include "AliStack.h"
32 #include "AliMCEvent.h"
33 #include "AliMCParticle.h"
34 #include "AliESDtrackCuts.h"
35 #include "AliESDInputHandler.h"
36 #include "AliESDpid.h"
37 #include "AliAODInputHandler.h"
38 #include "AliAODEvent.h"
39 #include "AliAODMCParticle.h"
40 #include "AliCentrality.h"
41 #include "AliTracker.h"
42
43 #include "AliEbyEPidRatioHelper.h"
44
45 using namespace std;
46
47
48 ClassImp(AliEbyEPidRatioHelper)
49 //________________________________________________________________________
50 AliEbyEPidRatioHelper::AliEbyEPidRatioHelper() :
51   fModeDistCreation(0),
52
53   fInputEventHandler(NULL),
54   fPIDResponse(NULL),
55   fESD(NULL),
56   fESDTrackCuts(NULL),
57   fAOD(NULL),
58   fAODtrackCutBit(1024),
59   fIsMC(kFALSE),
60   fMCEvent(NULL),
61   fStack(NULL),
62
63   fCentralityBin(-1),
64   fCentralityPercentile(-1.),
65
66   fCentralityBinMax(11),
67   fVertexZMax(10.),
68   fRapidityMax(0.9),
69   fPhiMin(0.),
70   fPhiMax(TMath::TwoPi()),
71   fMinTrackLengthMC(70.),
72   fNSigmaMaxCdd(3.),
73   fNSigmaMaxCzz(3.),
74
75   fPIDStrategy(0),
76   fNSigmaMaxITS(4.),
77   fNSigmaMaxTPC(4.),
78   fNSigmaMaxTPClow(3.),
79   fNSigmaMaxTOF(4.),
80   fMinPtForTOFRequired(0.69),
81   fMaxPtForTPClow(0.69),
82
83   fHEventStat0(NULL),
84   fHEventStat1(NULL),
85   fHEventStatMax(6),
86
87   fHTriggerStat(NULL),
88   fNTriggers(5),
89
90   fHCentralityStat(NULL),
91   fHCentralityPer(NULL),
92   fHCentralityPerAll(NULL),
93   fNCentralityBins(11),
94   
95   fRandom(NULL),
96   fIsRatio(kFALSE), 
97   fIsPtBin(kFALSE) {
98   // Constructor   
99   
100   AliLog::SetClassDebugLevel("AliEbyEPidRatioHelper",10);
101 }
102
103 const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthRap = 0.1;
104 const Float_t AliEbyEPidRatioHelper::fgkfHistBinWitdthPt  = 0.3; // 0.08 // 300 MeV  // was 80 MeV
105
106 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeCent[]  = {-0.5, 10.5};
107 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsCent    = 11 ;
108
109 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeEta[]   = {-0.9, 0.9};
110 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsEta     = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeEta[1] -
111                                                                    AliEbyEPidRatioHelper::fgkfHistRangeEta[0]) / 
112                                                                   AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
113
114 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeRap[]   = {-0.8, 0.8};
115 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsRap     = Int_t((AliEbyEPidRatioHelper::fgkfHistRangeRap[1] - AliEbyEPidRatioHelper::fgkfHistRangeRap[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthRap) +1;
116
117 const Float_t AliEbyEPidRatioHelper::fgkfHistRangePhi[]   = {0.0, TMath::TwoPi()};
118 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsPhi     = 21 ;
119
120 const Float_t AliEbyEPidRatioHelper::fgkfHistRangePt[]    = {0.2, 2.9}; // {0.2, 5.}; // was {0.3, 2.22}
121 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsPt      = Int_t((AliEbyEPidRatioHelper::fgkfHistRangePt[1] - AliEbyEPidRatioHelper::fgkfHistRangePt[0]) / AliEbyEPidRatioHelper::fgkfHistBinWitdthPt); 
122
123 const Float_t AliEbyEPidRatioHelper::fgkfHistRangeSign[]  = {-1.5, 1.5};
124 const Int_t   AliEbyEPidRatioHelper::fgkfHistNBinsSign    =  3;
125
126 const Char_t* AliEbyEPidRatioHelper::fgkEventNames[]         = {"All", "IsTriggered", "HasVertex", "Vz<Vz_{Max}", "Centrality [0,100]%"};
127 const Char_t* AliEbyEPidRatioHelper::fgkCentralityMaxNames[] = {"5", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"};
128 const Char_t* AliEbyEPidRatioHelper::fgkTriggerNames[]       = {"kMB", "kCentral", "kSemiCentral", "kEMCEJE", "kEMCEGA" }; 
129 const Char_t* AliEbyEPidRatioHelper::fgkCentralityNames[]    = {"0-5%", "5-10%", "10-20%", "20-30%", "30-40%", "40-50%","50-60%", "60-70%", "70-80%", "80-90%", "90-100%"};
130
131 const Char_t* AliEbyEPidRatioHelper::fgkPidName[4]      = {"Nch","Npi","Nka","Npr"};
132 const Char_t* AliEbyEPidRatioHelper::fgkPidLatex[4][2]  = {{"N_{-}","N_{+}"}, {"N_{#pi^{-}}","N_{#pi^{+}}"},{"N_{K^{-}}","N_{K^{+}}"}, {"N_{#bar{p}}","N_{p}"}};
133 const Char_t* AliEbyEPidRatioHelper::fgkPidTitles[4][2] = {{"Negative","Positive"},{"Anti-Pions","Pions"},{"Anti-Kaons","Kaons"}, {"Anti-Protons","Protons"}};
134
135
136
137 //________________________________________________________________________
138 AliEbyEPidRatioHelper::~AliEbyEPidRatioHelper() {
139   // Destructor
140
141   if (fRandom)
142     delete fRandom;
143
144   return;
145 }
146
147 //________________________________________________________________________
148 void AliEbyEPidRatioHelper::SetPhiRange(Float_t f1, Float_t f2) {
149   // -- Set phi range and adopt to phi-histogram
150   
151   fPhiMin = f1; 
152   fPhiMax = (f1 < f2) ? f2 : f2+TMath::TwoPi();
153
154   Float_t phiMin = fPhiMin;
155   Float_t phiMax = fPhiMax;
156   
157   // -- Update Ranges
158   Float_t binWidth = (AliEbyEPidRatioHelper::fgkfHistRangePhi[1] - AliEbyEPidRatioHelper::fgkfHistRangePhi[0]) / 
159     Float_t(AliEbyEPidRatioHelper::fgkfHistNBinsPhi);
160
161   Float_t lowEdge  = AliEbyEPidRatioHelper::fgkfHistRangePhi[0] - binWidth;
162   Float_t highEdge = AliEbyEPidRatioHelper::fgkfHistRangePhi[0];
163
164   for (Int_t ii = 1; ii <= AliEbyEPidRatioHelper::fgkfHistNBinsPhi; ++ii) {
165     lowEdge += binWidth;
166     highEdge += binWidth;
167
168     if (phiMin >= lowEdge && phiMin < highEdge ) 
169       phiMin = lowEdge;
170     if (phiMax > lowEdge && phiMax <= highEdge ) 
171       phiMax = highEdge;
172   }
173   
174   printf(">>>> Update Phi Range : [%f,%f] -> [%f,%f]\n", fPhiMin, fPhiMax, phiMin, phiMax);
175   fPhiMin = phiMin;
176   fPhiMax = phiMax;
177 }
178
179
180
181 //________________________________________________________________________
182 Int_t AliEbyEPidRatioHelper::Initialize(AliESDtrackCuts *cuts, Bool_t isMC, Bool_t isRatio, Bool_t isPtBin, Int_t trackCutBit, Int_t modeDistCreation) {
183   // -- Initialize helper
184
185   Int_t iResult = 0;
186
187   // -- ESD track cuts
188   fESDTrackCuts     = cuts;
189
190   
191   // -- Is MC
192   fIsMC             = isMC;
193   fIsRatio          = isRatio;
194   fIsPtBin          = isPtBin;
195
196
197   
198
199   // -- AOD track filter bit
200   fAODtrackCutBit   = trackCutBit;
201     
202   // -- mode Distribution creation
203   fModeDistCreation = modeDistCreation;
204
205   // -- Setup event cut statistics 
206   InitializeEventStats();
207
208   // -- Setup trigger statistics 
209   InitializeTriggerStats();
210
211   // -- Setup centrality statistics 
212   InitializeCentralityStats();
213
214   // -- PRINT PID Strategy
215   //    0 :   TPC(TPClow+TPCHigh)
216   //    1 :   ITS
217   //    2 :   TOF
218   //    3 :   ITS+TPC(TPClow+TPCHigh)
219   //    4 :   TPC(TPClow+TPCHigh)+TOF
220   //    5 :   TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
221   //    6 :   TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
222   //    7 :   TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
223   //    8 :   TPC(TPClow+TPCHigh)+ITS+TOF 
224   printf(">>>>  PID STRATEGY: %d || sigmaMax: ITS %.2f TPC %.2f TOF %.2f \n", fPIDStrategy, fNSigmaMaxITS, fNSigmaMaxTPC, fNSigmaMaxTOF);            
225                       
226   // -- Initialize random number generator
227   fRandom = new TRandom3();
228   fRandom->SetSeed();
229                       
230   return iResult;
231 }
232
233 //________________________________________________________________________
234 Int_t AliEbyEPidRatioHelper::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler, AliMCEvent *mcEvent) {
235   // -- Setup Event
236   
237   // -- Get ESD objects
238   if(esdHandler){
239     fInputEventHandler = static_cast<AliInputEventHandler*>(esdHandler);
240     fESD               = dynamic_cast<AliESDEvent*>(fInputEventHandler->GetEvent());
241     if (!fESD) {
242       AliError("ESD event handler not available");
243       return -1;
244     }
245   }
246
247   // -- Get AOD objects
248   else if(aodHandler){
249     fInputEventHandler = static_cast<AliInputEventHandler*>(aodHandler);
250     fAOD               = dynamic_cast<AliAODEvent*>(fInputEventHandler->GetEvent());
251     if (!fAOD) {
252       AliError("AOD event handler not available");
253       return -1;
254     }
255   }
256
257   // -- Get Common objects
258   fPIDResponse = fInputEventHandler->GetPIDResponse();
259
260   // -- Get MC objects
261   fMCEvent     = mcEvent;
262   if (fMCEvent)
263     fStack     = fMCEvent->Stack();
264
265   // -- Get event centrality
266   // >  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90|90-100 --> 11 bins
267   // >   0   1     2     3     4     5     6     7     8     9     10
268
269   AliCentrality *centrality = NULL;
270
271   if(esdHandler)
272     centrality = fESD->GetCentrality();
273   else if(aodHandler)
274     centrality = fAOD->GetHeader()->GetCentralityP();
275
276   if (!centrality) {
277     AliError("Centrality not available");
278     return -1;
279   }
280
281   Int_t centBin = centrality->GetCentralityClass10("V0M");
282   if (centBin == 0) { fCentralityBin = centrality->GetCentralityClass5("V0M"); }
283   else if (centBin == 11 || centBin == -1.)           { fCentralityBin = -1; }
284   else if (centBin > 0 && centBin < fNCentralityBins) { fCentralityBin = centBin + 1; }
285   else {  fCentralityBin = -2; }
286   
287   if (fCentralityBin >= fCentralityBinMax)
288     fCentralityBin = -2;
289
290   fCentralityPercentile = centrality->GetCentralityPercentile("V0M");
291   
292   return 0;
293 }
294 //________________________________________________________________________
295 Bool_t AliEbyEPidRatioHelper::IsEventTriggered() {
296   // -- Check if Event is triggered and fill Trigger Histogram
297   
298   Bool_t *aTriggerFired = new Bool_t[fNTriggers];
299   for (Int_t ii = 0; ii < fNTriggers; ++ii)
300     aTriggerFired[ii] = kFALSE;
301
302   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kMB))          aTriggerFired[0] = kTRUE;
303   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kCentral))     aTriggerFired[1] = kTRUE;
304   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kSemiCentral)) aTriggerFired[2] = kTRUE;
305   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEJE))      aTriggerFired[3] = kTRUE;
306   if ((fInputEventHandler->IsEventSelected() & AliVEvent::kEMCEGA))      aTriggerFired[4] = kTRUE;
307
308   Bool_t isTriggered = kFALSE;
309
310   for (Int_t ii=0; ii<fNTriggers; ++ii) {
311     if(aTriggerFired[ii]) {
312       isTriggered = kTRUE;
313       fHTriggerStat->Fill(ii);
314     }
315   }
316   
317   delete[] aTriggerFired;
318
319   return isTriggered;
320 }
321
322 //________________________________________________________________________
323 Bool_t AliEbyEPidRatioHelper::IsEventRejected() {
324   // -- Evaluate event statistics histograms
325     
326   Int_t *aEventCuts = new Int_t[fHEventStatMax];
327   // set aEventCuts[ii] to 1 in case of reject
328   
329   for (Int_t ii=0;ii<fHEventStatMax; ++ii)
330     aEventCuts[ii] = 0;
331
332   Int_t iCut = 0;
333
334   // -- 0 - Before Physics Selection   
335   aEventCuts[iCut] = 0;
336
337   // -- 1 - No Trigger fired
338   ++iCut;
339   if (!IsEventTriggered())
340     aEventCuts[iCut] = 1;
341
342   // -- 2 - No Vertex 
343   ++iCut;
344   const AliESDVertex* vtxESD = NULL;
345   const AliAODVertex* vtxAOD = NULL;
346   if (fESD){
347     vtxESD = fESD->GetPrimaryVertexTracks();
348     if (!vtxESD)
349       aEventCuts[iCut] = 1;
350   }
351   else if (fAOD){
352     vtxAOD = fAOD->GetPrimaryVertex();
353     if (!vtxAOD)
354       aEventCuts[iCut] = 1;
355   }
356
357   // -- 3 - Vertex z outside cut window
358   ++iCut;
359   if (vtxESD){
360     if(TMath::Abs(vtxESD->GetZv()) > fVertexZMax) 
361       aEventCuts[iCut] = 1;
362   }
363   else if(vtxAOD){
364     if(TMath::Abs(vtxAOD->GetZ()) > fVertexZMax) 
365       aEventCuts[iCut] = 1;
366   }
367   else
368     aEventCuts[iCut] = 1;
369
370   // -- 4 - Centrality = -1  (no centrality or not hadronic)
371   ++iCut;
372   if(fCentralityBin == -1.) 
373     aEventCuts[iCut] = 1;
374
375   // -- 5 - Centrality < fCentralityMax
376   ++iCut;
377   if(fCentralityBin == -2.) 
378     aEventCuts[iCut] = 1;
379
380   // -- Fill statistics / reject event
381   Bool_t isRejected = FillEventStats(aEventCuts);
382
383   // -- Cleanup 
384   delete[] aEventCuts;
385
386   //cout << isRejected << endl;
387
388   return isRejected;
389 }
390
391 //________________________________________________________________________
392 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicCharged(AliVParticle *particle, Int_t idxMC) {
393   // -- Check if MC particle is accepted for basic parameters
394   
395   if (!particle) 
396     return kFALSE;
397
398   // -- check if charged
399   if (particle->Charge() == 0.0) 
400     return kFALSE;
401   
402   // -- check if physical primary - ESD
403   if (fESD) {
404     if(!fStack->IsPhysicalPrimary(idxMC)) 
405       return kFALSE;
406   }
407   // -- check if physical primary - AOD
408   else {
409     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary()) 
410       return kFALSE;
411   }
412   
413   return kTRUE;
414 }
415
416 //________________________________________________________________________
417 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedBasicNeutral(AliVParticle *particle, Int_t idxMC) {
418   // -- Check if MC particle is accepted for basic parameters
419   
420   if (!particle) 
421     return kFALSE;
422
423   // -- check if charged
424   if (particle->Charge() != 0.0) 
425     return kFALSE;
426   
427   // -- check if physical primary - ESD
428   if (fESD) {
429     if(!fStack->IsPhysicalPrimary(idxMC)) 
430       return kFALSE;
431   }
432   // -- check if physical primary - AOD
433   else {
434     if(!(static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary()) 
435       return kFALSE;
436   }
437   
438   return kTRUE;
439 }
440
441 //________________________________________________________________________
442 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedRapidity(AliVParticle *particle, Double_t &yP, Int_t gCurPid) {
443   
444  if (gCurPid == 0) {
445    yP = particle->Eta();
446    return kTRUE;
447  }
448   
449  Double_t mP = AliPID::ParticleMass(AliPID::kPion);
450  if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
451  else if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
452  else if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
453  
454  // -- Calculate rapidities and kinematics
455   Double_t p  = particle->P();
456   Double_t pz = particle->Pz();
457
458   Double_t eP = TMath::Sqrt(p*p + mP*mP);
459   yP          = 0.5 * TMath::Log((eP + pz) / (eP - pz));  
460
461   // -- Check Rapidity window
462   if (TMath::Abs(yP) > fRapidityMax)
463     return kFALSE;
464   
465   return kTRUE;
466 }
467
468 //________________________________________________________________________
469 Bool_t AliEbyEPidRatioHelper::IsParticleAcceptedPhi(AliVParticle *particle) {
470   // -- Check if particle is accepted
471   // > in phi
472   // > return 0 if not accepted
473   
474   if (particle->Phi() > fPhiMin && particle->Phi() <= fPhiMax)
475     return kTRUE;
476   else if (particle->Phi() < fPhiMin && (particle->Phi() + TMath::TwoPi()) <= fPhiMax)
477     return kTRUE;
478   else
479     return kFALSE;
480 }
481
482 //_____________________________________________________________________________
483 Bool_t AliEbyEPidRatioHelper::IsParticleFindable(Int_t label) {
484   // -- Check if MC particle is findable tracks
485
486   AliMCParticle *mcParticle = static_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
487   if(!mcParticle) 
488     return kFALSE;
489   
490   Int_t counter; 
491   Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(), 0.05, counter, 3.0); 
492
493   return (tpcTrackLength > fMinTrackLengthMC);    
494 }
495 //________________________________________________________________________
496 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedBasicCharged(AliVTrack* track) {
497   // -- Check if track is accepted 
498   // > for basic parameters
499
500   if (!track)
501     return kFALSE;
502   
503   if (track->Charge() == 0) 
504     return kFALSE;
505   
506   return kTRUE;
507
508  
509 //________________________________________________________________________
510 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedRapidity(AliVTrack *track, Double_t &yP, Int_t gCurPid) {
511    if (gCurPid == 0) {
512     yP = track->Eta();
513     return kTRUE;
514   }
515   
516   Double_t mP = AliPID::ParticleMass(AliPID::kPion);
517   if(gCurPid == 1) mP = AliPID::ParticleMass(AliPID::kPion);
518   else if(gCurPid == 2) mP = AliPID::ParticleMass(AliPID::kKaon);
519   else if(gCurPid == 3) mP = AliPID::ParticleMass(AliPID::kProton);
520
521   // -- Calculate rapidities and kinematics
522   Double_t pvec[3];
523   track->GetPxPyPz(pvec);
524
525   Double_t p  = track->P();
526   Double_t eP = TMath::Sqrt(p*p + mP*mP);
527            yP = 0.5 * TMath::Log((eP + pvec[2]) / (eP - pvec[2]));
528   
529   // -- Check Rapidity window
530   if (TMath::Abs(yP) > fRapidityMax)
531     return kFALSE;
532   
533   return kTRUE;
534 }
535
536 //________________________________________________________________________
537 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedDCA(AliVTrack *vTrack) {
538    Bool_t isAccepted = kTRUE;
539
540   if (!fESD)
541     return isAccepted;
542
543   AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
544   if (!track)
545     return kFALSE;
546   
547   // -- Get nHits SPD
548   if (track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) {
549
550     // -- Get DCA nSigmas
551     Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
552     track->GetImpactParameters(dca,cov);
553
554     Float_t nSigmaCdd = (cov[0] != 0.) ? dca[0]/TMath::Sqrt(cov[0]) : -9.99; 
555     Float_t nSigmaCzz = (cov[2] != 0.) ? dca[1]/TMath::Sqrt(cov[2]) : -9.99; 
556     
557     if (fNSigmaMaxCdd != 0.) {
558       if (TMath::Abs(nSigmaCdd) > fNSigmaMaxCdd)
559         isAccepted = kFALSE;
560     }
561
562     if (fNSigmaMaxCzz != 0.) {
563       if (TMath::Abs(nSigmaCzz) > fNSigmaMaxCzz)
564         isAccepted = kFALSE;
565     }
566   }
567
568   return isAccepted;
569 }
570
571 //________________________________________________________________________
572 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedPID(AliVTrack *track, Double_t* pid, AliPID::EParticleType gCurPid) {
573
574   Bool_t isAcceptedITS    = kFALSE;
575   Bool_t isAcceptedTPC    = kFALSE;
576   Bool_t isAcceptedTPClow = kFALSE;
577   Bool_t isAcceptedTOF    = kFALSE;
578   Bool_t isAccepted       = kFALSE;
579
580   // -- In case not PID is used
581   if (gCurPid == 0) {
582     pid[0] = 10.;
583     pid[1] = 10.;
584     pid[2] = 10.;
585     return kTRUE;
586   }
587   
588   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kITS, track, gCurPid, pid[0]) == AliPIDResponse::kDetPidOk) {
589     if (TMath::Abs(pid[0]) < fNSigmaMaxITS) 
590       isAcceptedITS = kTRUE;
591   }
592
593   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, gCurPid, pid[1]) == AliPIDResponse::kDetPidOk) {
594     if (TMath::Abs(pid[1]) < fNSigmaMaxTPC) 
595       isAcceptedTPC = kTRUE;
596     if (TMath::Abs(pid[1]) < fNSigmaMaxTPClow) 
597       isAcceptedTPClow = kTRUE;
598     if (track->Pt() < fMaxPtForTPClow)
599       isAcceptedTPC = isAcceptedTPClow;
600   }
601
602   Bool_t hasPIDTOF = kFALSE;
603   if (fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, gCurPid, pid[2]) == AliPIDResponse::kDetPidOk) {
604     hasPIDTOF = kTRUE;
605     if (TMath::Abs(pid[2]) < fNSigmaMaxTOF) 
606       isAcceptedTOF = kTRUE;
607   }
608
609   // -- Check TOF missmatch for MC
610   
611   //if (ESD)
612   if (fIsMC && isAcceptedTOF) {
613     Int_t tofLabel[3];                                                                                                                                        
614     //    AliESDtrack* track = dynamic_cast<AliESDtrack*>(vTrack);
615     // TODO add code for AOD 
616
617     (dynamic_cast<AliESDtrack*>(track))->GetTOFLabel(tofLabel);
618    
619     Bool_t hasMatchTOF = kTRUE;
620     if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) 
621       hasMatchTOF = kFALSE;
622
623     TParticle *matchedTrack = fStack->Particle(TMath::Abs(tofLabel[0]));
624     if (TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel())) 
625       hasMatchTOF = kTRUE;
626
627     isAcceptedTOF = hasMatchTOF;
628   }
629
630   //    0 :   TPC(TPClow+TPCHigh)
631   //    1 :   ITS
632   //    2 :   TOF
633   //    3 :   ITS+TPC(TPClow+TPCHigh)
634   //    4 :   TPC(TPClow+TPCHigh)+TOF
635   //    5 :   TPC(TPClow+TPCHigh)+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
636   //    6 :   TPC(TPClow+TPCHigh)+ITS+TOF with TOF only for those tracks which have TOF information
637   //    7 :   TPC(TPClow+TPCHigh)+ITS+TOF for pT >= fMinPtForTOFRequired TOF is required, below, only used if there
638   //    8 :   TPC(TPClow+TPCHigh)+ITS+TOF 
639   if (fPIDStrategy == 0) {             //  TPC PID
640     isAccepted = isAcceptedTPC;
641   }
642   else if (fPIDStrategy == 1) {        //  ITS PID
643     isAccepted = isAcceptedITS;
644   }
645   else if (fPIDStrategy == 2) {        //  TOF PID
646     isAccepted = isAcceptedTOF;
647   }
648   else if (fPIDStrategy == 3) {        //  TPC+ITS PID
649     isAccepted = isAcceptedTPC && isAcceptedITS;
650   }
651   else if (fPIDStrategy == 4) {        //  TPC+TOF PID
652     isAccepted = isAcceptedTPC && isAcceptedTOF;
653   }
654   else if (fPIDStrategy == 5) {        //  TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
655     if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired) 
656       isAcceptedTOF = kTRUE;
657     isAccepted = isAcceptedTPC && isAcceptedTOF;
658   }
659   else if (fPIDStrategy == 6) {        //  ITS+TPC+TOF PID -- TOF only for those tracks which have TOF information
660     isAccepted = isAcceptedTPC && isAcceptedITS;
661     if (hasPIDTOF)
662       isAccepted = isAccepted && isAcceptedTOF;
663   }
664   else if (fPIDStrategy == 7) {        //  ITS+TPC+TOF PID -- for pT >= fMinPtForTOFRequired TOF is required
665     if (!hasPIDTOF && track->Pt() < fMinPtForTOFRequired)
666       isAcceptedTOF = kTRUE;
667     isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
668   }
669   else if (fPIDStrategy == 8) {        //  ITS+TPC+TOF PID
670     isAccepted = isAcceptedITS && isAcceptedTPC && isAcceptedTOF;
671   }
672
673   return isAccepted;
674 }
675
676 //________________________________________________________________________
677 Bool_t AliEbyEPidRatioHelper::IsTrackAcceptedPhi(AliVTrack *track) {
678   // -- Check if track is accepted
679   // > in phi
680   // > return 0 if not accepted
681   
682   if (track->Phi() > fPhiMin && track->Phi() <= fPhiMax)
683     return kTRUE;
684   else if (track->Phi() < fPhiMin && (track->Phi() + TMath::TwoPi()) <= fPhiMax)
685     return kTRUE;
686   else
687     return kFALSE;
688 }
689
690 //________________________________________________________________________
691 void AliEbyEPidRatioHelper::BinLogAxis(const THnBase *hn, Int_t axisNumber, AliESDtrackCuts* cuts) {
692     AliESDtrackCuts* esdTrackCuts = (cuts) ? cuts : fESDTrackCuts;
693
694   // -- Make logarithmic binning 
695   TAxis *axis = hn->GetAxis(axisNumber);
696   Int_t  nBins = axis->GetNbins();
697
698   Double_t from  = axis->GetXmin();
699   Double_t to    = axis->GetXmax();
700   Double_t *newBins = new Double_t[nBins + 1];
701    
702   newBins[0] = from;
703   Double_t factor = TMath::Power(to/from, 1./nBins);
704   
705   for (int ii = 1; ii <= nBins; ii++)
706    newBins[ii] = factor * newBins[ii-1];
707   
708   axis->Set(nBins, newBins);
709
710   delete [] newBins;
711
712   // -- Update Ranges
713   // ------------------
714   Float_t ptRange[2];
715   Float_t oldPtRange[2];
716   esdTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
717   esdTrackCuts->GetPtRange(oldPtRange[0],oldPtRange[1]);
718
719   Float_t minPtForTOFRequired = fMinPtForTOFRequired;
720   Float_t maxPtForTPClow      = fMaxPtForTPClow;
721
722   // -- Update minPt Cut
723   Int_t bin = axis->FindBin(ptRange[0]+10e-7);
724   ptRange[0] = axis->GetBinLowEdge(bin); 
725
726   // -- Update maxPt Cut
727   bin = axis->FindBin(ptRange[1]-10e-7);
728   ptRange[1] = axis->GetBinUpEdge(bin); 
729
730   if (ptRange[0] != oldPtRange[0] || ptRange[1] != oldPtRange[1]) {
731     printf(">>>> Update Pt Range : [%f,%f] -> [%f,%f]\n", oldPtRange[0], oldPtRange[1], ptRange[0], ptRange[1]);
732     esdTrackCuts->SetPtRange(ptRange[0],ptRange[1]);
733   }
734
735   // -- Update MinPtForTOFRequired
736   bin = axis->FindBin(minPtForTOFRequired-10e-7);
737   minPtForTOFRequired = axis->GetBinUpEdge(bin); 
738
739   if (minPtForTOFRequired != fMinPtForTOFRequired) {
740     printf(">>>> Update Min Pt for TOF : %f -> %f\n", fMinPtForTOFRequired, minPtForTOFRequired);
741     fMinPtForTOFRequired = minPtForTOFRequired;
742   }
743
744   // -- Update MaxPtForTPClow
745   bin = axis->FindBin(maxPtForTPClow-10e-7);
746   maxPtForTPClow = axis->GetBinUpEdge(bin); 
747
748   if (maxPtForTPClow != fMaxPtForTPClow) {
749     printf(">>>> Update Max Pt for TPC Low : %f -> %f\n", fMaxPtForTPClow, maxPtForTPClow);
750     fMaxPtForTPClow = maxPtForTPClow;
751   }
752 }
753
754 //________________________________________________________________________
755 void AliEbyEPidRatioHelper::InitializeEventStats() {
756   // -- Initialize event statistics histograms
757
758   fHEventStat0 = new TH1F("hEventStat0","Event cut statistics 0;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
759   fHEventStat1 = new TH1F("hEventStat1","Event cut statistics 1;Event Cuts;Events", fHEventStatMax,-0.5,fHEventStatMax-0.5);
760
761   for ( Int_t ii=0; ii < fHEventStatMax-1; ii++ ) {
762     fHEventStat0->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkEventNames[ii]);
763     fHEventStat1->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkEventNames[ii]);
764   }
765
766   fHEventStat0->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
767   fHEventStat1->GetXaxis()->SetBinLabel(fHEventStatMax, Form("Centrality [0-%s]%%", AliEbyEPidRatioHelper::fgkCentralityMaxNames[fCentralityBinMax-1]));
768 }
769
770 //________________________________________________________________________
771 void AliEbyEPidRatioHelper::InitializeTriggerStats() {
772   // -- Initialize trigger statistics histograms
773
774   fHTriggerStat = new TH1F("hTriggerStat","Trigger statistics;Trigger;Events", fNTriggers,-0.5,fNTriggers-0.5);
775
776   for ( Int_t ii=0; ii < fNTriggers; ii++ )
777     fHTriggerStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkTriggerNames[ii]);
778 }
779
780 //________________________________________________________________________
781 void AliEbyEPidRatioHelper::InitializeCentralityStats() {
782   // -- Initialize trigger statistics histograms
783
784   fHCentralityStat = new TH1F("hCentralityStat","Centrality statistics;Centrality Bins;Events", 
785                               fNCentralityBins,-0.5,fNCentralityBins-0.5);
786   
787   for ( Int_t ii=0; ii < fNCentralityBins; ii++ )
788     fHCentralityStat->GetXaxis()->SetBinLabel(ii+1, AliEbyEPidRatioHelper::fgkCentralityNames[ii]);
789   
790   fHCentralityPer = new TH1F("hCentralityPercentileAccepted","Centrality Percentile statistics;Centrality Bins;Events", 
791                              100,-0.5,99.5);
792
793   fHCentralityPerAll = new TH1F("hCentralityPercentileAll","Centrality Percentile statistics;Centrality Bins;Events", 
794                              100,-0.5,99.5);
795
796
797 }
798 //________________________________________________________________________
799 Bool_t AliEbyEPidRatioHelper::FillEventStats(Int_t *aEventCuts) {
800   // -- Fill event / centrality statistics 
801
802   Bool_t isRejected = kFALSE;
803
804   // -- Fill event statistics
805   for (Int_t idx = 0; idx < fHEventStatMax ; ++idx) {
806
807     if (aEventCuts[idx])
808       isRejected = kTRUE;
809     else
810       fHEventStat0->Fill(idx);
811   }
812   
813   for (Int_t idx = 0; idx < fHEventStatMax; ++idx) {
814     if (aEventCuts[idx])
815       break;
816     fHEventStat1->Fill(idx);
817   }
818
819   // -- Fill centrality statistics of accepted events
820   if (!isRejected) {
821     fHCentralityStat->Fill(fCentralityBin);
822     fHCentralityPer->Fill(fCentralityPercentile);
823   }
824
825   fHCentralityPerAll->Fill(fCentralityPercentile);
826
827
828
829   return isRejected;
830 }