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