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