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