Migrating PWG2/SPECTRA/Fit to new PWG structure
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / multPbPb / AliAnalysisMultPbTrackHistoManager.cxx
1 #include "AliAnalysisMultPbTrackHistoManager.h"
2 #include "AliLog.h"
3 #include "TH1.h"
4 #include "TH3D.h"
5 #include "TH1I.h"
6 #include "TROOT.h"
7 #include "TMCProcess.h"
8 #include "AliMCParticle.h"
9
10 #include <iostream>
11 #include "TH2D.h"
12
13 using namespace std;
14
15 ClassImp(AliAnalysisMultPbTrackHistoManager)
16
17 const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[]     = { "All Events", "After centrality selection",  "After physics Selection", "With Vertex (quality cuts)", "After ZDC cut", "After vertex range cut (10 cm)" };
18 const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim", 
19                                                                           "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
20 const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[]     = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
21 const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[]     = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake", "hRecHighestMeanPt", "hRecLowestMeanPt"};
22 const char * AliAnalysisMultPbTrackHistoManager::kSpeciesName[]     = { "pi+", "K+", "p", "l+",  "pi-", "K-", "barp", "l-", "Others"};
23
24
25 AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager() : AliHistoListWrapper(), fHNameSuffix(""){ 
26   // standard ctor
27
28 }
29
30 AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const char * name, const char * title): AliHistoListWrapper(name,title), fHNameSuffix("")  {
31   //named ctor
32 };
33
34 AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) : AliHistoListWrapper (obj) {
35   // copy ctor
36   AliError("Not Implemented");
37 };
38
39 AliAnalysisMultPbTrackHistoManager::~AliAnalysisMultPbTrackHistoManager() {
40   // dtor
41
42 }
43
44 TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id, Int_t particle) {
45   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
46
47   TString name = kHistoPtEtaVzNames[id];
48
49   if(particle >= 0) name += kSpeciesName[particle];
50
51   TH3D * h = (TH3D*) GetHisto(name.Data());
52   if (!h) {
53     h = BookHistoPtEtaVz(name.Data(), Form("Pt Eta Vz distribution (%s)",kHistoPtEtaVzNames[id]));
54   }
55
56   return h;
57
58 }
59
60 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
61   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
62
63   TH2D * h = (TH2D*) GetHisto(kHistoDCANames[id]);
64   if (!h) {
65     h = BookHistoDCA(kHistoDCANames[id], Form("DCA vs pt distribution (%s)",kHistoDCANames[id]));
66   }
67
68   return h;
69
70 }
71
72 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoV0vsNtracks(Histo_t id) {
73   // Returns a histo of V0 vs ntracks
74   
75
76   TString name = TString(kHistoPrefix[id])+"_V0vsNtracks";
77
78   TH2D * h =  (TH2D*) GetHisto(name);
79   if (!h) {
80     name+=fHNameSuffix;
81     Bool_t oldStatus = TH1::AddDirectoryStatus();
82     TH1::AddDirectory(kFALSE);
83
84     AliInfo(Form("Booking histo %s",name.Data()));
85
86     h = new TH2D (name.Data(), Form("V0 vs Ntracks (%s)",kHistoPrefix[id]), 300,0,3000, 300, 0, 20000);                  
87     h->SetXTitle("N_{Tracks}");
88     h->SetYTitle("V0 Amplitude");
89     TH1::AddDirectory(oldStatus);
90     fList->Add(h);
91
92
93   }
94   return h;
95   
96
97 }
98
99
100 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMult(Histo_t id) {
101   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
102
103   TString name = kHistoPrefix[id];
104   name += "Mult";
105   TH1D * h = (TH1D*) GetHisto(name.Data());
106   if (!h) {
107     h = BookHistoMult(name.Data(), Form("Multiplicity distribution (%s)",kHistoPrefix[id]));
108   }
109
110   return h;
111
112 }
113
114 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id, 
115                                                        Float_t minEta, Float_t maxEta, 
116                                                        Float_t minVz, Float_t maxVz, 
117                                                        Bool_t scaleWidth) {
118   // Returns a projection of the 3D histo pt/eta/vz.
119   // WARNING: since that is a histo, the requested range will be discretized to the binning.
120   // Always avoids under (over) flows
121   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
122
123   TH3D * h3D = GetHistoPtEtaVz(id);
124
125   // Get range in terms of bin numners.  If the float range is
126   // less than -11111 take the range from the first to the last bin (i.e. no
127   // under/over-flows)
128
129   Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
130   Int_t min2  = minVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;
131
132   Int_t max1 = maxEta  < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
133   Int_t max2  = maxVz  < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
134   // Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
135   // Int_t max2  = maxVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
136
137
138   TString hname = h3D->GetName();
139   hname = hname +  "_pt_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
140
141   
142   if (gROOT->FindObjectAny(hname.Data())){
143     AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
144     hname += "_2";
145   }
146
147   TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E");
148   if(scaleWidth) h ->Scale(1.,"width");
149
150   return h;
151
152 }
153
154
155 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoPtVz (Histo_t id, 
156                                                          Float_t minEta, Float_t maxEta, 
157                                                          Bool_t scaleWidth) {
158   // Returns a projection of the 3D histo pt/eta/vz.
159   // WARNING: since that is a histo, the requested range will be discretized to the binning.
160   // Always avoids under (over) flows
161   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
162
163   // FIXME: what do I do here for the scaling?
164
165   TH3D * h3D = GetHistoPtEtaVz(id);
166
167   // Get range in terms of bin numners.  If the float range is
168   // less than -11111 take the range from the first to the last bin (i.e. no
169   // under/over-flows)
170   Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
171   Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
172
173
174   TString hname = h3D->GetName();
175   hname = hname +  "_ptvz_" + long (min1)  + "_" + long(max1);
176   
177   if (gROOT->FindObjectAny(hname.Data())){
178     AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
179     hname += "_2";
180   }
181
182   h3D->GetYaxis()->SetRange(min1,max1);
183   
184   TH2D * h =  (TH2D*) h3D->Project3D("zxe");
185   if(scaleWidth) h ->Scale(1.,"width");
186
187   return h;
188
189 }
190
191
192 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id, 
193                                                        Float_t minPt, Float_t maxPt,
194                                                        Float_t minEta, Float_t maxEta,
195                                                        Bool_t scaleWidth) { 
196   // Returns a projection of the 3D histo pt/eta/vz.
197   // WARNING: since that is a histo, the requested range will be discretized to the binning.
198   // Always avoids under (over) flows
199   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
200
201   TH3D * h3D = GetHistoPtEtaVz(id);
202
203   // Get range in terms of bin numners.  If the float range is
204   // less than -11111 take the range from the first to the last bin (i.e. no
205   // under/over-flows)
206   Int_t min1  = minPt  < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
207   Int_t min2  = minEta < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minEta);
208
209   Int_t max1  = maxPt  < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
210   Int_t max2  = maxEta < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
211
212
213   TString hname = h3D->GetName();
214   hname = hname +  "_Vz_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
215
216   if (gROOT->FindObjectAny(hname.Data())){
217     AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
218     hname+="_2";
219   }
220
221   TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E");
222   if(scaleWidth) h ->Scale(1.,"width");
223   return h;
224
225
226 }
227
228 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id, 
229                                                         Float_t minPt, Float_t maxPt, 
230                                                         Float_t minVz, Float_t maxVz,
231                                                         Bool_t scaleWidth) {
232   // Returns a projection of the 3D histo pt/eta/vz.
233   // WARNING: since that is a histo, the requested range will be discretized to the binning.
234   // Always avoids under (over) flows
235   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
236
237   TH3D * h3D = GetHistoPtEtaVz(id);
238
239   // Get range in terms of bin numners.  If the float range is
240   // less than -11111 take the range from the first to the last bin (i.e. no
241   // under/over-flows)
242   Int_t min1 = minPt < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
243   Int_t min2 = minVz < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minVz);
244
245   Int_t max1 = maxPt < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
246   Int_t max2 = maxVz < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxVz-0.00001);
247
248   TString hname = h3D->GetName();
249   hname = hname +  "_Eta_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
250
251   if (gROOT->FindObjectAny(hname.Data())){
252     AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
253     hname+="_2";
254   }
255
256   TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E");
257   if(scaleWidth) h ->Scale(1.,"width");
258   return h;
259 }
260
261 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoProcess(Histo_t id) {
262
263   // Returns histogram with particle specties
264
265   TString name = TString(kHistoPrefix[id])+"_Process";
266
267   TH1D * h =  (TH1D*) GetHisto(name);
268   if (!h) {
269     name+=fHNameSuffix;
270     Bool_t oldStatus = TH1::AddDirectoryStatus();
271     TH1::AddDirectory(kFALSE);
272
273     AliInfo(Form("Booking histo %s",name.Data()));
274
275     h = new TH1D (name.Data(), Form("Particle production process (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);                       
276     Int_t nbin = kPNoProcess+1;
277     for(Int_t ibin = 0; ibin < nbin; ibin++){
278       h->GetXaxis()->SetBinLabel(ibin+1,TMCProcessName[ibin]);      
279     }
280     TH1::AddDirectory(oldStatus);
281     fList->Add(h);
282
283
284   }
285   return h;
286   
287
288 }
289
290
291 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoSpecies(Histo_t id) {
292
293   // Returns histogram with particle specties
294
295   TString name = TString(kHistoPrefix[id])+"_Species";
296
297   TH1D * h =  (TH1D*) GetHisto(name);
298   if (!h) {
299     name+=fHNameSuffix;
300     Bool_t oldStatus = TH1::AddDirectoryStatus();
301     TH1::AddDirectory(kFALSE);
302
303     AliInfo(Form("Booking histo %s",name.Data()));
304
305     h = new TH1D (name.Data(), Form("Particle species (%s)",kHistoPrefix[id]), kNPart+1, -0.5, kNPart+1-0.5);                    
306     Int_t nbin = kNPart+1;
307     for(Int_t ibin = 0; ibin < nbin; ibin++){
308       h->GetXaxis()->SetBinLabel(ibin+1,kSpeciesName[ibin]);      
309     }
310     TH1::AddDirectory(oldStatus);
311     fList->Add(h);
312
313
314   }
315   return h;
316   
317
318 }
319
320
321 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMeanPt(Histo_t id) {
322   // mean pt computed event by event
323   TString name = TString(kHistoPrefix[id])+"_MeanPt";
324
325   TH1D * h =  (TH1D*) GetHisto(name);
326   if (!h) {
327     name+=fHNameSuffix;
328     Bool_t oldStatus = TH1::AddDirectoryStatus();
329     TH1::AddDirectory(kFALSE);
330
331     AliInfo(Form("Booking histo %s",name.Data()));
332
333     h = new TH1D (name.Data(), Form("Pt Event by Event(%s)",kHistoPrefix[id]), 100, 0, 1);                       
334     h->Sumw2();
335     TH1::AddDirectory(oldStatus);
336     fList->Add(h);
337
338
339   }
340   return h;
341
342 }
343
344 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEvent(Histo_t id) {
345
346   // Returns of dNdpt, used to compute mean pt event by event
347
348   TString name = TString(kHistoPrefix[id])+"_PtEvent";
349
350   TH1D * h =  (TH1D*) GetHisto(name);
351   if (!h) {
352     name+=fHNameSuffix;
353     Bool_t oldStatus = TH1::AddDirectoryStatus();
354     TH1::AddDirectory(kFALSE);
355
356     AliInfo(Form("Booking histo %s",name.Data()));
357     const Int_t nptbins = 68;
358     const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
359
360     h = new TH1D (name.Data(), Form("Pt Event by Event(%s)",kHistoPrefix[id]), nptbins,binsPt);                  
361     h->SetBit(kKeepMaxMean,1);
362     if(id == kHistoRecLowestMeanPt ) h->SetBit(kKeepMinMean,1);
363     h->Sumw2();
364     TH1::AddDirectory(oldStatus);
365     fList->Add(h);
366
367
368   }
369   return h;
370   
371
372 }
373
374
375 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVzEvent(Histo_t id) {
376
377   // Returns histogram with Vz of the event
378
379   TString name = TString(kHistoPrefix[id])+"_VzEvent";
380
381   TH1D * h =  (TH1D*) GetHisto(name);
382   if (!h) {
383     name+=fHNameSuffix;
384     Bool_t oldStatus = TH1::AddDirectoryStatus();
385     TH1::AddDirectory(kFALSE);
386
387     AliInfo(Form("Booking histo %s",name.Data()));
388
389     h = new TH1D (name.Data(), Form("Vz of the event (%s)",kHistoPrefix[id]), 10, -10, 10);                      
390     h->Sumw2();
391     h->SetXTitle("V_{z}^{event} (cm)");
392     TH1::AddDirectory(oldStatus);
393     fList->Add(h);
394
395
396   }
397   return h;
398   
399
400 }
401
402
403
404 TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() {
405   // Returns histogram with event statistiscs (processed events at each step)
406
407   TH1I * h =  (TH1I*) GetHisto("hStats");
408   if (!h) h = BookHistoStats();
409   return h;
410   
411
412
413
414
415 TH1 * AliAnalysisMultPbTrackHistoManager::GetHisto(const char * name) {
416   //Search list for histo
417   // TODO: keep track of histo id rather than searching by name?
418   return (TH1*) fList->FindObject(TString(name)+fHNameSuffix);
419
420 }
421
422 void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * option) {
423   // Scales all histos in the list for nev
424   // option can be used to pass further options to TH1::Scale
425   TH1 * h = 0;
426   TIter iter = fList->MakeIterator();
427   while ((h = (TH1*) iter.Next())) {
428     if (!h->InheritsFrom("TH1")) {
429       AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName()));
430     }
431     if (h->InheritsFrom("TH1I")) {
432       AliInfo (Form("Not scaling integer histo %s",h->GetName()));
433       continue;
434     }
435     AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
436     
437     h->Scale(1./nev,option);
438   }
439
440 }
441
442 TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, const char * title) {
443   // Books a 3D histo of Pt/eta/vtx
444   // TODO: make the binning settable, variable binning?
445
446   Bool_t oldStatus = TH1::AddDirectoryStatus();
447   TH1::AddDirectory(kFALSE);
448
449   TString hname = name;
450   hname+=fHNameSuffix;
451
452   AliInfo(Form("Booking %s",hname.Data()));
453   
454   // Binning from Jacek task
455   // const Int_t nptbins = 49;
456   // const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0};//,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
457   const Int_t nptbins = 68;
458   const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
459   
460   const Int_t netabins=20;
461   Double_t binsEta[netabins+1];
462   Float_t minEta = -1;
463   Float_t maxEta =  1;
464   Float_t etaStep = (maxEta-minEta)/netabins;
465   for(Int_t ibin = 0; ibin < netabins; ibin++){    
466     binsEta[ibin]   = minEta + ibin*etaStep;
467     binsEta[ibin+1] = minEta + ibin*etaStep + etaStep;
468   }
469
470   const Int_t nvzbins=10;
471   Double_t binsVz[nvzbins+1];
472   Float_t minVz = -10;
473   Float_t maxVz =  10;
474   Float_t vzStep = (maxVz-minVz)/nvzbins;
475   for(Int_t ibin = 0; ibin < nvzbins; ibin++){    
476     binsVz[ibin]   = minVz + ibin*vzStep;
477     binsVz[ibin+1] = minVz + ibin*vzStep + vzStep;
478   }
479  
480
481   TH3D * h = new TH3D (hname,title, 
482                        nptbins,  binsPt,
483                        netabins, binsEta,
484                        nvzbins,  binsVz
485                        );
486
487   h->SetXTitle("p_{T} (GeV)");
488   h->SetYTitle("#eta");
489   h->SetZTitle("V_{z}^{tracks} (cm)");
490   h->Sumw2();
491   
492   fList->Add(h);
493
494   TH1::AddDirectory(oldStatus);
495   return h;
496 }
497
498 TH2D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
499   // Books a DCA histo 
500
501   const Int_t nptbins = 68;
502   const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
503   const Int_t ndcabins=100;
504   Double_t binsDca[ndcabins+1];
505   Float_t minDca = -3;
506   Float_t maxDca =  3;
507   //  const Double_t binsDCA[] = {-3,-2.8,-2.6,-2.4,-2.2,-2.0,-1.9,};
508   Float_t dcaStep = (maxDca-minDca)/ndcabins;
509   for(Int_t ibin = 0; ibin < ndcabins; ibin++){    
510     binsDca[ibin]   = minDca + ibin*dcaStep;
511     binsDca[ibin+1] = minDca + ibin*dcaStep + dcaStep;
512   }
513  
514
515
516
517   Bool_t oldStatus = TH1::AddDirectoryStatus();
518   TH1::AddDirectory(kFALSE);
519
520   TString hname = name;
521   hname+=fHNameSuffix;
522
523   AliInfo(Form("Booking %s",hname.Data()));
524   
525
526 #if defined WEIGHTED_DCA
527   TH1D * h = new TH1D (hname,title, 200,0,200);
528   h->SetXTitle("#Delta DCA");
529 #elif defined TRANSVERSE_DCA
530   TH2D * h = new TH2D (hname,title, ndcabins,binsDca,nptbins,binsPt);
531   h->SetYTitle("p_{T} (GeV)");
532   h->SetXTitle("d_{0} r#phi (cm)");
533 #endif
534   h->Sumw2();
535   
536   fList->Add(h);
537
538   TH1::AddDirectory(oldStatus);
539   return h;
540 }
541 TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoMult(const char * name, const char * title) {
542   // Books a multiplicity histo 
543
544   Bool_t oldStatus = TH1::AddDirectoryStatus();
545   TH1::AddDirectory(kFALSE);
546
547   TString hname = name;
548   hname+=fHNameSuffix;
549
550   AliInfo(Form("Booking %s",hname.Data()));
551   
552
553   TH1D * h = new TH1D (hname,title, 600, 0,6000);
554
555   h->SetXTitle("N tracks");
556   h->Sumw2();
557   
558   fList->Add(h);
559
560   TH1::AddDirectory(oldStatus);
561   return h;
562 }
563
564 TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() {
565   // Books histogram with event statistiscs (processed events at each step)
566
567   Bool_t oldStatus = TH1::AddDirectoryStatus();
568   TH1::AddDirectory(kFALSE);
569
570   AliInfo(Form("Booking Stat histo"));
571
572   TH1I * h = new TH1I (TString("hStats")+fHNameSuffix, "Number of processed events", kNStatBins, -0.5, kNStatBins-0.5);
573   for(Int_t istatbin = 0; istatbin < kNStatBins; istatbin++){
574     h->GetXaxis()->SetBinLabel(istatbin+1,kStatStepNames[istatbin]);
575   }
576   TH1::AddDirectory(oldStatus);
577   fList->Add(h);
578   return h;
579 }
580
581 Int_t AliAnalysisMultPbTrackHistoManager::GetLocalParticleID(AliMCParticle * part) {
582   // returns the local code (Part_t)
583   
584   Int_t pdgcode = part->PdgCode();
585   switch(pdgcode) {
586
587   case 211:
588     return kPartPiPlus;
589     break;
590   case -211:
591     return kPartPiMinus;
592     break;
593   case 2212:
594     return kPartP;
595     break;
596   case -2212:
597     return kPartPBar;
598     break;
599   case 321:
600     return kPartKPlus;
601     break;
602   case -321:
603     return kPartKMinus;
604     break;
605   case -11:
606     return kPartLMinus;
607     break;
608   case 11:
609     return kPartLPlus;
610     break;
611   case -13:
612     return kPartLMinus;
613     break;
614   case 13:
615     return kPartLPlus;
616     break;
617   default:
618     return kPartOther;
619   }
620 }
621
622
623  
624 Long64_t AliAnalysisMultPbTrackHistoManager::Merge(TCollection* list)
625 {
626   // Merge a list of AliHistoListWrapper objects with this.
627   // Returns the number of merged objects (including this).
628
629   // We have to make sure that all the list contain the same histos in
630   // the same order. We thus also have to sort the list (sorting is
631   // done by name in TList).
632
633   //  AliInfo("Merging");
634
635   if (!list)
636     return 0;
637
638   if (list->IsEmpty())
639     return 1;
640
641   TIterator* iter = list->MakeIterator();
642   TObject* obj;
643
644   // collections of all histograms
645   TList collections;
646
647   Int_t count = 0;
648
649   while ((obj = iter->Next())) {
650     Bool_t foundDiffinThisIterStep = kFALSE;
651
652     //    Printf("%d - %s",count, obj->GetName());
653     AliHistoListWrapper* entry = dynamic_cast<AliHistoListWrapper*> (obj);
654     if (entry == 0) 
655       continue;
656
657     TList * hlist = entry->GetList();
658
659     // Check if all histos in this fList are also in the one from entry and viceversa
660     // Use getters to automatically book non defined histos    
661
662     Bool_t areListsDifferent=kTRUE;
663     Int_t iloop = 0;
664     Int_t max_loops = hlist->GetSize() + fList->GetSize(); // In the worst case all of the histos will be different...
665     while(areListsDifferent) {
666       if(iloop>max_loops) AliFatal("Infinite Loop?");
667       iloop++;
668       // sort
669       hlist->Sort();
670       fList->Sort();
671       // loop over the largest 
672       TObject * hist =0;
673       TIterator * iterlist = 0;
674       TList * thislist  = 0; // the list over which I'm iterating
675       TList * otherlist = 0; // the other
676
677       if (hlist->GetSize() >= fList->GetSize()) { 
678         thislist  = hlist;
679         otherlist = fList;
680       }
681       else{
682         thislist  = fList;
683         otherlist = hlist;      
684       }
685       iterlist = thislist->MakeIterator();
686
687       while ((hist= iterlist->Next())){ 
688         if(!otherlist->FindObject(hist->GetName())){
689           AliInfo(Form("Adding object %s",hist->GetName()));      
690           TH1 * hclone =  (TH1*) hist->Clone();
691           if (!hclone->InheritsFrom("TH1")) AliFatal(Form("Found a %s. This class only supports objects inheriting from TH1",hclone->ClassName()));
692           hclone->Reset();
693           otherlist->Add(hclone);
694           foundDiffinThisIterStep=kTRUE;
695         }
696       }
697
698       // re-sort before checking
699       hlist->Sort();
700       fList->Sort();
701
702       // check if everything is fine    
703       areListsDifferent=kFALSE;
704       if (hlist->GetSize() == fList->GetSize()) {       
705         Int_t nhist =  fList->GetSize();
706         for(Int_t ihist = 0; ihist < nhist; ihist++){
707           if(strcmp(fList->At(ihist)->GetName(),hlist->At(ihist)->GetName())) areListsDifferent = kTRUE;
708         }
709       } else {
710         areListsDifferent=kTRUE;
711       }
712     }
713
714     // last check: if something is not ok die loudly 
715     if (hlist->GetSize() != fList->GetSize()) {
716       AliFatal("Mismatching size!");
717     }
718     Int_t nhist =  fList->GetSize();
719     for(Int_t ihist = 0; ihist < nhist; ihist++){
720       if(strcmp(fList->At(ihist)->GetName(),hlist->At(ihist)->GetName())){
721         AliFatal(Form("Mismatching histos: %s -> %s", fList->At(ihist)->GetName(),hlist->At(ihist)->GetName()));
722       } else {
723         // Specific merging strategies, according to the bits...
724         if (fList->At(ihist)->TestBits(kKeepMinMean|kKeepMaxMean)){
725           TH1 * h1 = (TH1*) fList->At(ihist);
726           TH1 * h2 = (TH1*) hlist->At(ihist);
727           if (h1->GetEntries()>0 && h2->GetEntries()>0) {
728             if (h1->TestBit(kKeepMinMean)) {
729               AliInfo(Form("Keeping only the histo with the lowest mean [%s][%f][%f]",h1->GetName(), h1->GetMean(), h2->GetMean()));
730               if(h1->GetMean() > h2->GetMean()) h1->Reset();
731               else h2->Reset();
732             }
733             if (h1->TestBit(kKeepMaxMean)) {
734               AliInfo(Form("Keeping only the histo with the highest mean [%s][%f][%f]",h1->GetName(), h1->GetMean(), h2->GetMean()));
735               if(h2->GetMean() > h1->GetMean()) h1->Reset();
736               else h2->Reset();
737             }
738           }
739         }
740       }
741     }
742     
743     if (foundDiffinThisIterStep){
744       iter->Reset(); // We found a difference: previous lists could
745                      // also be affected... We start from scratch
746       collections.Clear();
747       count = 0;
748     }
749     else {
750       
751       collections.Add(hlist);
752       
753       count++;
754     }
755   }
756
757   fList->Merge(&collections);
758   
759   delete iter;
760
761   return count+1;
762 }