]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/ChargedHadrons/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[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"
5ad99723 9#include "TMap.h"
a23f7c97 10#include <iostream>
aa6151eb 11#include "TH2D.h"
5ad99723 12#include "TDatabasePDG.h"
13#include "TParameter.h"
3b8cbf2d 14
a23f7c97 15using namespace std;
3b8cbf2d 16
a23f7c97 17ClassImp(AliAnalysisMultPbTrackHistoManager)
18
abd808b9 19const 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 20const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim",
21 "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
22const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[] = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
abd808b9 23const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[] = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake", "hRecHighestMeanPt", "hRecLowestMeanPt"};
4d0aa70f 24const char * AliAnalysisMultPbTrackHistoManager::kSpeciesName[] = { "pi+", "K+", "p", "l+", "pi-", "K-", "barp", "l-", "Others"};
a23f7c97 25
26
5ad99723 27AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager() : AliHistoListWrapper(), fHNameSuffix(""), fParticleSpecies(0){
a23f7c97 28 // standard ctor
29
30}
31
5ad99723 32AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const char * name, const char * title): AliHistoListWrapper(name,title), fHNameSuffix(""), fParticleSpecies(0) {
a23f7c97 33 //named ctor
34};
35
36AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) : AliHistoListWrapper (obj) {
37 // copy ctor
38 AliError("Not Implemented");
39};
40
41AliAnalysisMultPbTrackHistoManager::~AliAnalysisMultPbTrackHistoManager() {
42 // dtor
43
44}
45
5ad99723 46TH2D * 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
4d0aa70f 66TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id, Int_t particle) {
a23f7c97 67 // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
68
4d0aa70f 69 TString name = kHistoPtEtaVzNames[id];
70
71 if(particle >= 0) name += kSpeciesName[particle];
72
73 TH3D * h = (TH3D*) GetHisto(name.Data());
a23f7c97 74 if (!h) {
4d0aa70f 75 h = BookHistoPtEtaVz(name.Data(), Form("Pt Eta Vz distribution (%s)",kHistoPtEtaVzNames[id]));
a23f7c97 76 }
77
78 return h;
79
80}
81
bcc49144 82TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
a23f7c97 83 // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
84
bcc49144 85 TH2D * h = (TH2D*) GetHisto(kHistoDCANames[id]);
a23f7c97 86 if (!h) {
bcf2601a 87 h = BookHistoDCA(kHistoDCANames[id], Form("DCA vs pt distribution (%s)",kHistoDCANames[id]));
a23f7c97 88 }
89
90 return h;
91
92}
93
bcf2601a 94TH2D * 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);
bfae2924 109 h->SetXTitle("N_{Tracks}");
110 h->SetYTitle("V0 Amplitude");
bcf2601a 111 TH1::AddDirectory(oldStatus);
112 fList->Add(h);
113
114
115 }
116 return h;
117
118
119}
120
121
e0376287 122TH1D * 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
a23f7c97 136TH1D * 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)
bcc49144 150
4d0aa70f 151 Int_t min1 = minEta < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
152 Int_t min2 = minVz < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;
a23f7c97 153
bcc49144 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);
a23f7c97 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())){
e0376287 165 AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
166 hname += "_2";
a23f7c97 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
aa6151eb 176
177TH2D * 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
a23f7c97 214TH1D * 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())){
e0376287 239 AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
240 hname+="_2";
a23f7c97 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
250TH1D * 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())){
e0376287 274 AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
275 hname+="_2";
a23f7c97 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
3b8cbf2d 283TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoProcess(Histo_t id) {
e0376287 284
285 // Returns histogram with particle specties
286
3b8cbf2d 287 TString name = TString(kHistoPrefix[id])+"_Process";
e0376287 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
3b8cbf2d 297 h = new TH1D (name.Data(), Form("Particle production process (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);
e0376287 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
3b8cbf2d 313TH1D * 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
abd808b9 342
343TH1D * 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
366TH1D * 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
4d0aa70f 397TH1D * 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
3b8cbf2d 424
a23f7c97 425
426TH1I * 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
437TH1 * 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
444void 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 }
9441cdaa 453 if (h->InheritsFrom("TH1I")) {
454 AliInfo (Form("Not scaling integer histo %s",h->GetName()));
455 continue;
456 }
e0376287 457 AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
a23f7c97 458
459 h->Scale(1./nev,option);
460 }
461
462}
463
464TH3D * 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
e0376287 476 // Binning from Jacek task
eef42d18 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};
e0376287 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
a23f7c97 502
503 TH3D * h = new TH3D (hname,title,
e0376287 504 nptbins, binsPt,
505 netabins, binsEta,
506 nvzbins, binsVz
a23f7c97 507 );
508
4d0aa70f 509 h->SetXTitle("p_{T} (GeV)");
e0376287 510 h->SetYTitle("#eta");
4d0aa70f 511 h->SetZTitle("V_{z}^{tracks} (cm)");
a23f7c97 512 h->Sumw2();
513
a23f7c97 514 TH1::AddDirectory(oldStatus);
5ad99723 515 fList->Add(h);
a23f7c97 516 return h;
517}
518
bcc49144 519TH2D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
a23f7c97 520 // Books a DCA histo
521
bcc49144 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
a23f7c97 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
bcc49144 547#if defined WEIGHTED_DCA
54b31d6a 548 TH1D * h = new TH1D (hname,title, 200,0,200);
a23f7c97 549 h->SetXTitle("#Delta DCA");
bcc49144 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
a23f7c97 555 h->Sumw2();
556
557 fList->Add(h);
558
559 TH1::AddDirectory(oldStatus);
560 return h;
561}
e0376287 562TH1D * 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}
a23f7c97 584
585TH1I * 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
3b8cbf2d 602Int_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}
a23f7c97 642
643
5ad99723 644void 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
a23f7c97 660
abd808b9 661Long64_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
5ad99723 670 AliInfo("Merging");
abd808b9 671 if (!list)
abd808b9 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())) {
5ad99723 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;
abd808b9 692
5ad99723 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);
abd808b9 729 if (entry == 0)
730 continue;
731
5ad99723 732 Bool_t foundDiffinThisIterStep = kFALSE;
733
734 // Printf("%d - %s",count, obj->GetName());
735
abd808b9 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
5ad99723 840 AliInfo("Merged");
841
abd808b9 842 return count+1;
843}