Moving PbPb multiplicity in the new directory structure
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / multPbPb / AliAnalysisMultPbTrackHistoManager.cxx
CommitLineData
a23f7c97 1#include "AliAnalysisMultPbTrackHistoManager.h"
2#include "AliLog.h"
3#include "TH1.h"
4#include "TH3D.h"
5#include "TH1I.h"
6#include "TROOT.h"
e0376287 7#include "TMCProcess.h"
3b8cbf2d 8#include "AliMCParticle.h"
e0376287 9
a23f7c97 10#include <iostream>
aa6151eb 11#include "TH2D.h"
3b8cbf2d 12
a23f7c97 13using namespace std;
3b8cbf2d 14
a23f7c97 15ClassImp(AliAnalysisMultPbTrackHistoManager)
16
abd808b9 17const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[] = { "All Events", "After centrality selection", "After physics Selection", "With Vertex (quality cuts)", "After ZDC cut", "After vertex range cut (10 cm)" };
a23f7c97 18const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim",
19 "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
20const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[] = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
abd808b9 21const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[] = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake", "hRecHighestMeanPt", "hRecLowestMeanPt"};
4d0aa70f 22const char * AliAnalysisMultPbTrackHistoManager::kSpeciesName[] = { "pi+", "K+", "p", "l+", "pi-", "K-", "barp", "l-", "Others"};
a23f7c97 23
24
25AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager() : AliHistoListWrapper(), fHNameSuffix(""){
26 // standard ctor
27
28}
29
30AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const char * name, const char * title): AliHistoListWrapper(name,title), fHNameSuffix("") {
31 //named ctor
32};
33
34AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) : AliHistoListWrapper (obj) {
35 // copy ctor
36 AliError("Not Implemented");
37};
38
39AliAnalysisMultPbTrackHistoManager::~AliAnalysisMultPbTrackHistoManager() {
40 // dtor
41
42}
43
4d0aa70f 44TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id, Int_t particle) {
a23f7c97 45 // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
46
4d0aa70f 47 TString name = kHistoPtEtaVzNames[id];
48
49 if(particle >= 0) name += kSpeciesName[particle];
50
51 TH3D * h = (TH3D*) GetHisto(name.Data());
a23f7c97 52 if (!h) {
4d0aa70f 53 h = BookHistoPtEtaVz(name.Data(), Form("Pt Eta Vz distribution (%s)",kHistoPtEtaVzNames[id]));
a23f7c97 54 }
55
56 return h;
57
58}
59
bcc49144 60TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
a23f7c97 61 // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
62
bcc49144 63 TH2D * h = (TH2D*) GetHisto(kHistoDCANames[id]);
a23f7c97 64 if (!h) {
bcf2601a 65 h = BookHistoDCA(kHistoDCANames[id], Form("DCA vs pt distribution (%s)",kHistoDCANames[id]));
a23f7c97 66 }
67
68 return h;
69
70}
71
bcf2601a 72TH2D * 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);
bfae2924 87 h->SetXTitle("N_{Tracks}");
88 h->SetYTitle("V0 Amplitude");
bcf2601a 89 TH1::AddDirectory(oldStatus);
90 fList->Add(h);
91
92
93 }
94 return h;
95
96
97}
98
99
e0376287 100TH1D * 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
a23f7c97 114TH1D * 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)
bcc49144 128
4d0aa70f 129 Int_t min1 = minEta < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
130 Int_t min2 = minVz < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;
a23f7c97 131
bcc49144 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);
a23f7c97 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())){
e0376287 143 AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
144 hname += "_2";
a23f7c97 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
aa6151eb 154
155TH2D * 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
a23f7c97 192TH1D * 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())){
e0376287 217 AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
218 hname+="_2";
a23f7c97 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
228TH1D * 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())){
e0376287 252 AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
253 hname+="_2";
a23f7c97 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
3b8cbf2d 261TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoProcess(Histo_t id) {
e0376287 262
263 // Returns histogram with particle specties
264
3b8cbf2d 265 TString name = TString(kHistoPrefix[id])+"_Process";
e0376287 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
3b8cbf2d 275 h = new TH1D (name.Data(), Form("Particle production process (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);
e0376287 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
3b8cbf2d 291TH1D * 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
abd808b9 320
321TH1D * 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
344TH1D * 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
4d0aa70f 375TH1D * 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
3b8cbf2d 402
a23f7c97 403
404TH1I * 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
415TH1 * 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
422void 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 }
9441cdaa 431 if (h->InheritsFrom("TH1I")) {
432 AliInfo (Form("Not scaling integer histo %s",h->GetName()));
433 continue;
434 }
e0376287 435 AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
a23f7c97 436
437 h->Scale(1./nev,option);
438 }
439
440}
441
442TH3D * 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
e0376287 454 // Binning from Jacek task
eef42d18 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};
e0376287 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
a23f7c97 480
481 TH3D * h = new TH3D (hname,title,
e0376287 482 nptbins, binsPt,
483 netabins, binsEta,
484 nvzbins, binsVz
a23f7c97 485 );
486
4d0aa70f 487 h->SetXTitle("p_{T} (GeV)");
e0376287 488 h->SetYTitle("#eta");
4d0aa70f 489 h->SetZTitle("V_{z}^{tracks} (cm)");
a23f7c97 490 h->Sumw2();
491
492 fList->Add(h);
493
494 TH1::AddDirectory(oldStatus);
495 return h;
496}
497
bcc49144 498TH2D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
a23f7c97 499 // Books a DCA histo
500
bcc49144 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
a23f7c97 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
bcc49144 526#if defined WEIGHTED_DCA
54b31d6a 527 TH1D * h = new TH1D (hname,title, 200,0,200);
a23f7c97 528 h->SetXTitle("#Delta DCA");
bcc49144 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
a23f7c97 534 h->Sumw2();
535
536 fList->Add(h);
537
538 TH1::AddDirectory(oldStatus);
539 return h;
540}
e0376287 541TH1D * 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}
a23f7c97 563
564TH1I * 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
3b8cbf2d 581Int_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}
a23f7c97 621
622
623
abd808b9 624Long64_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}