]>
Commit | Line | Data |
---|---|---|
ce85db45 | 1 | //==================================================================== |
2 | #include "AliBasedNdetaTask.h" | |
3 | #include <TMath.h> | |
4 | #include <TH2D.h> | |
5 | #include <TH1D.h> | |
6 | #include <THStack.h> | |
7 | #include <TList.h> | |
8 | #include <AliAnalysisManager.h> | |
9 | #include <AliAODEvent.h> | |
10 | #include <AliAODHandler.h> | |
11 | #include <AliAODInputHandler.h> | |
12 | #include "AliForwardUtil.h" | |
13 | #include "AliAODForwardMult.h" | |
3846ec25 | 14 | #include <TFile.h> |
5bb5d1f6 | 15 | #include <TStyle.h> |
c8b1a7db | 16 | #include <TROOT.h> |
ce85db45 | 17 | |
18 | //____________________________________________________________________ | |
19 | AliBasedNdetaTask::AliBasedNdetaTask() | |
c8b1a7db | 20 | : AliBaseAODTask(), |
ce85db45 | 21 | fRebin(0), // Rebinning factor |
22 | fCutEdges(false), | |
23 | fSymmetrice(true), | |
e58000b7 | 24 | fCorrEmpty(true), |
c25b5e1b | 25 | fUseROOTProj(false), |
e58000b7 | 26 | fTriggerEff(1), |
e12f8c42 | 27 | fTriggerEff0(1), |
bfab35d9 | 28 | fShapeCorr(0), |
29 | fListOfCentralities(0), | |
bfab35d9 | 30 | fNormalizationScheme(kFull), |
3d9b0442 | 31 | fFinalMCCorrFile(""), |
bfab35d9 | 32 | fSatelliteVertices(0), |
cd49c50b | 33 | fglobalempiricalcorrection(0), |
bfab35d9 | 34 | fmeabsignalvscentr(0) |
fb3430ac | 35 | { |
36 | // | |
37 | // Constructor | |
38 | // | |
5ca83fee | 39 | DGUARD(fDebug,3,"Default CTOR of AliBasedNdetaTask"); |
fb3430ac | 40 | } |
ce85db45 | 41 | |
42 | //____________________________________________________________________ | |
43 | AliBasedNdetaTask::AliBasedNdetaTask(const char* name) | |
c8b1a7db | 44 | : AliBaseAODTask(Form("%sdNdeta", name)), |
ce85db45 | 45 | fRebin(5), // Rebinning factor |
46 | fCutEdges(false), | |
47 | fSymmetrice(true), | |
e58000b7 | 48 | fCorrEmpty(true), |
c25b5e1b | 49 | fUseROOTProj(false), |
e58000b7 | 50 | fTriggerEff(1), |
66cf95f2 | 51 | fTriggerEff0(1), |
3846ec25 | 52 | fShapeCorr(0), |
e1f47419 | 53 | fListOfCentralities(0), |
0be6c8cd | 54 | fNormalizationScheme(kFull), |
3d9b0442 | 55 | fFinalMCCorrFile(""), |
bfab35d9 | 56 | fSatelliteVertices(0), |
cd49c50b | 57 | fglobalempiricalcorrection(0), |
58 | fmeabsignalvscentr(0) | |
ce85db45 | 59 | { |
fb3430ac | 60 | // |
61 | // Constructor | |
62 | // | |
5ca83fee | 63 | DGUARD(fDebug, 3,"Named CTOR of AliBasedNdetaTask: %s", name); |
c8b1a7db | 64 | |
65 | fTriggerMask = AliAODForwardMult::kInel; | |
66 | fMinIpZ = -10; | |
67 | fMaxIpZ = +10; | |
e308a636 | 68 | fListOfCentralities = new TObjArray(1); |
0be6c8cd | 69 | |
70 | // Set the normalisation scheme | |
71 | SetNormalizationScheme(kFull); | |
72 | ||
ce85db45 | 73 | } |
74 | ||
ce85db45 | 75 | |
76 | //____________________________________________________________________ | |
77 | AliBasedNdetaTask::~AliBasedNdetaTask() | |
78 | { | |
fb3430ac | 79 | // |
80 | // Destructor | |
81 | // | |
f53fb4f6 | 82 | DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask"); |
ce85db45 | 83 | } |
84 | ||
f53fb4f6 | 85 | //________________________________________________________________________ |
86 | void | |
87 | AliBasedNdetaTask::SetDebugLevel(Int_t lvl) | |
88 | { | |
89 | AliAnalysisTaskSE::SetDebugLevel(lvl); | |
90 | for (Int_t i = 0; i < fListOfCentralities->GetEntries(); i++) { | |
91 | CentralityBin* bin = | |
92 | static_cast<CentralityBin*>(fListOfCentralities->At(i)); | |
93 | bin->SetDebugLevel(lvl); | |
94 | } | |
95 | } | |
96 | ||
e308a636 | 97 | //________________________________________________________________________ |
98 | void | |
99 | AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high) | |
e1f47419 | 100 | { |
101 | // | |
102 | // Add a centrality bin | |
103 | // | |
104 | // Parameters: | |
105 | // low Low cut | |
106 | // high High cut | |
107 | // | |
1a7e2e15 | 108 | DGUARD(fDebug,3,"Add a centrality bin [%d,%d] @ %d", low, high, at); |
e308a636 | 109 | CentralityBin* bin = MakeCentralityBin(GetName(), low, high); |
1a7e2e15 | 110 | if (!bin) { |
111 | Error("AddCentralityBin", | |
112 | "Failed to create centrality bin for %s [%d,%d] @ %d", | |
113 | GetName(), low, high, at); | |
114 | return; | |
115 | } | |
bfab35d9 | 116 | bin->SetSatelliteVertices(fSatelliteVertices); |
f53fb4f6 | 117 | bin->SetDebugLevel(fDebug); |
e308a636 | 118 | fListOfCentralities->AddAtAndExpand(bin, at); |
e1f47419 | 119 | } |
120 | ||
121 | //________________________________________________________________________ | |
122 | AliBasedNdetaTask::CentralityBin* | |
123 | AliBasedNdetaTask::MakeCentralityBin(const char* name, | |
124 | Short_t low, Short_t high) const | |
125 | { | |
126 | // | |
127 | // Make a centrality bin | |
128 | // | |
129 | // Parameters: | |
130 | // name Name used for histograms | |
131 | // low Low cut in percent | |
132 | // high High cut in percent | |
133 | // | |
134 | // Return: | |
135 | // A newly created centrality bin | |
136 | // | |
1a7e2e15 | 137 | DGUARD(fDebug,3,"Make a centrality bin %s [%d,%d]", name, low, high); |
e1f47419 | 138 | return new CentralityBin(name, low, high); |
139 | } | |
a76fb27d | 140 | |
141 | #define TESTAPPEND(SCHEME,BIT,STRING) \ | |
142 | do { if (!(SCHEME & BIT)) break; \ | |
143 | if (!s.IsNull()) s.Append(","); s.Append(STRING); } while(false) | |
144 | ||
ffca499d | 145 | //________________________________________________________________________ |
a76fb27d | 146 | const Char_t* |
147 | AliBasedNdetaTask::NormalizationSchemeString(UShort_t scheme) | |
148 | { | |
149 | // Create a string from normalization scheme bits | |
150 | static TString s; | |
151 | s = ""; | |
152 | ||
153 | if (scheme == kNone) | |
154 | return s.Data(); | |
155 | if (scheme == kFull) { | |
156 | s = "FULL"; | |
157 | return s.Data(); | |
158 | } | |
159 | TESTAPPEND(scheme, kEventLevel, "EVENT"); | |
160 | TESTAPPEND(scheme, kShape, "SHAPE"); | |
161 | TESTAPPEND(scheme, kBackground, "BACKGROUND"); | |
162 | TESTAPPEND(scheme, kTriggerEfficiency, "TRIGGER"); | |
163 | TESTAPPEND(scheme, kZeroBin, "ZEROBIN"); | |
164 | ||
165 | return s.Data(); | |
166 | } | |
167 | //________________________________________________________________________ | |
168 | UShort_t | |
169 | AliBasedNdetaTask::ParseNormalizationScheme(const char* what) | |
ffca499d | 170 | { |
ffca499d | 171 | UShort_t scheme = 0; |
172 | TString twhat(what); | |
173 | twhat.ToUpper(); | |
174 | TObjString* opt; | |
09d5920f | 175 | TObjArray* token = twhat.Tokenize(" ,|"); |
176 | TIter next(token); | |
ffca499d | 177 | while ((opt = static_cast<TObjString*>(next()))) { |
178 | TString s(opt->GetString()); | |
179 | if (s.IsNull()) continue; | |
0be6c8cd | 180 | Bool_t add = true; |
181 | switch (s[0]) { | |
182 | case '-': add = false; // Fall through | |
183 | case '+': s.Remove(0,1); // Remove character | |
184 | } | |
185 | UShort_t bit = 0; | |
186 | if (s.CompareTo("EVENT") == 0) bit = kEventLevel; | |
187 | else if (s.CompareTo("SHAPE") == 0) bit = kShape; | |
188 | else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground; | |
189 | else if (s.CompareTo("TRIGGER") == 0) bit = kTriggerEfficiency; | |
190 | else if (s.CompareTo("FULL") == 0) bit = kFull; | |
191 | else if (s.CompareTo("NONE") == 0) bit = kNone; | |
4fa8d795 | 192 | else if (s.CompareTo("ZEROBIN") == 0) bit = kZeroBin; |
ffca499d | 193 | else |
a76fb27d | 194 | ::Warning("SetNormalizationScheme", "Unknown option %s", s.Data()); |
0be6c8cd | 195 | if (add) scheme |= bit; |
196 | else scheme ^= bit; | |
ffca499d | 197 | } |
09d5920f | 198 | delete token; |
a76fb27d | 199 | return scheme; |
200 | } | |
201 | //________________________________________________________________________ | |
202 | void | |
203 | AliBasedNdetaTask::SetNormalizationScheme(const char* what) | |
204 | { | |
205 | // | |
206 | // Set normalisation scheme | |
207 | // | |
208 | DGUARD(fDebug,3,"Set the normalization scheme: %s", what); | |
209 | SetNormalizationScheme(ParseNormalizationScheme(what)); | |
ffca499d | 210 | } |
0be6c8cd | 211 | //________________________________________________________________________ |
212 | void | |
213 | AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme) | |
214 | { | |
f53fb4f6 | 215 | DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme); |
0be6c8cd | 216 | fNormalizationScheme = scheme; |
0be6c8cd | 217 | } |
e58000b7 | 218 | //________________________________________________________________________ |
219 | void | |
f67d699c | 220 | AliBasedNdetaTask::SetShapeCorrection(const TH2F* c) |
e58000b7 | 221 | { |
fb3430ac | 222 | // |
223 | // Set the shape correction (a.k.a., track correction) for selected | |
224 | // trigger(s) | |
225 | // | |
226 | // Parameters: | |
227 | // h Correction | |
228 | // | |
f53fb4f6 | 229 | DGUARD(fDebug,3,"Set the shape correction: %p", c); |
e58000b7 | 230 | if (!c) return; |
f67d699c | 231 | fShapeCorr = static_cast<TH2F*>(c->Clone()); |
e58000b7 | 232 | fShapeCorr->SetDirectory(0); |
233 | } | |
f67d699c | 234 | //________________________________________________________________________ |
235 | void | |
236 | AliBasedNdetaTask::InitializeCentBins() | |
237 | { | |
238 | if (fListOfCentralities->GetEntries() > 0) return; | |
239 | ||
240 | // Automatically add 'all' centrality bin if nothing has been defined. | |
241 | AddCentralityBin(0, 0, 0); | |
c8b1a7db | 242 | if (HasCentrality() && fCentAxis.GetXbins()) { |
243 | const TArrayD* bins = fCentAxis.GetXbins(); | |
244 | Int_t nbin = fCentAxis.GetNbins(); | |
f67d699c | 245 | for (Int_t i = 0; i < nbin; i++) |
246 | AddCentralityBin(i+1, Short_t((*bins)[i]), Short_t((*bins)[i+1])); | |
247 | } | |
248 | } | |
249 | ||
ce85db45 | 250 | //________________________________________________________________________ |
c8b1a7db | 251 | Bool_t |
252 | AliBasedNdetaTask::Book() | |
ce85db45 | 253 | { |
fb3430ac | 254 | // |
255 | // Create output objects. | |
256 | // | |
257 | // This is called once per slave process | |
258 | // | |
f53fb4f6 | 259 | DGUARD(fDebug,1,"Create user ouput object"); |
ce85db45 | 260 | |
bfab35d9 | 261 | fSums->Add(AliForwardUtil::MakeParameter("empirical", |
262 | fglobalempiricalcorrection != 0)); | |
c8b1a7db | 263 | fSums->Add(AliForwardUtil::MakeParameter("scheme", fNormalizationScheme)); |
264 | ||
265 | // Make our centrality bins | |
266 | InitializeCentBins(); | |
bfab35d9 | 267 | |
e1f47419 | 268 | // Loop over centrality bins |
269 | TIter next(fListOfCentralities); | |
270 | CentralityBin* bin = 0; | |
271 | while ((bin = static_cast<CentralityBin*>(next()))) | |
1a7e2e15 | 272 | bin->CreateOutputObjects(fSums, fTriggerMask); |
f67d699c | 273 | |
bfab35d9 | 274 | fmeabsignalvscentr=new TH2D("meanAbsSignalVsCentr", |
275 | "Mean absolute signal versus centrality", | |
276 | 400, 0, 20, 100, 0, 100); | |
277 | fSums->Add(fmeabsignalvscentr); | |
278 | ||
c8b1a7db | 279 | return true; |
ce85db45 | 280 | } |
c8b1a7db | 281 | |
ce85db45 | 282 | //____________________________________________________________________ |
c8b1a7db | 283 | Bool_t |
284 | AliBasedNdetaTask::CheckEvent(const AliAODForwardMult& fwd) | |
285 | { | |
286 | AliBaseAODTask::CheckEvent(fwd); | |
287 | // Here, we always return true, as the centrality bins will do | |
288 | // their own checks on the events - this is needed for event | |
289 | // normalization etc. | |
290 | return true; | |
291 | } | |
292 | ||
293 | //____________________________________________________________________ | |
294 | Bool_t | |
295 | AliBasedNdetaTask::Event(AliAODEvent& aod) | |
ce85db45 | 296 | { |
fb3430ac | 297 | // |
298 | // Process a single event | |
299 | // | |
300 | // Parameters: | |
301 | // option Not used | |
302 | // | |
ce85db45 | 303 | // Main loop |
f53fb4f6 | 304 | DGUARD(fDebug,1,"Analyse the AOD event"); |
bfab35d9 | 305 | |
c8b1a7db | 306 | AliAODForwardMult* forward = GetForward(aod); |
307 | if (!forward) return false; | |
e308a636 | 308 | |
e1f47419 | 309 | // Fill centrality histogram |
bfab35d9 | 310 | |
c8b1a7db | 311 | Double_t vtx = forward->GetIpZ(); |
312 | TH2D* data = GetHistogram(aod, false); | |
313 | TH2D* dataMC = GetHistogram(aod, true); | |
314 | if (!data) return false; | |
bfab35d9 | 315 | |
316 | CheckEventData(vtx, data, dataMC); | |
317 | ||
318 | if (!ApplyEmpiricalCorrection(forward,data)) | |
c8b1a7db | 319 | return false; |
cd49c50b | 320 | |
cd49c50b | 321 | |
4fa8d795 | 322 | Bool_t isZero = ((fNormalizationScheme & kZeroBin) && |
323 | !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0)); | |
bfab35d9 | 324 | Bool_t taken = false; |
325 | ||
e1f47419 | 326 | // Loop over centrality bins |
ffca499d | 327 | CentralityBin* allBin = |
328 | static_cast<CentralityBin*>(fListOfCentralities->At(0)); | |
bfab35d9 | 329 | if (allBin->ProcessEvent(forward, fTriggerMask, isZero, |
c8b1a7db | 330 | fMinIpZ, fMaxIpZ, data, dataMC)) taken = true; |
e308a636 | 331 | |
332 | // Find this centrality bin | |
c8b1a7db | 333 | if (HasCentrality()) { |
334 | Double_t cent = forward->GetCentrality(); | |
335 | Int_t icent = fCentAxis.FindBin(cent); | |
e308a636 | 336 | CentralityBin* thisBin = 0; |
c8b1a7db | 337 | if (icent >= 1 && icent <= fCentAxis.GetNbins()) |
e308a636 | 338 | thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent)); |
339 | if (thisBin) | |
c8b1a7db | 340 | if (thisBin->ProcessEvent(forward, fTriggerMask, isZero, fMinIpZ, |
341 | fMaxIpZ, data, dataMC)) taken = true; | |
e308a636 | 342 | } |
c8b1a7db | 343 | |
344 | return taken; | |
ce85db45 | 345 | } |
346 | ||
1f7aa5c7 | 347 | //________________________________________________________________________ |
348 | void AliBasedNdetaTask::CheckEventData(Double_t, | |
349 | TH2*, | |
350 | TH2*) | |
351 | { | |
352 | } | |
353 | ||
ce85db45 | 354 | //________________________________________________________________________ |
355 | void | |
356 | AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker, | |
357 | const char* title, const char* ytitle) | |
358 | { | |
fb3430ac | 359 | // |
360 | // Set histogram graphical options, etc. | |
361 | // | |
362 | // Parameters: | |
363 | // h Histogram to modify | |
364 | // colour Marker color | |
365 | // marker Marker style | |
366 | // title Title of histogram | |
367 | // ytitle Title on y-axis. | |
368 | // | |
ce85db45 | 369 | h->SetTitle(title); |
370 | h->SetMarkerColor(colour); | |
371 | h->SetMarkerStyle(marker); | |
9ecab72f | 372 | h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1); |
ce85db45 | 373 | h->SetFillStyle(0); |
c8b1a7db | 374 | TString ytit; |
375 | if (ytitle && ytitle[0] != '\0') ytit = ytitle; | |
376 | ytit = "#frac{1}{#it{N}}#frac{d#it{N}_{ch}}{d#it{#eta}}"; | |
377 | h->SetYTitle(ytit); | |
ce85db45 | 378 | h->GetXaxis()->SetTitleFont(132); |
379 | h->GetXaxis()->SetLabelFont(132); | |
380 | h->GetXaxis()->SetNdivisions(10); | |
381 | h->GetYaxis()->SetTitleFont(132); | |
382 | h->GetYaxis()->SetLabelFont(132); | |
383 | h->GetYaxis()->SetNdivisions(10); | |
384 | h->GetYaxis()->SetDecimals(); | |
385 | h->SetStats(0); | |
386 | } | |
387 | ||
4fa8d795 | 388 | //________________________________________________________________________ |
389 | void | |
390 | AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm) | |
391 | { | |
392 | // Normalize to the acceptance - | |
393 | // dndeta->Divide(accNorm); | |
394 | for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { | |
395 | Double_t a = norm->GetBinContent(i); | |
396 | for (Int_t j = 1; j <= copy->GetNbinsY(); j++) { | |
397 | if (a <= 0) { | |
398 | copy->SetBinContent(i,j,0); | |
399 | copy->SetBinError(i,j,0); | |
400 | continue; | |
401 | } | |
402 | Double_t c = copy->GetBinContent(i, j); | |
403 | Double_t e = copy->GetBinError(i, j); | |
404 | copy->SetBinContent(i, j, c / a); | |
405 | copy->SetBinError(i, j, e / a); | |
406 | } | |
407 | } | |
408 | } | |
5ca83fee | 409 | //________________________________________________________________________ |
410 | void | |
411 | AliBasedNdetaTask::ScaleToCoverage(TH1D* copy, const TH1D* norm) | |
412 | { | |
413 | // Normalize to the acceptance - | |
414 | // dndeta->Divide(accNorm); | |
415 | for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { | |
416 | Double_t a = norm->GetBinContent(i); | |
417 | if (a <= 0) { | |
418 | copy->SetBinContent(i,0); | |
419 | copy->SetBinError(i,0); | |
420 | continue; | |
421 | } | |
422 | Double_t c = copy->GetBinContent(i); | |
423 | Double_t e = copy->GetBinError(i); | |
424 | copy->SetBinContent(i, c / a); | |
425 | copy->SetBinError(i, e / a); | |
426 | } | |
427 | } | |
4fa8d795 | 428 | |
ce85db45 | 429 | //________________________________________________________________________ |
430 | TH1D* | |
431 | AliBasedNdetaTask::ProjectX(const TH2D* h, | |
432 | const char* name, | |
433 | Int_t firstbin, | |
434 | Int_t lastbin, | |
c25b5e1b | 435 | bool useRoot, |
ce85db45 | 436 | bool corr, |
e1f47419 | 437 | bool error) |
ce85db45 | 438 | { |
fb3430ac | 439 | // |
440 | // Project onto the X axis | |
441 | // | |
442 | // Parameters: | |
443 | // h 2D histogram | |
444 | // name New name | |
445 | // firstbin First bin to use | |
446 | // lastbin Last bin to use | |
447 | // error Whether to calculate errors | |
448 | // | |
449 | // Return: | |
450 | // Newly created histogram or null | |
451 | // | |
ce85db45 | 452 | if (!h) return 0; |
c25b5e1b | 453 | if (useRoot) |
454 | return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : "")); | |
ce85db45 | 455 | |
456 | TAxis* xaxis = h->GetXaxis(); | |
457 | TAxis* yaxis = h->GetYaxis(); | |
458 | TH1D* ret = new TH1D(name, h->GetTitle(), xaxis->GetNbins(), | |
459 | xaxis->GetXmin(), xaxis->GetXmax()); | |
460 | static_cast<const TAttLine*>(h)->Copy(*ret); | |
461 | static_cast<const TAttFill*>(h)->Copy(*ret); | |
462 | static_cast<const TAttMarker*>(h)->Copy(*ret); | |
463 | ret->GetXaxis()->ImportAttributes(xaxis); | |
464 | ||
465 | Int_t first = firstbin; | |
466 | Int_t last = lastbin; | |
5ca83fee | 467 | if (first < 0) first = 1; |
468 | else if (first >= yaxis->GetNbins()+2) first = yaxis->GetNbins()+1; | |
ce85db45 | 469 | if (last < 0) last = yaxis->GetNbins(); |
5ca83fee | 470 | else if (last >= yaxis->GetNbins()+2) last = yaxis->GetNbins()+1; |
ce85db45 | 471 | if (last-first < 0) { |
e1f47419 | 472 | AliWarningGeneral("AliBasedNdetaTask", |
473 | Form("Nothing to project [%d,%d]", first, last)); | |
ce85db45 | 474 | return 0; |
475 | ||
476 | } | |
477 | ||
478 | // Loop over X bins | |
c8b1a7db | 479 | //DMSG(fDebug,3,"Projecting bins [%d,%d] of %s", first, last, h->GetName())); |
ce85db45 | 480 | Int_t ybins = (last-first+1); |
481 | for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) { | |
482 | Double_t content = 0; | |
483 | Double_t error2 = 0; | |
484 | Int_t nbins = 0; | |
485 | ||
486 | ||
487 | for (Int_t ybin = first; ybin <= last; ybin++) { | |
b8f92f9d | 488 | Double_t c1 = h->GetBinContent(h->GetBin(xbin, ybin)); |
489 | Double_t e1 = h->GetBinError(h->GetBin(xbin, ybin)); | |
ce85db45 | 490 | |
491 | // Ignore empty bins | |
492 | if (c1 < 1e-12) continue; | |
493 | if (e1 < 1e-12) { | |
494 | if (error) continue; | |
495 | e1 = 1; | |
496 | } | |
497 | ||
498 | content += c1; | |
499 | error2 += e1*e1; | |
500 | nbins++; | |
501 | } // for (ybin) | |
502 | if(content > 0 && nbins > 0) { | |
503 | Double_t factor = (corr ? Double_t(ybins) / nbins : 1); | |
c25b5e1b | 504 | #if 0 |
505 | AliWarningGeneral(ret->GetName(), | |
506 | Form("factor @ %d is %d/%d -> %f", | |
507 | xbin, ybins, nbins, factor)); | |
508 | #endif | |
ce85db45 | 509 | if (error) { |
510 | // Calculate weighted average | |
511 | ret->SetBinContent(xbin, content * factor); | |
512 | ret->SetBinError(xbin, factor * TMath::Sqrt(error2)); | |
513 | } | |
514 | else | |
515 | ret->SetBinContent(xbin, factor * content); | |
516 | } | |
517 | } // for (xbin) | |
518 | ||
519 | return ret; | |
520 | } | |
bfab35d9 | 521 | |
ce85db45 | 522 | //________________________________________________________________________ |
c8b1a7db | 523 | Bool_t |
524 | AliBasedNdetaTask::Finalize() | |
ce85db45 | 525 | { |
fb3430ac | 526 | // |
527 | // Called at end of event processing.. | |
528 | // | |
529 | // This is called once in the master | |
530 | // | |
531 | // Parameters: | |
532 | // option Not used | |
533 | // | |
ce85db45 | 534 | // Draw result to screen, or perform fitting, normalizations Called |
535 | // once at the end of the query | |
f53fb4f6 | 536 | DGUARD(fDebug,1,"Process final merged results"); |
0be6c8cd | 537 | |
c8b1a7db | 538 | UShort_t sNN; |
539 | UShort_t sys; | |
540 | ULong_t trig; | |
541 | UShort_t scheme; | |
542 | AliForwardUtil::GetParameter(fSums->FindObject("sNN"), sNN); | |
543 | AliForwardUtil::GetParameter(fSums->FindObject("sys"), sys); | |
544 | AliForwardUtil::GetParameter(fSums->FindObject("trigger"), trig); | |
545 | AliForwardUtil::GetParameter(fSums->FindObject("scheme"), scheme); | |
546 | ||
547 | TAxis* cA = static_cast<TAxis*>(fSums->FindObject("centAxis")); | |
548 | if (cA) cA->Copy(fCentAxis); | |
ce85db45 | 549 | |
c8b1a7db | 550 | if(sNN > 0 && sys == AliForwardUtil::kPP) |
551 | LoadNormalizationData(sys, sNN); | |
f67d699c | 552 | |
553 | InitializeCentBins(); | |
554 | ||
0be6c8cd | 555 | // Print before we loop |
556 | Print(); | |
ce85db45 | 557 | |
e1f47419 | 558 | // Loop over centrality bins |
559 | TIter next(fListOfCentralities); | |
560 | CentralityBin* bin = 0; | |
5bb5d1f6 | 561 | gStyle->SetPalette(1); |
562 | THStack* dndetaStack = new THStack("dndeta", "dN/d#eta"); | |
563 | THStack* dndetaStackRebin = new THStack(Form("dndeta_rebin%02d", fRebin), | |
564 | "dN_{ch}/d#eta"); | |
565 | THStack* dndetaMCStack = new THStack("dndetaMC", "dN_{ch}/d#eta"); | |
566 | THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin), | |
567 | "dN_{ch}/d#eta"); | |
568 | ||
c89b9ac1 | 569 | TList* mclist = 0; |
570 | TList* truthlist = 0; | |
571 | ||
797161e8 | 572 | if (fFinalMCCorrFile.Contains(".root")) { |
c89b9ac1 | 573 | TFile* ftest = TFile::Open(fFinalMCCorrFile.Data()); |
574 | if(ftest) { | |
797161e8 | 575 | mclist = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName()))); |
576 | truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults")); | |
c89b9ac1 | 577 | } |
797161e8 | 578 | else |
579 | AliWarning("MC analysis file invalid - no final MC correction possible"); | |
c89b9ac1 | 580 | } |
9ecab72f | 581 | Int_t style = GetMarker(); |
582 | Int_t color = GetColor(); | |
c89b9ac1 | 583 | |
c8b1a7db | 584 | DMSG(fDebug,3,"Marker style=%d, color=%d", style, color); |
5bb5d1f6 | 585 | while ((bin = static_cast<CentralityBin*>(next()))) { |
c8b1a7db | 586 | bin->End(fSums, fResults, fNormalizationScheme, fShapeCorr, |
66cf95f2 | 587 | fTriggerEff, fTriggerEff0, |
c25b5e1b | 588 | fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges, |
c89b9ac1 | 589 | fTriggerMask, style, color, mclist, truthlist); |
c8b1a7db | 590 | if (HasCentrality() && bin->IsAllBin()) |
591 | // If we have centrality bins, do not add the min-bias | |
592 | // distribution to the output stack. | |
593 | continue; | |
594 | TH1* dndeta = bin->GetResult(0, false, ""); | |
595 | TH1* dndetaSym = fSymmetrice ? bin->GetResult(0, true, "") : 0; | |
596 | TH1* dndetaMC = bin->GetResult(0, false, "MC", false); | |
597 | TH1* dndetaMCSym = fSymmetrice ? bin->GetResult(0, true, "MC", false) : 0; | |
598 | DMSG(fDebug,2,"Results: bare=%p sym=%p mcbare=%p mcsym=%p", | |
599 | dndeta, dndetaSym, dndetaMC, dndetaMCSym); | |
5bb5d1f6 | 600 | if (dndeta) dndetaStack->Add(dndeta); |
601 | if (dndetaSym) dndetaStack->Add(dndetaSym); | |
602 | if (dndetaMC) dndetaMCStack->Add(dndetaMC); | |
603 | if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym); | |
604 | if (fRebin > 1) { | |
c8b1a7db | 605 | dndeta = bin->GetResult(fRebin, false, ""); |
606 | dndetaSym = fSymmetrice ? bin->GetResult(fRebin, true, "") : 0; | |
607 | dndetaMC = bin->GetResult(fRebin, false, "MC", false); | |
608 | dndetaMCSym = fSymmetrice ? bin->GetResult(fRebin, true, "MC", false): 0; | |
5bb5d1f6 | 609 | if (dndeta) dndetaStackRebin->Add(dndeta); |
610 | if (dndetaSym) dndetaStackRebin->Add(dndetaSym); | |
611 | if (dndetaMC) dndetaMCStackRebin->Add(dndetaMC); | |
612 | if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym); | |
613 | } | |
614 | } | |
615 | // Output the stack | |
c8b1a7db | 616 | fResults->Add(dndetaStack); |
5bb5d1f6 | 617 | |
618 | // If available output rebinned stack | |
c8b1a7db | 619 | if (fRebin > 0 && (!dndetaStackRebin->GetHists() || |
620 | dndetaStackRebin->GetHists()->GetEntries() <= 0)) { | |
5bb5d1f6 | 621 | AliWarning("No rebinned histograms found"); |
622 | delete dndetaStackRebin; | |
623 | dndetaStackRebin = 0; | |
624 | } | |
c8b1a7db | 625 | if (dndetaStackRebin) fResults->Add(dndetaStackRebin); |
5bb5d1f6 | 626 | |
627 | // If available, output track-ref stack | |
628 | if (!dndetaMCStack->GetHists() || | |
629 | dndetaMCStack->GetHists()->GetEntries() <= 0) { | |
c8b1a7db | 630 | // AliWarning("No MC histograms found"); |
5bb5d1f6 | 631 | delete dndetaMCStack; |
632 | dndetaMCStack = 0; | |
633 | } | |
c8b1a7db | 634 | if (dndetaMCStack) fResults->Add(dndetaMCStack); |
5bb5d1f6 | 635 | |
636 | // If available, output rebinned track-ref stack | |
c8b1a7db | 637 | if ((fRebin > 0 && dndetaMCStack) |
638 | && (!dndetaMCStackRebin->GetHists() || | |
639 | dndetaMCStackRebin->GetHists()->GetEntries() <= 0)) { | |
5bb5d1f6 | 640 | AliWarning("No rebinned MC histograms found"); |
641 | delete dndetaMCStackRebin; | |
642 | dndetaMCStackRebin = 0; | |
643 | } | |
c8b1a7db | 644 | if (dndetaMCStackRebin) fResults->Add(dndetaMCStackRebin); |
ce85db45 | 645 | |
e1f47419 | 646 | // Output collision energy string |
c8b1a7db | 647 | if (sNN > 0) { |
648 | TNamed* sNNObj = new TNamed("sNN", | |
a76fb27d | 649 | AliForwardUtil::CenterOfMassEnergyString(sNN)); |
650 | sNNObj->SetUniqueID(sNN); | |
c8b1a7db | 651 | fResults->Add(sNNObj); // fSNNString->Clone()); |
a76fb27d | 652 | } |
ce85db45 | 653 | |
e1f47419 | 654 | // Output collision system string |
c8b1a7db | 655 | if (sys > 0) { |
656 | TNamed* sysObj = new TNamed("sys", | |
a76fb27d | 657 | AliForwardUtil::CollisionSystemString(sys)); |
658 | sysObj->SetUniqueID(sys); | |
c8b1a7db | 659 | fResults->Add(sysObj); // fSysString->Clone()); |
a76fb27d | 660 | } |
ce85db45 | 661 | |
e308a636 | 662 | // Output centrality axis |
c8b1a7db | 663 | fResults->Add(&fCentAxis); |
e308a636 | 664 | |
e1f47419 | 665 | // Output trigger string |
c8b1a7db | 666 | if (trig) { |
667 | TNamed* maskObj = new TNamed("trigger", | |
668 | AliAODForwardMult::GetTriggerString(trig)); | |
669 | maskObj->SetUniqueID(trig); | |
670 | fResults->Add(maskObj); // fTriggerString->Clone()); | |
a76fb27d | 671 | } |
0be6c8cd | 672 | |
673 | // Normalization string | |
c8b1a7db | 674 | if (scheme > 0) { |
675 | TNamed* schemeObj = new TNamed("scheme", | |
a76fb27d | 676 | NormalizationSchemeString(scheme)); |
677 | schemeObj->SetUniqueID(scheme); | |
c8b1a7db | 678 | fResults->Add(schemeObj); // fSchemeString->Clone()); |
a76fb27d | 679 | } |
ce85db45 | 680 | |
e1f47419 | 681 | // Output vertex axis |
c8b1a7db | 682 | TAxis* vtxAxis = new TAxis(1,fMinIpZ,fMaxIpZ); |
ce85db45 | 683 | vtxAxis->SetName("vtxAxis"); |
c8b1a7db | 684 | vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fMinIpZ,fMaxIpZ)); |
685 | fResults->Add(vtxAxis); | |
e1f47419 | 686 | |
bfab35d9 | 687 | // Output trigger efficiency and shape correction |
c8b1a7db | 688 | fResults->Add(AliForwardUtil::MakeParameter("triggerEff", fTriggerEff)); |
689 | fResults->Add(AliForwardUtil::MakeParameter("triggerEff0", fTriggerEff0)); | |
690 | if (fShapeCorr) fResults->Add(fShapeCorr); | |
ffca499d | 691 | |
c25b5e1b | 692 | TNamed* options = new TNamed("options",""); |
693 | TString str; | |
694 | str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not ")); | |
695 | str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not ")); | |
696 | str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not ")); | |
697 | options->SetTitle(str); | |
c8b1a7db | 698 | fResults->Add(options); |
c25b5e1b | 699 | |
c8b1a7db | 700 | return true; |
ce85db45 | 701 | } |
3846ec25 | 702 | //________________________________________________________________________ |
703 | void | |
704 | AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy) | |
705 | { | |
e1f47419 | 706 | // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD |
f53fb4f6 | 707 | DGUARD(fDebug,1,"Load the normalization data for sys=%d, energy=%d", |
708 | sys, energy); | |
3846ec25 | 709 | TString type("pp"); |
3846ec25 | 710 | TString snn("900"); |
711 | if(energy == 7000) snn.Form("7000"); | |
712 | if(energy == 2750) snn.Form("2750"); | |
713 | ||
d43c6cd0 | 714 | // Check if shape correction/trigger efficiency was requsted and not |
715 | // already set | |
716 | Bool_t needShape = ((fNormalizationScheme & kShape) && !fShapeCorr); | |
717 | Bool_t needEff = ((fNormalizationScheme & kTriggerEfficiency) && | |
718 | ((1 - fTriggerEff) < 1e-6) && fTriggerEff > 0); | |
c8b1a7db | 719 | if (needShape) DMSG(fDebug, 0, "Will load shape correction"); |
720 | if (needEff) DMSG(fDebug, 0, "Will load trigger efficiency, was=%f, %f", | |
721 | fTriggerEff, fTriggerEff0); | |
722 | if(!needShape) { // && !needEff) { | |
723 | DMSG(fDebug,1,"Objects already set for normalization - no action taken"); | |
e1f47419 | 724 | return; |
725 | } | |
d43c6cd0 | 726 | |
727 | TString fname(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/" | |
728 | "Normalization/normalizationHists_%s_%s.root", | |
729 | type.Data(),snn.Data())); | |
730 | AliWarningF("Using old-style corrections from %s", fname.Data()); | |
731 | TFile* fin = TFile::Open(fname, "READ"); | |
e1f47419 | 732 | if(!fin) { |
d43c6cd0 | 733 | AliWarningF("no file for normalization of %d/%d (%s)", |
734 | sys, energy, fname.Data()); | |
e1f47419 | 735 | return; |
736 | } | |
3846ec25 | 737 | |
4fa8d795 | 738 | // Shape correction |
d43c6cd0 | 739 | if (needShape) { |
4fa8d795 | 740 | TString trigName("All"); |
741 | if (fTriggerMask == AliAODForwardMult::kInel || | |
742 | fTriggerMask == AliAODForwardMult::kNClusterGt0) | |
743 | trigName = "Inel"; | |
744 | else if (fTriggerMask == AliAODForwardMult::kNSD) | |
745 | trigName = "NSD"; | |
746 | else if (fTriggerMask == AliAODForwardMult::kInelGt0) | |
747 | trigName = "InelGt0"; | |
748 | else { | |
749 | AliWarning(Form("Normalization for trigger %s not known, using all", | |
750 | AliAODForwardMult::GetTriggerString(fTriggerMask))); | |
751 | } | |
752 | ||
753 | TString shapeCorName(Form("h%sNormalization", trigName.Data())); | |
754 | TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName)); | |
755 | if (shapeCor) SetShapeCorrection(shapeCor); | |
756 | else { | |
757 | AliWarning(Form("No shape correction found for %s", trigName.Data())); | |
758 | } | |
759 | } | |
ce85db45 | 760 | |
e1f47419 | 761 | // Trigger efficiency |
d43c6cd0 | 762 | if (needEff) { |
763 | TString effName(Form("%sTriggerEff", | |
764 | fTriggerMask == AliAODForwardMult::kInel ? "inel" : | |
765 | fTriggerMask == AliAODForwardMult::kNSD ? "nsd" : | |
766 | fTriggerMask == AliAODForwardMult::kInelGt0 ? | |
767 | "inelgt0" : "all")); | |
768 | Double_t trigEff = 1; | |
241cca4d | 769 | TObject* eff = fin->Get(effName); |
770 | if (eff) AliForwardUtil::GetParameter(eff, trigEff); | |
3846ec25 | 771 | |
d43c6cd0 | 772 | if (trigEff <= 0) |
773 | AliWarningF("Retrieved trigger efficiency %s is %f<=0, ignoring", | |
774 | effName.Data(), trigEff); | |
775 | else | |
776 | SetTriggerEff(trigEff); | |
777 | ||
778 | // Trigger efficiency | |
779 | TString eff0Name(effName); | |
780 | eff0Name.Append("0"); | |
781 | ||
782 | Double_t trigEff0 = 1; | |
783 | TObject* eff0 = fin->Get(eff0Name); | |
784 | if (eff0) AliForwardUtil::GetParameter(eff, trigEff0); | |
785 | if (trigEff0 < 0) | |
786 | AliWarningF("Retrieved trigger efficiency %s is %f<0, ignoring", | |
787 | eff0Name.Data(), trigEff0); | |
788 | else | |
789 | SetTriggerEff0(trigEff0); | |
66cf95f2 | 790 | } |
d43c6cd0 | 791 | |
b30dee70 | 792 | // TEMPORARY FIX |
793 | // Rescale the shape correction by the trigger efficiency | |
4fa8d795 | 794 | if (fShapeCorr) { |
795 | AliWarning(Form("Rescaling shape correction by trigger efficiency: " | |
d43c6cd0 | 796 | "1/E_X=1/%f", fTriggerEff)); |
797 | fShapeCorr->Scale(1. / fTriggerEff); | |
4fa8d795 | 798 | } |
d43c6cd0 | 799 | if (fin) fin->Close(); |
b30dee70 | 800 | |
e1f47419 | 801 | // Print - out |
c8b1a7db | 802 | if (fDebug > 1 && fShapeCorr && fTriggerEff) |
803 | DMSG(fDebug, 1, "Loaded objects for normalization."); | |
3846ec25 | 804 | } |
0be6c8cd | 805 | |
806 | ||
c8b1a7db | 807 | #define PF(N,V,...) \ |
808 | AliForwardUtil::PrintField(N,V, ## __VA_ARGS__) | |
809 | #define PFB(N,FLAG) \ | |
810 | do { \ | |
811 | AliForwardUtil::PrintName(N); \ | |
812 | std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \ | |
813 | } while(false) | |
814 | #define PFV(N,VALUE) \ | |
815 | do { \ | |
816 | AliForwardUtil::PrintName(N); \ | |
817 | std::cout << (VALUE) << std::endl; } while(false) | |
818 | ||
0be6c8cd | 819 | //________________________________________________________________________ |
820 | void | |
c8b1a7db | 821 | AliBasedNdetaTask::Print(Option_t* option) const |
0be6c8cd | 822 | { |
823 | // | |
824 | // Print information | |
825 | // | |
c8b1a7db | 826 | AliBaseAODTask::Print(option); |
827 | TString schemeString(NormalizationSchemeString(fNormalizationScheme)); | |
828 | ||
829 | gROOT->IncreaseDirLevel(); | |
830 | PFV("Rebin factor", fRebin); | |
831 | PFB("Cut edges", fCutEdges); | |
832 | PFB("Symmertrice", fSymmetrice); | |
833 | PFB("Use TH2::ProjectionX", fUseROOTProj); | |
834 | PFB("Correct for empty", fCorrEmpty); | |
835 | PFV("Normalization scheme", schemeString ); | |
836 | PFV("Trigger efficiency", fTriggerEff); | |
837 | PFV("Bin-0 Trigger efficiency", fTriggerEff0); | |
838 | PFV("Shape correction", (fShapeCorr?fShapeCorr->GetName():"none"));; | |
839 | gROOT->DecreaseDirLevel(); | |
0be6c8cd | 840 | } |
841 | ||
ce85db45 | 842 | //________________________________________________________________________ |
843 | TH1D* | |
e1f47419 | 844 | AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges) |
ce85db45 | 845 | { |
fb3430ac | 846 | // |
847 | // Make a copy of the input histogram and rebin that histogram | |
848 | // | |
849 | // Parameters: | |
850 | // h Histogram to rebin | |
851 | // | |
852 | // Return: | |
853 | // New (rebinned) histogram | |
854 | // | |
e1f47419 | 855 | if (rebin <= 1) return 0; |
ce85db45 | 856 | |
857 | Int_t nBins = h->GetNbinsX(); | |
e1f47419 | 858 | if(nBins % rebin != 0) { |
859 | AliWarningGeneral("AliBasedNdetaTask", | |
860 | Form("Rebin factor %d is not a devisor of current number " | |
861 | "of bins %d in the histogram %s", | |
862 | rebin, nBins, h->GetName())); | |
ce85db45 | 863 | return 0; |
864 | } | |
865 | ||
866 | // Make a copy | |
867 | TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d", | |
e1f47419 | 868 | h->GetName(), rebin))); |
869 | tmp->Rebin(rebin); | |
bfab35d9 | 870 | // Hist should be reset, as it otherwise messes with the cutEdges option |
1f7aa5c7 | 871 | tmp->Reset(); |
ce85db45 | 872 | tmp->SetDirectory(0); |
873 | ||
874 | // The new number of bins | |
e1f47419 | 875 | Int_t nBinsNew = nBins / rebin; |
ce85db45 | 876 | for(Int_t i = 1;i<= nBinsNew; i++) { |
877 | Double_t content = 0; | |
878 | Double_t sumw = 0; | |
879 | Double_t wsum = 0; | |
880 | Int_t nbins = 0; | |
e1f47419 | 881 | for(Int_t j = 1; j<=rebin;j++) { |
882 | Int_t bin = (i-1)*rebin + j; | |
ce85db45 | 883 | Double_t c = h->GetBinContent(bin); |
bfab35d9 | 884 | if (c <= 0) { |
885 | continue; // old TODO: check | |
886 | //content = -1; // new | |
887 | //break; // also new | |
888 | } | |
889 | ||
e1f47419 | 890 | if (cutEdges) { |
ce85db45 | 891 | if (h->GetBinContent(bin+1)<=0 || |
3846ec25 | 892 | h->GetBinContent(bin-1)<=0) { |
e1f47419 | 893 | #if 0 |
894 | AliWarningGeneral("AliBasedNdetaTask", | |
895 | Form("removing bin %d=%f of %s (%d=%f,%d=%f)", | |
896 | bin, c, h->GetName(), | |
897 | bin+1, h->GetBinContent(bin+1), | |
898 | bin-1, h->GetBinContent(bin-1))); | |
899 | #endif | |
ce85db45 | 900 | continue; |
901 | } | |
902 | } | |
903 | Double_t e = h->GetBinError(bin); | |
904 | Double_t w = 1 / (e*e); // 1/c/c | |
905 | content += c; | |
906 | sumw += w; | |
907 | wsum += w * c; | |
908 | nbins++; | |
909 | } | |
910 | ||
911 | if(content > 0 && nbins > 0) { | |
912 | tmp->SetBinContent(i, wsum / sumw); | |
913 | tmp->SetBinError(i,1./TMath::Sqrt(sumw)); | |
914 | } | |
915 | } | |
916 | ||
917 | return tmp; | |
918 | } | |
919 | ||
920 | //__________________________________________________________________ | |
ce85db45 | 921 | TH1* |
e1f47419 | 922 | AliBasedNdetaTask::Symmetrice(const TH1* h) |
ce85db45 | 923 | { |
fb3430ac | 924 | // |
925 | // Make an extension of @a h to make it symmetric about 0 | |
926 | // | |
927 | // Parameters: | |
928 | // h Histogram to symmertrice | |
929 | // | |
930 | // Return: | |
931 | // Symmetric extension of @a h | |
932 | // | |
ce85db45 | 933 | Int_t nBins = h->GetNbinsX(); |
934 | TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName()))); | |
935 | s->SetTitle(Form("%s (mirrored)", h->GetTitle())); | |
936 | s->Reset(); | |
937 | s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin()); | |
938 | s->SetMarkerColor(h->GetMarkerColor()); | |
939 | s->SetMarkerSize(h->GetMarkerSize()); | |
9ecab72f | 940 | s->SetMarkerStyle(FlipHollowStyle(h->GetMarkerStyle())); |
ce85db45 | 941 | s->SetFillColor(h->GetFillColor()); |
942 | s->SetFillStyle(h->GetFillStyle()); | |
943 | s->SetDirectory(0); | |
944 | ||
945 | // Find the first and last bin with data | |
946 | Int_t first = nBins+1; | |
947 | Int_t last = 0; | |
948 | for (Int_t i = 1; i <= nBins; i++) { | |
949 | if (h->GetBinContent(i) <= 0) continue; | |
950 | first = TMath::Min(first, i); | |
951 | last = TMath::Max(last, i); | |
952 | } | |
953 | ||
954 | Double_t xfirst = h->GetBinCenter(first-1); | |
955 | Int_t f1 = h->GetXaxis()->FindBin(-xfirst); | |
956 | Int_t l2 = s->GetXaxis()->FindBin(xfirst); | |
957 | for (Int_t i = f1, j=l2; i <= last; i++,j--) { | |
958 | s->SetBinContent(j, h->GetBinContent(i)); | |
959 | s->SetBinError(j, h->GetBinError(i)); | |
960 | } | |
961 | // Fill in overlap bin | |
962 | s->SetBinContent(l2+1, h->GetBinContent(first)); | |
963 | s->SetBinError(l2+1, h->GetBinError(first)); | |
964 | return s; | |
965 | } | |
e1f47419 | 966 | |
9ecab72f | 967 | //__________________________________________________________________ |
968 | Int_t | |
969 | AliBasedNdetaTask::GetMarkerStyle(UShort_t bits) | |
970 | { | |
971 | Int_t base = bits & (0xFE); | |
972 | Bool_t hollow = bits & kHollow; | |
973 | switch (base) { | |
974 | case kCircle: return (hollow ? 24 : 20); | |
975 | case kSquare: return (hollow ? 25 : 21); | |
976 | case kUpTriangle: return (hollow ? 26 : 22); | |
977 | case kDownTriangle: return (hollow ? 32 : 23); | |
978 | case kDiamond: return (hollow ? 27 : 33); | |
979 | case kCross: return (hollow ? 28 : 34); | |
980 | case kStar: return (hollow ? 30 : 29); | |
981 | } | |
982 | return 1; | |
983 | } | |
984 | //__________________________________________________________________ | |
985 | UShort_t | |
986 | AliBasedNdetaTask::GetMarkerBits(Int_t style) | |
987 | { | |
988 | UShort_t bits = 0; | |
989 | switch (style) { | |
990 | case 24: case 25: case 26: case 27: case 28: case 30: case 32: | |
991 | bits |= kHollow; break; | |
992 | } | |
993 | switch (style) { | |
994 | case 20: case 24: bits |= kCircle; break; | |
995 | case 21: case 25: bits |= kSquare; break; | |
996 | case 22: case 26: bits |= kUpTriangle; break; | |
997 | case 23: case 32: bits |= kDownTriangle; break; | |
998 | case 27: case 33: bits |= kDiamond; break; | |
999 | case 28: case 34: bits |= kCross; break; | |
1000 | case 29: case 30: bits |= kStar; break; | |
1001 | } | |
1002 | return bits; | |
1003 | } | |
1004 | //__________________________________________________________________ | |
1005 | Int_t | |
1006 | AliBasedNdetaTask::FlipHollowStyle(Int_t style) | |
1007 | { | |
1008 | UShort_t bits = GetMarkerBits(style); | |
1009 | Int_t ret = GetMarkerStyle(bits ^ kHollow); | |
1010 | return ret; | |
1011 | } | |
1012 | ||
4fa8d795 | 1013 | //==================================================================== |
1014 | void | |
1015 | AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col) | |
1016 | { | |
f53fb4f6 | 1017 | DGUARD(fDebug,1,"Initializing sums with %s", data->GetName()); |
4fa8d795 | 1018 | TString n(GetHistName(0)); |
1019 | TString n0(GetHistName(1)); | |
1020 | const char* postfix = GetTitle(); | |
1021 | ||
1022 | fSum = static_cast<TH2D*>(data->Clone(n)); | |
1023 | if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix)); | |
1024 | fSum->SetDirectory(0); | |
1025 | fSum->SetMarkerColor(col); | |
9ecab72f | 1026 | fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid)); |
4fa8d795 | 1027 | fSum->Reset(); |
1028 | list->Add(fSum); | |
1029 | ||
1030 | fSum0 = static_cast<TH2D*>(data->Clone(n0)); | |
1031 | if (postfix) | |
1032 | fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix)); | |
1033 | else | |
1034 | fSum0->SetTitle(Form("%s 0-bin", data->GetTitle())); | |
1035 | fSum0->SetDirectory(0); | |
1036 | fSum0->SetMarkerColor(col); | |
9ecab72f | 1037 | fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow)); |
4fa8d795 | 1038 | fSum0->Reset(); |
1039 | list->Add(fSum0); | |
1040 | ||
1041 | fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5); | |
1042 | fEvents->SetDirectory(0); | |
1043 | fEvents->GetXaxis()->SetBinLabel(1, "Non-zero"); | |
1044 | fEvents->GetXaxis()->SetBinLabel(2, "Zero"); | |
1045 | list->Add(fEvents); | |
1046 | } | |
1047 | ||
1048 | //____________________________________________________________________ | |
1049 | TString | |
c8b1a7db | 1050 | AliBasedNdetaTask::Sum::GetHistName(const char* /*name*/, |
1051 | Int_t what, const char* post) | |
4fa8d795 | 1052 | { |
c8b1a7db | 1053 | TString n; |
1054 | switch (what) { | |
1055 | case 0: n = "sum"; break; | |
1056 | case 1: n = "sum0"; break; | |
1057 | case 2: n = "events"; break; | |
1058 | } | |
f67d699c | 1059 | if (post && post[0] != '\0') n.Append(post); |
4fa8d795 | 1060 | return n; |
1061 | } | |
1062 | ||
f67d699c | 1063 | //____________________________________________________________________ |
1064 | TString | |
1065 | AliBasedNdetaTask::Sum::GetHistName(Int_t what) const | |
1066 | { | |
1067 | return GetHistName(GetName(), what, GetTitle()); | |
1068 | } | |
1069 | ||
4fa8d795 | 1070 | //____________________________________________________________________ |
1071 | void | |
1072 | AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero) | |
1073 | { | |
f53fb4f6 | 1074 | DGUARD(fDebug,2,"Adding %s to sums", data->GetName()); |
4fa8d795 | 1075 | if (isZero) fSum0->Add(data); |
1076 | else fSum->Add(data); | |
1077 | fEvents->Fill(isZero ? 1 : 0); | |
1078 | } | |
1079 | ||
1080 | //____________________________________________________________________ | |
1081 | TH2D* | |
f67d699c | 1082 | AliBasedNdetaTask::Sum::CalcSum(TList* output, |
1083 | Double_t& ntotal, | |
1084 | Double_t epsilon0, | |
1085 | Double_t epsilon, | |
1086 | Int_t marker, | |
1087 | Bool_t rootProj, | |
1088 | Bool_t corrEmpty) const | |
4fa8d795 | 1089 | { |
f67d699c | 1090 | DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName()); |
9ecab72f | 1091 | |
e12f8c42 | 1092 | // The return value `ret' is not scaled in anyway |
f67d699c | 1093 | TH2D* ret = static_cast<TH2D*>(fSum->Clone(fSum->GetName())); |
4fa8d795 | 1094 | ret->SetDirectory(0); |
f67d699c | 1095 | Int_t n = Int_t(fEvents->GetBinContent(1)); |
1096 | Int_t n0 = Int_t(fEvents->GetBinContent(2)); | |
c8b1a7db | 1097 | ntotal = n; |
1098 | if (n0 > 0) { | |
1099 | ret->Reset(); | |
1100 | DMSG(fDebug,1, | |
1101 | "Adding histograms %s(%d) and %s(%d) w/weights %f and %f resp.", | |
1102 | fSum0->GetName(), n, fSum->GetName(), n0, 1./epsilon,1./epsilon0); | |
1103 | ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon); | |
1104 | ntotal = n / epsilon + n0 / epsilon0; | |
1105 | } | |
4fa8d795 | 1106 | |
1107 | TList* out = new TList; | |
1108 | out->SetOwner(); | |
1109 | const char* postfix = GetTitle(); | |
1110 | if (!postfix) postfix = ""; | |
1111 | out->SetName(Form("partial%s", postfix)); | |
1112 | output->Add(out); | |
1113 | ||
1114 | // Now make copies, normalize them, and store in output list | |
e12f8c42 | 1115 | // Note, these are the only ones normalized here |
1116 | // These are mainly for diagnostics | |
f67d699c | 1117 | TH2D* sumCopy = static_cast<TH2D*>(fSum->Clone("sum")); |
1118 | TH2D* sum0Copy = static_cast<TH2D*>(fSum0->Clone("sum0")); | |
4fa8d795 | 1119 | TH2D* retCopy = static_cast<TH2D*>(ret->Clone("sumAll")); |
9ecab72f | 1120 | sumCopy->SetMarkerStyle(FlipHollowStyle(marker)); |
4fa8d795 | 1121 | sumCopy->SetDirectory(0); |
9ecab72f | 1122 | sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4)); |
4fa8d795 | 1123 | sum0Copy->SetDirectory(0); |
1124 | retCopy->SetMarkerStyle(marker); | |
1125 | retCopy->SetDirectory(0); | |
1126 | ||
5ca83fee | 1127 | Int_t nY = fSum->GetNbinsY(); |
1128 | Int_t o = 0; // nY+1; | |
1129 | TH1D* norm = ProjectX(fSum, "norm", o, o, rootProj, corrEmpty, false); | |
1130 | TH1D* norm0 = ProjectX(fSum0, "norm0", o, o, rootProj, corrEmpty, false); | |
1131 | TH1D* normAll = ProjectX(ret, "normAll", o, o, rootProj, corrEmpty, false); | |
e12f8c42 | 1132 | norm->SetTitle("#eta coverage - >0-bin"); |
1133 | norm0->SetTitle("#eta coverage - 0-bin"); | |
1134 | normAll->SetTitle("#eta coverage"); | |
4fa8d795 | 1135 | norm->SetDirectory(0); |
1136 | norm0->SetDirectory(0); | |
1137 | normAll->SetDirectory(0); | |
1138 | ||
c8b1a7db | 1139 | TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1,nY,rootProj,corrEmpty); |
1140 | TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1,nY,rootProj,corrEmpty); | |
1141 | TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1,nY,rootProj,corrEmpty); | |
1142 | sumCopyPx-> SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sumCopy->GetTitle())); | |
e12f8c42 | 1143 | sum0CopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sum0Copy->GetTitle())); |
c8b1a7db | 1144 | retCopyPx-> SetTitle(Form("#sum_{i}^{N_{#phi}}%s", retCopy->GetTitle())); |
1145 | sumCopyPx-> SetDirectory(0); | |
4fa8d795 | 1146 | sum0CopyPx->SetDirectory(0); |
c8b1a7db | 1147 | retCopyPx-> SetDirectory(0); |
4fa8d795 | 1148 | |
c8b1a7db | 1149 | TH1D* phi = ProjectX(fSum, "phi", nY+1,nY+1,rootProj,corrEmpty,false); |
1150 | TH1D* phi0 = ProjectX(fSum0, "phi0", nY+1,nY+1,rootProj,corrEmpty,false); | |
1151 | TH1D* phiAll = ProjectX(ret, "phiAll", nY+1,nY+1,rootProj,corrEmpty,false); | |
1152 | phi ->SetTitle("#phi acceptance from dead strips - >0-bin"); | |
1153 | phi0 ->SetTitle("#phi acceptance from dead strips - 0-bin"); | |
e12f8c42 | 1154 | phiAll->SetTitle("#phi acceptance from dead strips"); |
c8b1a7db | 1155 | phi ->SetDirectory(0); |
1156 | phi0 ->SetDirectory(0); | |
5ca83fee | 1157 | phiAll->SetDirectory(0); |
1158 | ||
e12f8c42 | 1159 | const TH1D* cov = (corrEmpty ? norm : phi); |
1160 | const TH1D* cov0 = (corrEmpty ? norm0 : phi0); | |
1161 | const TH1D* covAll = (corrEmpty ? normAll : phiAll); | |
1162 | ||
1163 | // Here, we scale to the coverage (or phi acceptance) | |
1164 | ScaleToCoverage(sumCopy, cov); | |
1165 | ScaleToCoverage(sum0Copy, cov0); | |
1166 | ScaleToCoverage(retCopy, covAll); | |
1167 | ||
4fa8d795 | 1168 | // Scale our 1D histograms |
c8b1a7db | 1169 | sumCopyPx ->Scale(1., "width"); |
4fa8d795 | 1170 | sum0CopyPx->Scale(1., "width"); |
c8b1a7db | 1171 | retCopyPx ->Scale(1., "width"); |
4fa8d795 | 1172 | |
c8b1a7db | 1173 | DMSG(fDebug,2,"Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(), |
1174 | sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()); | |
4fa8d795 | 1175 | |
1176 | // Scale the normalization - they should be 1 at the maximum | |
c8b1a7db | 1177 | norm ->Scale(n > 0 ? 1. / n : 1); |
1178 | norm0 ->Scale(n0 > 0 ? 1. / n0 : 1); | |
4fa8d795 | 1179 | normAll->Scale(ntotal > 0 ? 1. / ntotal : 1); |
1180 | ||
e12f8c42 | 1181 | // Scale the normalization - they should be 1 at the maximum |
c8b1a7db | 1182 | phi ->Scale(n > 0 ? 1. / n : 1); |
1183 | phi0 ->Scale(n0 > 0 ? 1. / n0 : 1); | |
e12f8c42 | 1184 | phiAll->Scale(ntotal > 0 ? 1. / ntotal : 1); |
1185 | ||
4fa8d795 | 1186 | out->Add(sumCopy); |
1187 | out->Add(sum0Copy); | |
1188 | out->Add(retCopy); | |
1189 | out->Add(sumCopyPx); | |
1190 | out->Add(sum0CopyPx); | |
1191 | out->Add(retCopyPx); | |
1192 | out->Add(norm); | |
1193 | out->Add(norm0); | |
1194 | out->Add(normAll); | |
5ca83fee | 1195 | out->Add(phi); |
1196 | out->Add(phi0); | |
1197 | out->Add(phiAll); | |
4fa8d795 | 1198 | |
c8b1a7db | 1199 | if (fDebug >= 1) { |
1200 | if (n0 > 0) | |
1201 | DMSG(fDebug,1,"Returning (1/%f * %s + 1/%f * %s), " | |
f67d699c | 1202 | "1/%f * %d + 1/%f * %d = %d", |
1203 | epsilon0, fSum0->GetName(), epsilon, fSum->GetName(), | |
1204 | epsilon0, n0, epsilon, n, int(ntotal)); | |
c8b1a7db | 1205 | else |
1206 | DMSG(fDebug,1, "Returning %s, %d", fSum->GetName(), int(ntotal)); | |
1207 | } | |
1208 | ||
4fa8d795 | 1209 | #if 0 |
1210 | for (Int_t i = 1; i <= ret->GetNbinsX(); i++) { | |
1211 | Double_t nc = sum->GetBinContent(i, 0); | |
1212 | Double_t nc0 = sum0->GetBinContent(i, 0); | |
1213 | ret->SetBinContent(i, 0, nc + nc0); // Just count events | |
1214 | } | |
1215 | #endif | |
1216 | ||
1217 | return ret; | |
1218 | } | |
1219 | ||
e1f47419 | 1220 | //==================================================================== |
1221 | AliBasedNdetaTask::CentralityBin::CentralityBin() | |
1222 | : TNamed("", ""), | |
1223 | fSums(0), | |
1224 | fOutput(0), | |
1225 | fSum(0), | |
1226 | fSumMC(0), | |
1227 | fTriggers(0), | |
bfab35d9 | 1228 | fStatus(0), |
e1f47419 | 1229 | fLow(0), |
c89b9ac1 | 1230 | fHigh(0), |
bfab35d9 | 1231 | fDoFinalMCCorrection(false), |
1232 | fSatelliteVertices(false), | |
f53fb4f6 | 1233 | fDebug(0) |
e1f47419 | 1234 | { |
1235 | // | |
1236 | // Constructor | |
1237 | // | |
5ca83fee | 1238 | DGUARD(fDebug,3,"Default CTOR of AliBasedNdeta::CentralityBin"); |
e1f47419 | 1239 | } |
1240 | //____________________________________________________________________ | |
1241 | AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name, | |
1242 | Short_t low, Short_t high) | |
1243 | : TNamed(name, ""), | |
1244 | fSums(0), | |
1245 | fOutput(0), | |
1246 | fSum(0), | |
1247 | fSumMC(0), | |
1248 | fTriggers(0), | |
bfab35d9 | 1249 | fStatus(0), |
e1f47419 | 1250 | fLow(low), |
c89b9ac1 | 1251 | fHigh(high), |
f53fb4f6 | 1252 | fDoFinalMCCorrection(false), |
bfab35d9 | 1253 | fSatelliteVertices(false), |
f53fb4f6 | 1254 | fDebug(0) |
e1f47419 | 1255 | { |
1256 | // | |
1257 | // Constructor | |
1258 | // | |
1259 | // Parameters: | |
1260 | // name Name used for histograms (e.g., Forward) | |
1261 | // low Lower centrality cut in percent | |
1262 | // high Upper centrality cut in percent | |
1263 | // | |
5ca83fee | 1264 | DGUARD(fDebug,3,"Named CTOR of AliBasedNdeta::CentralityBin: %s [%3d,%3d]", |
1265 | name,low,high); | |
e1f47419 | 1266 | if (low <= 0 && high <= 0) { |
1267 | fLow = 0; | |
1268 | fHigh = 0; | |
1269 | SetTitle("All centralities"); | |
1270 | } | |
1271 | else { | |
1272 | fLow = low; | |
1273 | fHigh = high; | |
1274 | SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high)); | |
1275 | } | |
1276 | } | |
1277 | //____________________________________________________________________ | |
1278 | AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o) | |
1279 | : TNamed(o), | |
1280 | fSums(o.fSums), | |
1281 | fOutput(o.fOutput), | |
1282 | fSum(o.fSum), | |
1283 | fSumMC(o.fSumMC), | |
1284 | fTriggers(o.fTriggers), | |
bfab35d9 | 1285 | fStatus(o.fStatus), |
e1f47419 | 1286 | fLow(o.fLow), |
c89b9ac1 | 1287 | fHigh(o.fHigh), |
bfab35d9 | 1288 | fDoFinalMCCorrection(o.fDoFinalMCCorrection), |
1289 | fSatelliteVertices(o.fSatelliteVertices), | |
f53fb4f6 | 1290 | fDebug(o.fDebug) |
e1f47419 | 1291 | { |
1292 | // | |
1293 | // Copy constructor | |
1294 | // | |
1295 | // Parameters: | |
1296 | // other Object to copy from | |
1297 | // | |
5ca83fee | 1298 | DGUARD(fDebug,3,"Copy CTOR of AliBasedNdeta::CentralityBin"); |
e1f47419 | 1299 | } |
1300 | //____________________________________________________________________ | |
1301 | AliBasedNdetaTask::CentralityBin::~CentralityBin() | |
1302 | { | |
1303 | // | |
1304 | // Destructor | |
1305 | // | |
5ca83fee | 1306 | DGUARD(fDebug,3,"DTOR of AliBasedNdeta::CentralityBin"); |
1307 | ||
1308 | // if (fSums) fSums->Delete(); | |
1309 | // if (fOutput) fOutput->Delete(); | |
e1f47419 | 1310 | } |
1311 | ||
1312 | //____________________________________________________________________ | |
1313 | AliBasedNdetaTask::CentralityBin& | |
1314 | AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o) | |
1315 | { | |
1316 | // | |
1317 | // Assignment operator | |
1318 | // | |
1319 | // Parameters: | |
1320 | // other Object to assign from | |
1321 | // | |
1322 | // Return: | |
1323 | // Reference to this | |
1324 | // | |
f53fb4f6 | 1325 | DGUARD(fDebug,3,"Centrality bin assignment"); |
d015ecfe | 1326 | if (&o == this) return *this; |
e1f47419 | 1327 | SetName(o.GetName()); |
1328 | SetTitle(o.GetTitle()); | |
1329 | fSums = o.fSums; | |
1330 | fOutput = o.fOutput; | |
1331 | fSum = o.fSum; | |
1332 | fSumMC = o.fSumMC; | |
1333 | fTriggers = o.fTriggers; | |
bfab35d9 | 1334 | fStatus = o.fStatus; |
e1f47419 | 1335 | fLow = o.fLow; |
1336 | fHigh = o.fHigh; | |
c89b9ac1 | 1337 | fDoFinalMCCorrection = o.fDoFinalMCCorrection; |
bfab35d9 | 1338 | fSatelliteVertices = o.fSatelliteVertices; |
e1f47419 | 1339 | |
1340 | return *this; | |
1341 | } | |
1342 | //____________________________________________________________________ | |
5bb5d1f6 | 1343 | Int_t |
9ecab72f | 1344 | AliBasedNdetaTask::CentralityBin::GetColor(Int_t fallback) const |
5bb5d1f6 | 1345 | { |
9ecab72f | 1346 | if (IsAllBin()) return fallback; |
5bb5d1f6 | 1347 | Float_t fc = (fLow+double(fHigh-fLow)/2) / 100; |
1348 | Int_t nCol = gStyle->GetNumberOfColors(); | |
1349 | Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5)); | |
1350 | Int_t col = gStyle->GetColorPalette(icol); | |
1351 | return col; | |
1352 | } | |
1353 | //____________________________________________________________________ | |
e1f47419 | 1354 | const char* |
1355 | AliBasedNdetaTask::CentralityBin::GetListName() const | |
1356 | { | |
1357 | // | |
1358 | // Get the list name | |
1359 | // | |
1360 | // Return: | |
1361 | // List Name | |
1362 | // | |
1363 | if (IsAllBin()) return "all"; | |
1364 | return Form("cent%03d_%03d", fLow, fHigh); | |
1365 | } | |
1366 | //____________________________________________________________________ | |
1367 | void | |
1a7e2e15 | 1368 | AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir, Int_t mask) |
e1f47419 | 1369 | { |
1370 | // | |
1371 | // Create output objects | |
1372 | // | |
1373 | // Parameters: | |
1374 | // dir Parent list | |
1375 | // | |
f53fb4f6 | 1376 | DGUARD(fDebug,1,"Create centrality bin output objects"); |
e1f47419 | 1377 | fSums = new TList; |
1378 | fSums->SetName(GetListName()); | |
1379 | fSums->SetOwner(); | |
1380 | dir->Add(fSums); | |
1381 | ||
1a7e2e15 | 1382 | fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers", mask); |
ffca499d | 1383 | fTriggers->SetDirectory(0); |
bfab35d9 | 1384 | |
1385 | fStatus = AliAODForwardMult::MakeStatusHistogram("status"); | |
1386 | fStatus->SetDirectory(0); | |
1387 | ||
e1f47419 | 1388 | fSums->Add(fTriggers); |
bfab35d9 | 1389 | fSums->Add(fStatus); |
e1f47419 | 1390 | } |
f53fb4f6 | 1391 | //____________________________________________________________________ |
1392 | void | |
1393 | AliBasedNdetaTask::CentralityBin::SetDebugLevel(Int_t lvl) | |
1394 | { | |
1395 | fDebug = lvl; | |
1396 | if (fSum) fSum->fDebug = lvl; | |
1397 | if (fSumMC) fSumMC->fDebug = lvl; | |
1398 | } | |
1399 | ||
f67d699c | 1400 | //____________________________________________________________________ |
1401 | Bool_t | |
1402 | AliBasedNdetaTask::CentralityBin::ReadSum(TList* list, bool mc) | |
1403 | { | |
1404 | const char* post = (mc ? "MC" : ""); | |
1405 | TString sn = Sum::GetHistName(GetName(),0,post); | |
1406 | TString sn0 = Sum::GetHistName(GetName(),1,post); | |
1407 | TString ev = Sum::GetHistName(GetName(),2,post); | |
c8b1a7db | 1408 | TH2D* sum = static_cast<TH2D*>(list->FindObject(sn)); |
1409 | TH2D* sum0 = static_cast<TH2D*>(list->FindObject(sn0)); | |
1410 | TH1I* events = static_cast<TH1I*>(list->FindObject(ev)); | |
f67d699c | 1411 | if (!sum || !sum0 || !events) { |
c8b1a7db | 1412 | if (!mc) |
1413 | AliWarningF("Failed to find one or more histograms: " | |
1414 | "%s (%p) %s (%p) %s (%p)", | |
1415 | sn.Data(), sum, | |
1416 | sn0.Data(), sum0, | |
1417 | ev.Data(), events); | |
f67d699c | 1418 | return false; |
1419 | } | |
1420 | Sum* ret = new Sum(GetName(), post); | |
1421 | ret->fSum = sum; | |
1422 | ret->fSum0 = sum0; | |
1423 | ret->fEvents = events; | |
1424 | ret->fDebug = fDebug; | |
1425 | if (mc) fSumMC = ret; | |
1426 | else fSum = ret; | |
1427 | ||
1428 | return true; | |
1429 | } | |
1430 | ||
e1f47419 | 1431 | //____________________________________________________________________ |
1432 | void | |
1433 | AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc) | |
1434 | { | |
1435 | // | |
1436 | // Create sum histogram | |
1437 | // | |
1438 | // Parameters: | |
1439 | // data Data histogram to clone | |
1440 | // mc (optional) MC histogram to clone | |
1441 | // | |
f67d699c | 1442 | DGUARD(fDebug,1,"Create centrality bin sums from %s", |
1443 | data ? data->GetName() : "(null)"); | |
9ecab72f | 1444 | if (data) { |
1445 | fSum = new Sum(GetName(),""); | |
1446 | fSum->Init(fSums, data, GetColor()); | |
f53fb4f6 | 1447 | fSum->fDebug = fDebug; |
9ecab72f | 1448 | } |
1449 | ||
e1f47419 | 1450 | // If no MC data is given, then do not create MC sum histogram |
1451 | if (!mc) return; | |
1452 | ||
4fa8d795 | 1453 | fSumMC = new Sum(GetName(), "MC"); |
1454 | fSumMC->Init(fSums, mc, GetColor()); | |
f53fb4f6 | 1455 | fSumMC->fDebug = fDebug; |
e1f47419 | 1456 | } |
1457 | ||
1458 | //____________________________________________________________________ | |
1459 | Bool_t | |
1460 | AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward, | |
1461 | Int_t triggerMask, | |
1462 | Double_t vzMin, Double_t vzMax) | |
1463 | { | |
1464 | // | |
1465 | // Check the trigger, vertex, and centrality | |
1466 | // | |
1467 | // Parameters: | |
1468 | // aod Event input | |
1469 | // | |
1470 | // Return: | |
1471 | // true if the event is to be used | |
1472 | // | |
1473 | if (!forward) return false; | |
e1f47419 | 1474 | |
f53fb4f6 | 1475 | DGUARD(fDebug,2,"Check the event"); |
ffca499d | 1476 | // We do not check for centrality here - it's already done |
bfab35d9 | 1477 | return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, |
1478 | fTriggers, fStatus); | |
e1f47419 | 1479 | } |
1480 | ||
1481 | ||
1482 | //____________________________________________________________________ | |
bfab35d9 | 1483 | Bool_t |
e1f47419 | 1484 | AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward, |
4fa8d795 | 1485 | Int_t triggerMask, Bool_t isZero, |
e1f47419 | 1486 | Double_t vzMin, Double_t vzMax, |
1487 | const TH2D* data, const TH2D* mc) | |
1488 | { | |
1489 | // | |
1490 | // Process an event | |
1491 | // | |
1492 | // Parameters: | |
1493 | // forward Forward data (for trigger, vertex, & centrality) | |
1494 | // triggerMask Trigger mask | |
1495 | // vzMin Minimum IP z coordinate | |
1496 | // vzMax Maximum IP z coordinate | |
1497 | // data Data histogram | |
1498 | // mc MC histogram | |
1499 | // | |
f53fb4f6 | 1500 | DGUARD(fDebug,1,"Process one event for %s a given centrality bin", |
f67d699c | 1501 | data ? data->GetName() : "(null)"); |
bfab35d9 | 1502 | if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return false; |
1503 | if (!data) return false; | |
e1f47419 | 1504 | if (!fSum) CreateSums(data, mc); |
4fa8d795 | 1505 | |
1506 | fSum->Add(data, isZero); | |
1507 | if (mc) fSumMC->Add(mc, isZero); | |
bfab35d9 | 1508 | |
1509 | return true; | |
e1f47419 | 1510 | } |
1511 | ||
1512 | //________________________________________________________________________ | |
b30dee70 | 1513 | Double_t |
1514 | AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t, | |
1515 | UShort_t scheme, | |
1516 | Double_t trigEff, | |
1a7e2e15 | 1517 | Double_t& ntotal, |
1518 | TString* text) const | |
e1f47419 | 1519 | { |
1520 | // | |
b30dee70 | 1521 | // Calculate normalization |
e1f47419 | 1522 | // |
b30dee70 | 1523 | // Parameters: |
1524 | // t Trigger histogram | |
1525 | // scheme Normaliztion scheme | |
1526 | // trigEff From MC | |
1527 | // ntotal On return, contains the number of events. | |
e1f47419 | 1528 | // |
e12f8c42 | 1529 | DGUARD(fDebug,1,"Normalize centrality bin %s [%3d-%3d%%] with %s", |
1530 | GetName(), fLow, fHigh, t.GetName()); | |
0be6c8cd | 1531 | Double_t nAll = t.GetBinContent(AliAODForwardMult::kBinAll); |
1532 | Double_t nB = t.GetBinContent(AliAODForwardMult::kBinB); | |
1533 | Double_t nA = t.GetBinContent(AliAODForwardMult::kBinA); | |
1534 | Double_t nC = t.GetBinContent(AliAODForwardMult::kBinC); | |
1535 | Double_t nE = t.GetBinContent(AliAODForwardMult::kBinE); | |
1536 | Double_t nOffline = t.GetBinContent(AliAODForwardMult::kBinOffline); | |
1537 | Double_t nTriggered = t.GetBinContent(AliAODForwardMult::kWithTrigger); | |
1538 | Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex); | |
66cf95f2 | 1539 | Double_t nAccepted = ntotal; |
4fa8d795 | 1540 | ntotal = 0; |
e1f47419 | 1541 | |
e308a636 | 1542 | if (nTriggered <= 0.1) { |
e1f47419 | 1543 | AliError("Number of triggered events <= 0"); |
b30dee70 | 1544 | return -1; |
e1f47419 | 1545 | } |
e308a636 | 1546 | if (nWithVertex <= 0.1) { |
1547 | AliError("Number of events with vertex <= 0"); | |
b30dee70 | 1548 | return -1; |
e1f47419 | 1549 | } |
b30dee70 | 1550 | ntotal = nAccepted; |
1551 | Double_t vtxEff = nWithVertex / nTriggered; | |
0be6c8cd | 1552 | Double_t scaler = 1; |
0be6c8cd | 1553 | Double_t beta = nA + nC - 2*nE; |
b30dee70 | 1554 | |
1a7e2e15 | 1555 | |
1556 | TString rhs("N = N_acc"); | |
1557 | if (!(scheme & kZeroBin)) { | |
1558 | if (scheme & kEventLevel) { | |
1559 | ntotal = nAccepted / vtxEff; | |
1560 | scaler = vtxEff; | |
c8b1a7db | 1561 | DMSG(fDebug,0,"Calculating event normalisation as\n" |
1562 | " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)", | |
1563 | Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex), | |
1564 | ntotal, scaler); | |
1a7e2e15 | 1565 | if (scheme & kBackground) { |
1566 | // 1 E_V E_V | |
1567 | // s = --------- = ------------- = ------------ | |
1568 | // 1 - beta 1 - beta E_V 1 - beta N_V | |
1569 | // --- ---- -------- ---- --- | |
1570 | // E_V N_V N_V N_V N_T | |
1571 | // | |
1572 | // E_V | |
1573 | // = ------------ | |
1574 | // 1 - beta | |
1575 | // ---- | |
1576 | // N_T | |
1577 | // | |
1578 | ntotal -= nAccepted * beta / nWithVertex; | |
1579 | // This one is direct and correct. | |
1580 | // scaler = 1. / (1. / vtxEff - beta / nWithVertex); | |
1581 | // A simpler expresion | |
1582 | scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689 | |
c8b1a7db | 1583 | DMSG(fDebug,0,"Calculating event normalisation as\n" |
1584 | " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n" | |
1585 | " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)", | |
1586 | Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta), | |
1587 | nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta), | |
1588 | Int_t(nWithVertex), ntotal, scaler); | |
1a7e2e15 | 1589 | rhs.Append("(1/eps_V - beta/N_vtx)"); |
1590 | } // Background | |
1591 | else | |
1592 | rhs.Append("/eps_V"); | |
1593 | } // Event-level | |
1594 | if (scheme & kTriggerEfficiency) { | |
77f97e3f | 1595 | Int_t old = Int_t(ntotal); |
c8b1a7db | 1596 | ntotal /= trigEff; |
1597 | scaler *= trigEff; | |
1598 | DMSG(fDebug,0,"Correcting for trigger efficiency:\n" | |
1599 | " N = 1 / E_X * N = 1 / %f * %d = %f (%f)", | |
1600 | trigEff, old, ntotal, scaler); | |
1a7e2e15 | 1601 | rhs.Append("/eps_T"); |
1602 | } // Trigger efficiency | |
1603 | } | |
1604 | else { | |
4fa8d795 | 1605 | // Calculate as |
1606 | // | |
1607 | // N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta) | |
1608 | // = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) | |
1a7e2e15 | 1609 | // = N_A (1 + 1/E_X (1/E_V - 1 - beta / N_V)) |
4fa8d795 | 1610 | // |
1611 | // s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) | |
1612 | // = N_V / (N_V + 1/E_X (N_T - N_V - beta)) | |
1613 | // | |
1614 | if (!(scheme & kBackground)) beta = 0; | |
1615 | ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1 | |
1616 | - beta / nWithVertex)); | |
1617 | scaler = nWithVertex / (nWithVertex + | |
1618 | 1/trigEff * (nTriggered-nWithVertex-beta)); | |
c8b1a7db | 1619 | DMSG(fDebug,0,"Calculating event normalisation as\n" |
1620 | " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n" | |
1621 | " N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = " | |
1622 | "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)", | |
1623 | Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta), | |
1624 | Int_t(nAccepted), trigEff, Int_t(nTriggered), | |
1625 | Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex), | |
1626 | ntotal, scaler); | |
1a7e2e15 | 1627 | rhs.Append("(1+1/eps_T(1/eps_V-1-beta/N_vtx))"); |
4fa8d795 | 1628 | } |
1a7e2e15 | 1629 | |
1630 | if (text) { | |
c8b1a7db | 1631 | text->Append(Form("%-40s = %d\n", "N_all", UInt_t(nAll))); |
1632 | text->Append(Form("%-40s = %d\n", "N_acc", UInt_t(nAccepted))); | |
1633 | text->Append(Form("%-40s = %d\n", "N_trg", UInt_t(nTriggered))); | |
1634 | text->Append(Form("%-40s = %d\n", "N_vtx", UInt_t(nWithVertex))); | |
1635 | text->Append(Form("%-40s = %d\n", "N_B", UInt_t(nB))); | |
1636 | text->Append(Form("%-40s = %d\n", "N_A", UInt_t(nA))); | |
1637 | text->Append(Form("%-40s = %d\n", "N_C", UInt_t(nC))); | |
1638 | text->Append(Form("%-40s = %d\n", "N_E", UInt_t(nE))); | |
1a7e2e15 | 1639 | text->Append(Form("%-40s = %d\n", "beta = N_A + N_C - 2N_E",UInt_t(beta))); |
1640 | text->Append(Form("%-40s = %f\n", "eps_V = N_vtx/N_trg",vtxEff)); | |
c8b1a7db | 1641 | text->Append(Form("%-40s = %f\n", "eps_T", trigEff)); |
1642 | text->Append(Form("%-40s = %f\n", rhs.Data(), ntotal)); | |
ffca499d | 1643 | } |
c8b1a7db | 1644 | TString tN = t.GetXaxis()->GetBinLabel(AliAODForwardMult::kWithTrigger); |
1645 | tN.ReplaceAll("w/Selected trigger (",""); | |
1646 | tN.ReplaceAll(")", ""); | |
1647 | DMSG(fDebug,0,"\n" | |
1648 | " Total of %9d events for %s\n" | |
1649 | " of these %9d have an offline trigger\n" | |
1650 | " of these N_T = %9d has the selected trigger (%s)\n" | |
1651 | " of these N_V = %9d has a vertex\n" | |
1652 | " of these N_A = %9d were in the selected range\n" | |
1653 | " Triggers by hardware type:\n" | |
1654 | " N_b = %9d\n" | |
1655 | " N_ac = %9d (%d+%d)\n" | |
0be6c8cd | 1656 | " N_e = %9d\n" |
c8b1a7db | 1657 | " Vertex efficiency: %f\n" |
1658 | " Trigger efficiency: %f\n" | |
1659 | " Total number of events: N = %f\n" | |
1660 | " Scaler (N_A/N): %f\n" | |
1661 | " %25s = %f", | |
1662 | Int_t(nAll), GetTitle(), Int_t(nOffline), | |
1663 | Int_t(nTriggered), tN.Data(), | |
1664 | Int_t(nWithVertex), Int_t(nAccepted), | |
1665 | Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE), | |
1666 | vtxEff, trigEff, ntotal, scaler, rhs.Data(), ntotal); | |
b30dee70 | 1667 | return scaler; |
1668 | } | |
1669 | ||
5bb5d1f6 | 1670 | //________________________________________________________________________ |
1671 | const char* | |
1672 | AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin, | |
1673 | Bool_t sym, | |
1674 | const char* postfix) const | |
1675 | { | |
1676 | static TString n; | |
c8b1a7db | 1677 | n = GetName(); |
1678 | n.ReplaceAll("dNdeta", ""); | |
1679 | n.Prepend("dndeta"); | |
1680 | n.Append(postfix); | |
1681 | // n = Form("dndeta%s%s",GetName(), postfix); | |
5bb5d1f6 | 1682 | if (rebin > 1) n.Append(Form("_rebin%02d", rebin)); |
1683 | if (sym) n.Append("_mirror"); | |
1684 | return n.Data(); | |
1685 | } | |
1686 | //________________________________________________________________________ | |
1687 | TH1* | |
c8b1a7db | 1688 | AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin, |
1689 | Bool_t sym, | |
1690 | const char* postfix, | |
1691 | Bool_t verbose) const | |
5bb5d1f6 | 1692 | { |
1693 | if (!fOutput) { | |
c8b1a7db | 1694 | AliWarningF("No output list defined in %s [%3d,%3d]", GetName(), |
1695 | fLow, fHigh); | |
5bb5d1f6 | 1696 | return 0; |
1697 | } | |
1698 | TString n = GetResultName(rebin, sym, postfix); | |
1699 | TObject* o = fOutput->FindObject(n.Data()); | |
1700 | if (!o) { | |
c8b1a7db | 1701 | if (verbose) |
1702 | AliWarningF("Object %s not found in output list of %s", | |
1703 | n.Data(), GetName()); | |
5bb5d1f6 | 1704 | return 0; |
1705 | } | |
1706 | return static_cast<TH1*>(o); | |
1707 | } | |
1708 | ||
c25b5e1b | 1709 | //________________________________________________________________________ |
1710 | void | |
1711 | AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, | |
1712 | const char* postfix, | |
1713 | bool rootProj, | |
1714 | bool corrEmpty, | |
f67d699c | 1715 | const TH2F* shapeCorr, |
c25b5e1b | 1716 | Double_t scaler, |
1717 | bool symmetrice, | |
1718 | Int_t rebin, | |
1719 | bool cutEdges, | |
9ecab72f | 1720 | Int_t marker, |
c89b9ac1 | 1721 | Int_t color, |
1722 | TList* mclist, | |
1723 | TList* truthlist) | |
c25b5e1b | 1724 | { |
1725 | // | |
1726 | // Generate the dN/deta result from input | |
1727 | // | |
1728 | // Parameters: | |
1729 | // sum Sum of 2D hists | |
1730 | // postfix Post fix on names | |
1731 | // rootProj Whether to use ROOT TH2::ProjectionX | |
1732 | // corrEmpty Correct for empty bins | |
1733 | // shapeCorr Shape correction to use | |
1734 | // scaler Event-level normalization scaler | |
1735 | // symmetrice Whether to make symmetric extensions | |
1736 | // rebin Whether to rebin | |
1737 | // cutEdges Whether to cut edges when rebinning | |
1738 | // | |
f53fb4f6 | 1739 | DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName()); |
c8b1a7db | 1740 | TString base(GetName()); |
1741 | base.ReplaceAll("dNdeta", ""); | |
1742 | base.Append(postfix); | |
1743 | TH2D* copy = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s", | |
1744 | base.Data()))); | |
bfab35d9 | 1745 | |
1746 | TH1D* accNorm = 0; | |
5ca83fee | 1747 | Int_t nY = sum->GetNbinsY(); |
bfab35d9 | 1748 | // Hack HHD Hans test code to manually remove FMD2 dead channel (but |
1749 | // it is on outer?) | |
1750 | // | |
1751 | // cholm comment: The original hack has been moved to | |
1752 | // AliForwarddNdetaTask::CheckEventData - this simplifies things a | |
1753 | // great deal, and we could in principle use the new phi acceptance. | |
1754 | // | |
1755 | // However, since we may have filtered out the dead sectors in the | |
1756 | // AOD production already, we can't be sure we can recalculate the | |
1757 | // phi acceptance correctly, so for now, we rely on fCorrEmpty being set. | |
5ca83fee | 1758 | Int_t o = (corrEmpty ? 0 : nY+1); |
c8b1a7db | 1759 | accNorm = ProjectX(sum, Form("norm%s",base.Data()), o, o, |
bfab35d9 | 1760 | rootProj, corrEmpty, false); |
c25b5e1b | 1761 | accNorm->SetDirectory(0); |
1762 | ||
1763 | // ---- Scale by shape correction ---------------------------------- | |
1764 | if (shapeCorr) copy->Divide(shapeCorr); | |
c8b1a7db | 1765 | else DMSG(fDebug,1,"No shape correction specified, or disabled"); |
c25b5e1b | 1766 | |
4fa8d795 | 1767 | // --- Normalize to the coverage ----------------------------------- |
5ca83fee | 1768 | if (corrEmpty) { |
1769 | ScaleToCoverage(copy, accNorm); | |
1770 | // --- Event-level normalization --------------------------------- | |
1771 | copy->Scale(scaler); | |
1772 | } | |
c25b5e1b | 1773 | |
1774 | // --- Project on X axis ------------------------------------------- | |
c8b1a7db | 1775 | TH1D* dndeta = ProjectX(copy, Form("dndeta%s",base.Data()), |
5ca83fee | 1776 | 1, nY, rootProj, corrEmpty); |
c25b5e1b | 1777 | dndeta->SetDirectory(0); |
1778 | // Event-level normalization | |
5ca83fee | 1779 | if (!corrEmpty) { |
1780 | ScaleToCoverage(dndeta, accNorm); | |
1781 | dndeta->Scale(scaler); | |
1782 | } | |
c25b5e1b | 1783 | dndeta->Scale(1., "width"); |
1784 | copy->Scale(1., "width"); | |
c89b9ac1 | 1785 | |
c8b1a7db | 1786 | TH1D* dndetaMCCorrection = 0; |
1787 | TH1D* dndetaMCtruth = 0; | |
1788 | TList* centlist = 0; | |
1789 | TList* truthcentlist = 0; | |
c89b9ac1 | 1790 | |
bfab35d9 | 1791 | // --- Possible final correction to <MC analysis> / <MC truth> ----- |
1792 | // we get the rebinned distribution for satellite to make the correction | |
1793 | TString rebinSuf(fSatelliteVertices ? "_rebin05" : ""); | |
c8b1a7db | 1794 | if(mclist) { |
c89b9ac1 | 1795 | centlist = static_cast<TList*> (mclist->FindObject(GetListName())); |
c8b1a7db | 1796 | if(centlist) |
1797 | dndetaMCCorrection = | |
1798 | static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s", | |
1799 | base.Data(), | |
1800 | rebinSuf.Data()))); | |
1801 | } | |
1802 | if (truthlist) { | |
1803 | truthcentlist = static_cast<TList*>(truthlist->FindObject(GetListName())); | |
1804 | if (truthcentlist) | |
1805 | // TODO here new is "dndetaTruth" | |
1806 | dndetaMCtruth = | |
1807 | static_cast<TH1D*>(truthcentlist->FindObject(Form("dndetaMCTruth%s", | |
1808 | rebinSuf.Data()))); | |
1809 | } | |
bfab35d9 | 1810 | |
1811 | if (dndetaMCCorrection && dndetaMCtruth) { | |
c89b9ac1 | 1812 | AliInfo("Correcting with final MC correction"); |
bfab35d9 | 1813 | TString testString(dndetaMCCorrection->GetName()); |
1814 | ||
1815 | // We take the measured MC dN/deta and divide with truth | |
c89b9ac1 | 1816 | dndetaMCCorrection->Divide(dndetaMCtruth); |
dd3fe5b6 | 1817 | dndetaMCCorrection->SetTitle("Final MC correction"); |
1818 | dndetaMCCorrection->SetName("finalMCCorr"); | |
bfab35d9 | 1819 | for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) { |
1820 | if(dndetaMCCorrection->GetBinContent(m) < 0.5 || | |
1821 | dndetaMCCorrection->GetBinContent(m) > 1.75) { | |
1822 | dndetaMCCorrection->SetBinContent(m,1.); | |
1823 | dndetaMCCorrection->SetBinError(m,0.1); | |
1824 | } | |
1825 | } | |
1826 | // Applying the correction | |
1827 | if (!fSatelliteVertices) | |
1828 | // For non-satellites we took the same binning, so we do a straight | |
1829 | // division | |
1830 | dndeta->Divide(dndetaMCCorrection); | |
1831 | else { | |
1832 | // For satelitte events, we took coarser binned histograms, so | |
1833 | // we need to do a bit more | |
1834 | for(Int_t m = 1; m <= dndeta->GetNbinsX(); m++) { | |
1835 | if(dndeta->GetBinContent(m) <= 0.01 ) continue; | |
1836 | ||
1837 | Double_t eta = dndeta->GetXaxis()->GetBinCenter(m); | |
1838 | Int_t bin = dndetaMCCorrection->GetXaxis()->FindBin(eta); | |
1839 | Double_t mccorr = dndetaMCCorrection->GetBinContent(bin); | |
1840 | Double_t mcerror = dndetaMCCorrection->GetBinError(bin); | |
1841 | if (mccorr < 1e-6) { | |
1842 | dndeta->SetBinContent(m, 0); | |
1843 | dndeta->SetBinError(m, 0); | |
1844 | } | |
1845 | Double_t value = dndeta->GetBinContent(m); | |
1846 | Double_t error = dndeta->GetBinError(m); | |
1847 | Double_t sumw2 = (error * error * mccorr * mccorr + | |
1848 | mcerror * mcerror * value * value); | |
1849 | dndeta->SetBinContent(m,value/mccorr) ; | |
1850 | dndeta->SetBinError(m,TMath::Sqrt(sumw2)/mccorr/mccorr); | |
1851 | } | |
1852 | } | |
c89b9ac1 | 1853 | } |
66cf95f2 | 1854 | else |
c8b1a7db | 1855 | DMSG(fDebug,1,"No final MC correction applied"); |
c89b9ac1 | 1856 | |
c25b5e1b | 1857 | // --- Set some histogram attributes ------------------------------- |
5bb5d1f6 | 1858 | TString post; |
9ecab72f | 1859 | Int_t rColor = GetColor(color); |
5bb5d1f6 | 1860 | if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix); |
9ecab72f | 1861 | SetHistogramAttributes(dndeta, rColor, marker, |
5bb5d1f6 | 1862 | Form("ALICE %s%s", GetName(), post.Data())); |
9ecab72f | 1863 | SetHistogramAttributes(accNorm, rColor, marker, |
5bb5d1f6 | 1864 | Form("ALICE %s normalisation%s", |
1865 | GetName(), post.Data())); | |
c25b5e1b | 1866 | |
1867 | // --- Make symmetric extensions and rebinnings -------------------- | |
bfab35d9 | 1868 | if (symmetrice) fOutput->Add(Symmetrice(dndeta)); |
c25b5e1b | 1869 | fOutput->Add(dndeta); |
1870 | fOutput->Add(accNorm); | |
1871 | fOutput->Add(copy); | |
1872 | fOutput->Add(Rebin(dndeta, rebin, cutEdges)); | |
bfab35d9 | 1873 | if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges))); |
dd3fe5b6 | 1874 | if (dndetaMCCorrection) fOutput->Add(dndetaMCCorrection); |
bfab35d9 | 1875 | |
1876 | // HHD Test of dN/deta in phi bins add flag later? | |
1877 | // | |
1878 | // cholm comment: We disable this for now | |
1879 | #if 0 | |
1880 | for (Int_t nn=1; nn <= sum->GetNbinsY(); nn++) { | |
c8b1a7db | 1881 | TH1D* dndeta_phi = ProjectX(copy, Form("dndeta%s_phibin%d", |
1882 | base.Data(), nn), | |
bfab35d9 | 1883 | nn, nn, rootProj, corrEmpty); |
1884 | dndeta_phi->SetDirectory(0); | |
1885 | // Event-level normalization | |
1886 | dndeta_phi->Scale(TMath::Pi()/10., "width"); | |
1887 | ||
1888 | if(centlist) | |
1889 | dndetaMCCorrection = | |
c8b1a7db | 1890 | static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s_phibin%d", |
1891 | base.Data(),nn))); | |
bfab35d9 | 1892 | if(truthcentlist) |
1893 | dndetaMCtruth | |
1894 | = static_cast<TH1D*>(truthcentlist->FindObject("dndetaMCTruth")); | |
1895 | ||
1896 | if (dndetaMCCorrection && dndetaMCtruth) { | |
1897 | AliInfo("Correcting with final MC correction"); | |
1898 | TString testString(dndetaMCCorrection->GetName()); | |
1899 | dndetaMCCorrection->Divide(dndetaMCtruth); | |
1900 | dndetaMCCorrection->SetTitle(Form("Final_MC_correction_phibin%d",nn)); | |
1901 | dndetaMCCorrection->SetName(Form("Final_MC_correction_phibin%d",nn)); | |
1902 | for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) { | |
1903 | if(dndetaMCCorrection->GetBinContent(m) < 0.25 || | |
1904 | dndetaMCCorrection->GetBinContent(m) > 1.75) { | |
1905 | dndetaMCCorrection->SetBinContent(m,1.); | |
1906 | dndetaMCCorrection->SetBinError(m,0.1); | |
1907 | } | |
1908 | } | |
1909 | //Applying the correction | |
1910 | dndeta_phi->Divide(dndetaMCCorrection); | |
1911 | } | |
1912 | fOutput->Add(dndeta_phi); | |
1913 | fOutput->Add(Rebin(dndeta_phi, rebin, cutEdges)); | |
1914 | if(dndetaMCCorrection) fOutput->Add(dndetaMCCorrection); | |
1915 | } // End of phi | |
1916 | #endif | |
c25b5e1b | 1917 | } |
1918 | ||
b30dee70 | 1919 | //________________________________________________________________________ |
1920 | void | |
1921 | AliBasedNdetaTask::CentralityBin::End(TList* sums, | |
1922 | TList* results, | |
1923 | UShort_t scheme, | |
f67d699c | 1924 | const TH2F* shapeCorr, |
b30dee70 | 1925 | Double_t trigEff, |
66cf95f2 | 1926 | Double_t trigEff0, |
b30dee70 | 1927 | Bool_t symmetrice, |
1928 | Int_t rebin, | |
c25b5e1b | 1929 | Bool_t rootProj, |
b30dee70 | 1930 | Bool_t corrEmpty, |
1931 | Bool_t cutEdges, | |
1932 | Int_t triggerMask, | |
9ecab72f | 1933 | Int_t marker, |
c89b9ac1 | 1934 | Int_t color, |
1935 | TList* mclist, | |
1936 | TList* truthlist) | |
b30dee70 | 1937 | { |
1938 | // | |
1939 | // End of processing | |
1940 | // | |
1941 | // Parameters: | |
1942 | // sums List of sums | |
1943 | // results Output list of results | |
1944 | // shapeCorr Shape correction or nil | |
1945 | // trigEff Trigger efficiency | |
1946 | // symmetrice Whether to symmetrice the results | |
1947 | // rebin Whether to rebin the results | |
1948 | // corrEmpty Whether to correct for empty bins | |
1949 | // cutEdges Whether to cut edges when rebinning | |
1950 | // triggerMask Trigger mask | |
1951 | // | |
f53fb4f6 | 1952 | DGUARD(fDebug,1,"End centrality bin procesing"); |
b30dee70 | 1953 | |
1954 | fSums = dynamic_cast<TList*>(sums->FindObject(GetListName())); | |
1955 | if(!fSums) { | |
1956 | AliError("Could not retrieve TList fSums"); | |
1957 | return; | |
1958 | } | |
e1f47419 | 1959 | |
b30dee70 | 1960 | fOutput = new TList; |
1961 | fOutput->SetName(GetListName()); | |
1962 | fOutput->SetOwner(); | |
1963 | results->Add(fOutput); | |
1964 | ||
9ecab72f | 1965 | if (!fSum) { |
f67d699c | 1966 | if (!ReadSum(fSums, false)) { |
1967 | AliInfo("This task did not produce any output"); | |
1968 | return; | |
1969 | } | |
9ecab72f | 1970 | } |
f67d699c | 1971 | if (!fSumMC) ReadSum(fSums, true); |
9ecab72f | 1972 | |
b30dee70 | 1973 | fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers")); |
1974 | ||
1975 | if (!fTriggers) { | |
1976 | AliError("Couldn't find histogram 'triggers' in list"); | |
1977 | return; | |
1978 | } | |
b30dee70 | 1979 | |
b30dee70 | 1980 | // --- Get normalization scaler ------------------------------------ |
4fa8d795 | 1981 | Double_t epsilonT = trigEff; |
66cf95f2 | 1982 | Double_t epsilonT0 = trigEff0; |
c8b1a7db | 1983 | DMSG(fDebug,2,"Using epsilonT=%f, epsilonT0=%f for 0x%x", |
1984 | epsilonT, epsilonT0, triggerMask); | |
4fa8d795 | 1985 | |
1986 | // Get our histograms | |
1987 | Double_t nSum = 0; | |
66cf95f2 | 1988 | TH2D* sum = fSum->CalcSum(fOutput, nSum, epsilonT0, epsilonT, |
c8b1a7db | 1989 | marker, rootProj, corrEmpty); |
4fa8d795 | 1990 | Double_t nSumMC = 0; |
1991 | TH2D* sumMC = 0; | |
f67d699c | 1992 | if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC, |
66cf95f2 | 1993 | epsilonT0, epsilonT, marker, |
f67d699c | 1994 | rootProj, corrEmpty); |
4fa8d795 | 1995 | if (!sum) { |
1996 | AliError("Failed to get sum from summer - bailing out"); | |
1997 | return; | |
1998 | } | |
1999 | ||
1a7e2e15 | 2000 | TString text; |
4fa8d795 | 2001 | Double_t ntotal = nSum; |
1a7e2e15 | 2002 | Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text); |
b30dee70 | 2003 | if (scaler < 0) { |
2004 | AliError("Failed to calculate normalization - bailing out"); | |
2005 | return; | |
2006 | } | |
e1f47419 | 2007 | fOutput->Add(fTriggers->Clone()); |
1a7e2e15 | 2008 | fOutput->Add(new TNamed("normCalc", text.Data())); |
c25b5e1b | 2009 | |
2010 | // --- Make result and store --------------------------------------- | |
4fa8d795 | 2011 | MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0, |
797161e8 | 2012 | scaler, symmetrice, rebin, cutEdges, marker, color, |
2013 | mclist, truthlist); | |
e1f47419 | 2014 | |
b30dee70 | 2015 | // --- Process result from TrackRefs ------------------------------- |
4fa8d795 | 2016 | if (sumMC) |
2017 | MakeResult(sumMC, "MC", rootProj, corrEmpty, | |
c25b5e1b | 2018 | (scheme & kShape) ? shapeCorr : 0, |
4fa8d795 | 2019 | scaler, symmetrice, rebin, cutEdges, |
797161e8 | 2020 | GetMarkerStyle(GetMarkerBits(marker)+4), color, |
2021 | mclist, truthlist); | |
4fa8d795 | 2022 | |
e308a636 | 2023 | // Temporary stuff |
5bb5d1f6 | 2024 | // if (!IsAllBin()) return; |
e308a636 | 2025 | |
e1f47419 | 2026 | } |
8449e3e0 | 2027 | //____________________________________________________________________ |
2028 | Bool_t | |
2029 | AliBasedNdetaTask::ApplyEmpiricalCorrection(const AliAODForwardMult* aod, | |
2030 | TH2D* data) | |
3d9b0442 | 2031 | { |
8449e3e0 | 2032 | if (!fglobalempiricalcorrection || !data) |
2033 | return true; | |
2034 | Float_t zvertex=aod->GetIpZ(); | |
2035 | Int_t binzvertex=fglobalempiricalcorrection->GetXaxis()->FindBin(zvertex); | |
2036 | if(binzvertex<1||binzvertex>fglobalempiricalcorrection->GetNbinsX()) | |
2037 | return false; | |
2038 | for (int i=1;i<=data->GetNbinsX();i++) { | |
2039 | Int_t bincorrection=fglobalempiricalcorrection->GetYaxis() | |
2040 | ->FindBin(data->GetXaxis()->GetBinCenter(i)); | |
2041 | if(bincorrection<1||bincorrection>fglobalempiricalcorrection->GetNbinsY()) | |
2042 | return false; | |
2043 | Float_t correction=fglobalempiricalcorrection | |
2044 | ->GetBinContent(binzvertex,bincorrection); | |
2045 | if(correction<0.001) { | |
2046 | data->SetBinContent(i,0,0); | |
2047 | data->SetBinContent(i,data->GetNbinsY()+1,0); | |
2048 | } | |
2049 | for(int j=1;j<=data->GetNbinsY();j++) { | |
2050 | if (data->GetBinContent(i,j)>0.0) { | |
2051 | data->SetBinContent(i,j,data->GetBinContent(i,j)*correction); | |
2052 | data->SetBinError(i,j,data->GetBinError(i,j)*correction); | |
2053 | } | |
2054 | } | |
2055 | } | |
2056 | return true; | |
3d9b0442 | 2057 | } |
e1f47419 | 2058 | |
fb3430ac | 2059 | // |
2060 | // EOF | |
2061 | // |