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