a5e0d04c5065262b805d1b5f44b5b3fc8ea50ed2
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisNetParticleDistribution.cxx
1 //-*- Mode: C++ -*-
2
3 #include "TMath.h"
4 #include "TAxis.h"
5 #include "TSystem.h" 
6 #include "TProfile.h" 
7 #include "TH2F.h" 
8 #include "TH3F.h" 
9 #include "TFile.h" 
10 #include "TPRegexp.h"
11
12 #include "AliStack.h"
13 #include "AliMCEvent.h"
14 #include "AliMCParticle.h"
15 #include "AliESDtrackCuts.h"
16 #include "AliESDInputHandler.h"
17 #include "AliESDpid.h"
18 #include "AliCentrality.h"
19 #include "AliTracker.h"
20
21 #include "AliAnalysisNetParticleDistribution.h"
22
23 using namespace std;
24
25 // Task for NetParticle checks
26 // Author: Jochen Thaeder <jochen@thaeder.de>
27
28 ClassImp(AliAnalysisNetParticleDistribution)
29
30 /*
31  * ---------------------------------------------------------------------------------
32  *                            Constructor / Destructor
33  * ---------------------------------------------------------------------------------
34  */
35
36 //________________________________________________________________________
37 AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() :
38   fHelper(NULL),
39
40   fOutList(NULL),
41
42   fESDHandler(NULL),
43   fPIDResponse(NULL),
44   fESD(NULL),
45   fIsMC(kFALSE),
46   fMCEvent(NULL),
47   fStack(NULL),
48   fESDTrackCuts(NULL),
49   fEtaMax(0.9),
50   fPtRange(NULL),
51   fNp(NULL),
52   fNCorrNp(2),
53   fCorrNp(NULL),
54   fNMCNp(5),
55   fMCNp(NULL),
56   fNControlMCNp(5),
57   fControlMCNp(NULL),
58
59   fHnTrackUnCorr(NULL) {
60   // Constructor   
61   
62   AliLog::SetClassDebugLevel("AliAnalysisNetParticleDistribution",10);
63 }
64
65 //________________________________________________________________________
66 AliAnalysisNetParticleDistribution::~AliAnalysisNetParticleDistribution() {
67   // Destructor
68
69   if (fNp) delete[] fNp;  
70
71   for (Int_t ii = 0; ii < fNCorrNp; ++ii) 
72     if (fCorrNp[ii]) delete[] fCorrNp[ii];
73   if (fCorrNp) delete[] fCorrNp;
74
75   for (Int_t ii = 0; ii < fNMCNp; ++ii) 
76     if (fMCNp[ii]) delete[] fMCNp[ii];
77   if (fMCNp) delete[] fMCNp;
78
79   for (Int_t ii = 0; ii < fNControlMCNp; ++ii) 
80     if (fControlMCNp[ii]) delete[] fControlMCNp[ii];
81   if (fControlMCNp) delete[] fControlMCNp;
82
83   return;
84 }
85
86 /*
87  * ---------------------------------------------------------------------------------
88  *                                 Public Methods
89  * ---------------------------------------------------------------------------------
90  */
91
92 //________________________________________________________________________
93 Int_t AliAnalysisNetParticleDistribution::Initialize(AliAnalysisNetParticleHelper* helper, AliESDtrackCuts* cuts, Bool_t isMC, Float_t *ptRange, Float_t etaMax) {
94   // -- Initialize
95   
96   fHelper = helper;
97   fESDTrackCuts = cuts;
98   fIsMC = isMC;
99   fPtRange = ptRange;
100   fEtaMax = etaMax;
101
102   // ------------------------------------------------------------------
103   // -- N particles / N anti-particles
104   // ------------------------------------------------------------------
105   //  Np          : arr[particle]
106   //  CorrNp      : arr[corrSet][particle]
107   //  MCNp        : arr[corrSet][particle] - MC
108   //  ControlMCNp : arr[corrSet][particle] - Control MC
109
110   fNp     = new Float_t[2];
111
112   fCorrNp = new Float_t*[fNCorrNp];
113   for (Int_t ii = 0 ; ii < fNCorrNp; ++ii)
114     fCorrNp[ii] = new Float_t[2];
115
116   fMCNp = new Float_t*[fNMCNp];
117   for (Int_t ii = 0 ; ii < fNMCNp; ++ii)
118     fMCNp[ii] = new Float_t[2];
119
120   fControlMCNp = new Float_t*[fNControlMCNp];
121   for (Int_t ii = 0 ; ii < fNControlMCNp; ++ii)
122     fControlMCNp[ii] = new Float_t[2];
123
124   ResetEvent();
125
126   return 0;
127 }
128
129 //________________________________________________________________________
130 void AliAnalysisNetParticleDistribution::CreateHistograms(TList* outList) {
131   // -- Add histograms to outlist
132
133   fOutList = outList;
134   
135   // ------------------------------------------------------------------
136   // -- Create net particle histograms
137   // ------------------------------------------------------------------
138
139   AddHistSet("fHDist", "Uncorrected");
140   AddHistSet("fHDistCorr0", "Corrected [without cross section correction]");
141   AddHistSet("fHDistCorr1", "Corrected [with cross section correction]");
142   
143   if (fIsMC) {
144     AddHistSet("fHMCrapidity", "MC primary in |y| < 0.5");
145     AddHistSet("fHMCptMin",    "MC primary in |y| + #it{p}_{T} > 0.1");
146     AddHistSet("fHMCpt",       Form("MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
147     AddHistSet("fHMCeta",      Form("MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
148     AddHistSet("fHMCetapt",    Form("MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
149     
150     AddHistSet("fHControlMCLambdarapidity", "Control MC Lambda primary in |y| < 0.5");
151     AddHistSet("fHControlMCLambdaptMin",    "Control MC Lambda primary in |y| + #it{p}_{T} > 0.1");
152     AddHistSet("fHControlMCLambdapt",       Form("Control MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
153     AddHistSet("fHControlMCLambdaeta",      Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
154     AddHistSet("fHControlMCLambdaetapt",    Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
155   }
156
157   // ------------------------------------------------------------------
158   // -- Get Probe Particle Container
159   // ------------------------------------------------------------------
160
161   Double_t dCent[2] = {-0.5, 8.5};
162   Int_t iCent       = 9;
163   
164   Double_t dEta[2]  = {-0.9, 0.9}; // -> 37 bins
165   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; 
166
167   Double_t dRap[2]  = {-0.5, 0.5}; 
168   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; 
169
170   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
171   Int_t iPhi        = 42;
172
173   Double_t dPt[2]   = {0.1, 3.0}; 
174   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
175
176   Double_t dSign[2] = {-1.5, 1.5};
177   Int_t iSign       = 3;
178
179   //                          cent:     pt:     sign:     eta:     phi:       y
180   Int_t    binShort[6] = {   iCent,    iPt,    iSign,    iEta,    iPhi,    iRap};
181   Double_t minShort[6] = {dCent[0], dPt[0], dSign[0], dEta[0], dPhi[0], dRap[0]};
182   Double_t maxShort[6] = {dCent[1], dPt[1], dSign[1], dEta[1], dPhi[1], dRap[1]};
183
184   // -- UnCorrected
185   fOutList->Add(new THnSparseF("fHnTrackUnCorr", "cent:pt:sign:eta:phi:y", 6, binShort, minShort, maxShort));  
186   fHnTrackUnCorr = static_cast<THnSparseF*>(fOutList->Last());
187   fHnTrackUnCorr->GetAxis(0)->SetTitle("centrality");
188   fHnTrackUnCorr->GetAxis(1)->SetTitle("#it{p}_{T} (GeV/#it{c})");
189   fHnTrackUnCorr->GetAxis(2)->SetTitle("sign");
190   fHnTrackUnCorr->GetAxis(3)->SetTitle("#eta");
191   fHnTrackUnCorr->GetAxis(4)->SetTitle("#varphi");
192   fHnTrackUnCorr->GetAxis(5)->SetTitle("#it{y}");
193
194   fHelper->BinLogAxis(fHnTrackUnCorr, 1);
195
196   // ------------------------------------------------------------------
197
198   return;
199 }
200
201 //________________________________________________________________________
202 Int_t AliAnalysisNetParticleDistribution::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
203   // -- Setup Event
204
205   ResetEvent();
206
207   // -- Get ESD objects
208   fESDHandler  = esdHandler;
209   fPIDResponse = esdHandler->GetPIDResponse();
210   fESD         = fESDHandler->GetEvent();
211
212   // -- Get MC objects
213   fMCEvent     = mcEvent;
214   if (fMCEvent)
215     fStack     = fMCEvent->Stack();
216
217   return 0;
218 }
219
220 //________________________________________________________________________
221 void AliAnalysisNetParticleDistribution::ResetEvent() {
222   // -- Reset event
223   
224   // -- Reset ESD Event
225   fESD       = NULL;
226
227   // -- Reset MC Event
228   if (fIsMC)
229     fMCEvent = NULL;
230
231   // -- Reset N particles/anti-particles
232   for (Int_t jj = 0; jj < 2; ++jj)
233     fNp[jj] = 0.;
234   
235   // -- Reset N corrected particles/anti-particles
236   for (Int_t ii = 0; ii < fNCorrNp; ++ii) 
237     for (Int_t jj = 0; jj < 2; ++jj)
238       fCorrNp[ii][jj] = 0.;
239
240   // -- Reset N MC particles/anti-particles
241   for (Int_t ii = 0; ii < fNMCNp; ++ii) 
242     for (Int_t jj = 0; jj < 2; ++jj)
243       fMCNp[ii][jj] = 0.;
244
245   // -- Reset N control MC particles/anti-particles
246   for (Int_t ii = 0; ii < fNControlMCNp; ++ii) 
247     for (Int_t jj = 0; jj < 2; ++jj)
248       fControlMCNp[ii][jj] = 0.;
249 }
250
251 //________________________________________________________________________
252 Int_t AliAnalysisNetParticleDistribution::Process() {
253   // -- Process NetParticle Distributions
254   
255   // -- Fill ESD tracks
256   ProcessESDTracks();
257     
258   // -- Fill MC truth particles
259   if (fIsMC)  {
260     ProcessStackParticles();
261     ProcessStackControlParticles();
262   }
263
264   return 0;
265 }
266
267 //________________________________________________________________________
268 Int_t AliAnalysisNetParticleDistribution::ProcessESDTracks() {
269   // -- Process ESD tracks and fill histograms
270
271   for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
272     AliESDtrack *track = fESD->GetTrack(idxTrack); 
273
274     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
275     // -- Check track cuts
276     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
277     
278     // -- Check if charged track is accepted for basic parameters
279     if (!fHelper->IsTrackAcceptedBasicCharged(track))
280       continue;
281     
282     // -- Check if accepted
283     if (!fESDTrackCuts->AcceptTrack(track)) 
284       continue;
285
286     // -- Check if accepted in rapidity window
287     Double_t yP;
288     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
289       continue;
290
291     // -- Check if accepted bt PID from TPC or TPC+TOF
292     Double_t pid[2];
293     if (!fHelper->IsTrackAcceptedPID(track, pid))
294       continue;
295
296     // -- Check if accepted with thighter DCA cuts
297     if (fHelper->IsTrackAcceptedDCA(track))
298       continue;
299     
300     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
301     // -- Fill Probe Particle
302     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
303
304     Double_t aTrack[6] = {
305       Double_t(fHelper->GetCentralityBin()),       //  0 centrality 
306       track->Pt(),                    //  1 pt
307       track->GetSign(),               //  2 sign
308       track->Eta(),                   //  3 eta
309       track->Phi(),                   //  4 phi
310       yP                              //  5 rapidity
311     };
312     
313     fHnTrackUnCorr->Fill(aTrack);
314     
315     // -- Count particle / anti-particle 
316     // ------------------------------------------------------------------
317     //  idxPart = 0 -> anti particle
318     //  idxPart = 1 -> particle
319
320     Int_t idxPart = (track->GetSign() < 0) ? 0 : 1;
321     fNp[idxPart] += 1.;
322     
323     for (Int_t ii = 0; ii < fNCorrNp; ++ii) 
324       fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii);      
325     
326   } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
327
328   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
329   // -- Fill Particle Fluctuation Histograms
330   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
331
332   // -- Uncorrected
333   FillHistSet("fHDist", fNp);
334     
335   // -- Corrected 
336   for (Int_t ii = 0 ; ii < fNCorrNp; ++ii) 
337     FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]); 
338
339   return 0;
340 }
341
342 //________________________________________________________________________
343 Int_t AliAnalysisNetParticleDistribution::ProcessStackParticles() {
344   // -- Process primary particles from the stack and fill histograms
345
346   Int_t pdgCode    = AliPID::ParticleCode(fHelper->GetParticleSpecies());
347   
348   for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
349     TParticle* particle = fStack->Particle(idxMC);
350     if (!particle) 
351       continue;
352
353     // -- Check basic MC properties -> charged physical primary
354     if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
355       continue;
356     
357     // -- Check if particle / anti-particle
358     if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
359       continue;
360     
361     // -- Get particle : 0 anti-particle / 1 particle
362     Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
363
364     // >> NOW only anti-particle / particle 
365     // >> With idxPart
366     
367     // -- Check rapidity window
368     Double_t yMC;
369     if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
370       continue;
371     
372     fMCNp[0][idxPart] += 1.;         // -> MCrapidity
373     
374     // -- Check acceptance
375     if (particle->Pt() > 0.1 )
376       fMCNp[1][idxPart] += 1.;       // -> MCptMin
377     
378     if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
379       fMCNp[2][idxPart] += 1.;       // -> MCpt
380     
381     if (TMath::Abs(particle->Eta()) <= fEtaMax)
382       fMCNp[3][idxPart] += 1.;       // -> MCeta
383     
384     if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
385       fMCNp[4][idxPart] += 1.;       // -> MCetapt
386     
387   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
388   
389   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
390   // -- Fill Particle Fluctuation Histograms
391   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
392
393   FillHistSet("fHMCrapidity", fMCNp[0]);
394   FillHistSet("fHMCptMin",    fMCNp[1]);
395   FillHistSet("fHMCpt",       fMCNp[2]);
396   FillHistSet("fHMCeta",      fMCNp[3]);
397   FillHistSet("fHMCetapt",    fMCNp[4]);
398
399   return 0;
400 }
401
402 //________________________________________________________________________
403 Int_t AliAnalysisNetParticleDistribution::ProcessStackControlParticles() {
404   // -- Process primary particles from the stack and fill histograms
405
406   Int_t pdgCode    = fHelper->GetControlParticleCode();
407   Bool_t isNeutral  = fHelper->IsControlParticleNeutral();
408   const Char_t* name = fHelper->GetControlParticleName().Data();
409
410   for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
411     TParticle* particle = fStack->Particle(idxMC);
412     if (!particle) 
413       continue;
414     
415     // -- Check basic MC properties -> neutral or charged physical primary
416     if (isNeutral) {
417       if (!fHelper->IsParticleAcceptedBasicNeutral(particle, idxMC))
418         continue;
419     }
420     else {
421       if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
422         continue;
423     }
424     
425     // -- Check if particle / anti-particle
426     if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
427       continue;
428
429     // -- Get particle : 0 anti-particle / 1 particle
430     Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
431
432     // >> NOW only anti-particle / particle 
433     // >> With idxPart
434     
435     // -- Check rapidity window
436     Double_t yMC;
437     if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
438       continue;
439
440     fControlMCNp[0][idxPart] += 1.;         // -> ControlMCrapidity
441
442     // -- Check acceptance
443     if (particle->Pt() > 0.1 )
444       fControlMCNp[1][idxPart] += 1.;       // -> ControlMCptMin
445     
446     if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
447       fControlMCNp[2][idxPart] += 1.;       // -> ControlMCpt
448     
449     if (TMath::Abs(particle->Eta()) <= fEtaMax)
450       fControlMCNp[3][idxPart] += 1.;       // -> ControlMCeta
451     
452     if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
453       fControlMCNp[4][idxPart] += 1.;       // -> ControlMCetapt
454   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
455
456   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
457   // -- Fill Particle Fluctuation Histograms
458   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
459
460   FillHistSet(Form("fHControlMC%srapidity", name), fControlMCNp[0], 0);
461   FillHistSet(Form("fHControlMC%sptMin",    name), fControlMCNp[1], 1);
462   FillHistSet(Form("fHControlMC%spt",       name), fControlMCNp[2], 2);
463   FillHistSet(Form("fHControlMC%seta",      name), fControlMCNp[3], 3);
464   FillHistSet(Form("fHControlMC%setapt",    name), fControlMCNp[4], 4);
465
466   return 0;
467 }
468
469 ////////////////////////////////////////////////////////////////////////////////////////////
470
471
472 //________________________________________________________________________
473 void  AliAnalysisNetParticleDistribution::AddHistSet(const Char_t *name, const Char_t *title)  {
474   // -- Add histogram sets for particle and anti-particle
475   
476   const Char_t* aPartNames[]             = {"pbar",          "p"};
477   const Char_t* aPartTitles[]            = {"Anti-Proton",   "Proton"};
478   const Char_t* aPartAxisTitles[]        = {"#bar{p}",       "p"};
479
480   const Char_t* aControlPartNames[]      = {"lambdabar",     "lambda"};
481   const Char_t* aControlPartTitles[]     = {"Anti-Lambda",   "Lambda"};
482
483   fOutList->Add(new TList);
484   TList *list = static_cast<TList*>(fOutList->Last());
485   list->SetName(name) ;
486   list->SetOwner(kTRUE);
487
488   TString sName(name);
489
490   TString sPtTitle("");
491   if (!sName.Contains("fHMC") || !sName.Contains("fHControlMC"))
492     sPtTitle += Form("#it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]);
493   
494   for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
495     list->Add(new TH2F(Form("%s%s", name, aPartNames[idxPart]), 
496                        Form("%s %s Dist %s;Centrality;N_{%s}", title, aPartTitles[idxPart], sPtTitle.Data(), aPartAxisTitles[idxPart]), 
497                        24, -0.5, 11.5, 501, -0.5, 500.49));
498   } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
499    
500   list->Add(new TH2F(Form("%sNet%s", name, aPartNames[1]), 
501                      Form("%s Net %s Dist %s;Centrality;N_{p} - N_{#bar{p}}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5, 201, -100.5, 100.49));
502
503   list->Add(new TProfile(Form("%sNet%sM", name, aPartNames[1]), 
504                          Form("%s Net %s Dist %s;Centrality;N_{p} - N_{#bar{p}}", title, aPartTitles[1], sPtTitle.Data()),  24, -0.5, 11.5));
505   list->Add(new TProfile(Form("%sNet%s2M", name, aPartNames[1]), 
506                          Form("%s (Net %s)^{2} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{2}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
507   list->Add(new TProfile(Form("%sNet%s3M", name, aPartNames[1]), 
508                          Form("%s (Net %s)^{3} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{3}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
509   list->Add(new TProfile(Form("%sNet%s4M", name, aPartNames[1]), 
510                          Form("%s (Net %s)^{4} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{4}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
511   list->Add(new TProfile(Form("%sNet%s5M", name, aPartNames[1]), 
512                          Form("%s (Net %s)^{5} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{5}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
513   list->Add(new TProfile(Form("%sNet%s6M", name, aPartNames[1]), 
514                          Form("%s (Net %s)^{6} Dist %s;Centrality;(N_{p} - N_{#bar{p}})^{6}", title, aPartTitles[1], sPtTitle.Data()), 24, -0.5, 11.5));
515                          
516   list->Add(new TH2F(Form("%sNet%sOverSum", name, aPartNames[1]), 
517                      Form("%s (Net %s)/ Sum Dist %s;Centrality;(N_{p} - N_{#bar{p}})/(N_{p} + N_{#bar{p}})", title, aPartTitles[1], sPtTitle.Data()), 
518                      24, -0.5, 11.5, 801, -50.5, 50.49));
519
520   if (sName.Contains("fHControlMC")) {
521     for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
522       list->Add(new TH2F(Form("%s%sOver%s", name, aPartNames[idxPart], aControlPartNames[idxPart]), 
523                          Form("%s %s / %s Dist %s;Centrality;N_{%s}/N_{Tracks}", 
524                               title, aPartTitles[idxPart], aControlPartTitles[idxPart], sPtTitle.Data(), aPartAxisTitles[idxPart]), 
525                          24, -0.5, 11.5, 101, 0., 1.));
526     } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
527   }
528  
529   //  if (sName.Contains("fHDist")) {
530   //    list->Add(new TH3F(Form("%sNet%s_RecVsMC", name, aPartNames[1]), 
531   //                   Form("%s Net %s Dist - Rec Vs MC - %s;Centrality;(N_{p} - N_{#bar{p}})_{rec};(N_{p} - N_{#bar{p}})_{MC}", 
532   //                        title, aPartTitles[1], sPtTitle.Data()), 
533   //                   24, -0.5, 11.5, 201, -100.5, 100.49, 201, -100.5, 100.49));
534   //  }
535   
536   return;
537 }
538
539 //________________________________________________________________________
540 void AliAnalysisNetParticleDistribution::FillHistSet(const Char_t *name, Float_t *np, Int_t controlIdx)  {
541   // -- Add histogram sets for particle and anti-particle
542
543   TList *list = static_cast<TList*>(fOutList->FindObject(name));
544
545   Float_t centralityBin = fHelper->GetCentralityBin();
546
547   (static_cast<TH2F*>(list->FindObject(Form("%spbar", name))))->Fill(centralityBin, np[0]);
548   (static_cast<TH2F*>(list->FindObject(Form("%sp",    name))))->Fill(centralityBin, np[1]);
549
550   Float_t sumNp    = np[1]+np[0];
551   Float_t deltaNp  = np[1]-np[0];
552   Float_t deltaNp2 = deltaNp * deltaNp;
553   Float_t deltaNp3 = deltaNp2 * deltaNp;
554
555   (static_cast<TH2F*>(list->FindObject(Form("%sNetp",  name))))->Fill(centralityBin, deltaNp);
556
557   (static_cast<TProfile*>(list->FindObject(Form("%sNetpM",  name))))->Fill(centralityBin, deltaNp);
558   (static_cast<TProfile*>(list->FindObject(Form("%sNetp2M", name))))->Fill(centralityBin, deltaNp2);
559   (static_cast<TProfile*>(list->FindObject(Form("%sNetp3M", name))))->Fill(centralityBin, deltaNp2*deltaNp);
560   (static_cast<TProfile*>(list->FindObject(Form("%sNetp4M", name))))->Fill(centralityBin, deltaNp2*deltaNp2);
561   (static_cast<TProfile*>(list->FindObject(Form("%sNetp5M", name))))->Fill(centralityBin, deltaNp3*deltaNp2);
562   (static_cast<TProfile*>(list->FindObject(Form("%sNetp6M", name))))->Fill(centralityBin, deltaNp3*deltaNp3);
563
564   (static_cast<TH2F*>(list->FindObject(Form("%sNetpOverSum", name))))->Fill(centralityBin, deltaNp/sumNp);
565
566   TString sName(name);
567   if (sName.Contains("fHControlMC") && controlIdx >= 0) {
568     (static_cast<TH2F*>(list->FindObject(Form("%spbarOverlambdabar", name))))->Fill(centralityBin, np[0]/fControlMCNp[controlIdx][0]);
569     (static_cast<TH2F*>(list->FindObject(Form("%spOverlambda",       name))))->Fill(centralityBin, np[1]/fControlMCNp[controlIdx][1]);
570   }
571
572   //  if (sName.Contains("fHDist")) {
573   //    Float_t deltaMCNp  = fMCNp[controlIdx][1]-fMCNp[controldx][0];
574   //    (static_cast<TH3F*>(list->FindObject(Form("%sNetp_RecVsMC", name))))->Fill(centralityBin, deltaNp, deltaMCNp);
575   //  }
576
577   return;
578 }
579