]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/AliForwardMultDists.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Attempt to enable monitor objects in Proof(Lite) - doesn't work
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardMultDists.cxx
CommitLineData
bfab35d9 1#include "AliForwardMultDists.h"
2#include <AliForwardUtil.h>
3#include <AliAODForwardMult.h>
4#include <AliAODCentralMult.h>
5#include <AliAODEvent.h>
6#include <AliLog.h>
7#include <TH1.h>
8#include <TH2.h>
9#include <TVector2.h>
10#include <THStack.h>
11#include <TParameter.h>
52b36573 12#include <TError.h>
bfab35d9 13#include <iostream>
14#include <iomanip>
15
c8b1a7db 16//====================================================================
17namespace {
18 /**
19 * Marker styles
20 */
21 enum {
22 kSolid = 0x000,
23 kHollow = 0x001,
24 kCircle = 0x002,
25 kSquare = 0x004,
26 kUpTriangle = 0x006,
27 kDownTriangle = 0x008,
28 kDiamond = 0x00a,
29 kCross = 0x00c,
30 kStar = 0x00e
31 };
32 /**
33 * Get a ROOT marker style from bit options
34 *
35 * @param bits Bit mask of options
36 *
37 * @return ROOT marker style
38 */
39 Int_t GetMarkerStyle(UShort_t bits)
40 {
41 Int_t base = bits & (0xFE);
42 Bool_t hollow = bits & kHollow;
43 switch (base) {
44 case kCircle: return (hollow ? 24 : 20);
45 case kSquare: return (hollow ? 25 : 21);
46 case kUpTriangle: return (hollow ? 26 : 22);
47 case kDownTriangle: return (hollow ? 32 : 23);
48 case kDiamond: return (hollow ? 27 : 33);
49 case kCross: return (hollow ? 28 : 34);
50 case kStar: return (hollow ? 30 : 29);
51 }
52 return 1;
53 }
54 /**
55 * Get the marker option bits from a ROOT style
56 *
57 * @param style ROOT style
58 *
59 * @return Pattern of marker options
60 */
61 UShort_t GetMarkerBits(Int_t style)
62 {
63 UShort_t bits = 0;
64 switch (style) {
65 case 24: case 25: case 26: case 27: case 28: case 30: case 32:
66 bits |= kHollow; break;
67 }
68 switch (style) {
69 case 20: case 24: bits |= kCircle; break;
70 case 21: case 25: bits |= kSquare; break;
71 case 22: case 26: bits |= kUpTriangle; break;
72 case 23: case 32: bits |= kDownTriangle; break;
73 case 27: case 33: bits |= kDiamond; break;
74 case 28: case 34: bits |= kCross; break;
75 case 29: case 30: bits |= kStar; break;
76 }
77 return bits;
78 }
79 static Int_t GetIndexMarker(Int_t idx)
80 {
81 const UShort_t nMarkers = 7;
82 UShort_t markers[] = {
83 kCircle,
84 kSquare,
85 kUpTriangle,
86 kDownTriangle,
87 kDiamond,
88 kCross,
89 kStar
90 };
91
92 return markers[idx % nMarkers];
93 }
94 /**
95 * Flip the 'hollow-ness' of a marker
96 *
97 * @param style ROOT Style
98 *
99 * @return ROOT style
100 */
101 Int_t FlipHollowStyle(Int_t style)
102 {
103 UShort_t bits = GetMarkerBits(style);
104 Int_t ret = GetMarkerStyle(bits ^ kHollow);
105 return ret;
106 }
107}
108
109//====================================================================
52b36573 110AliForwardMultDists::BinSpec::BinSpec(Double_t etaMin,
111 Double_t etaMax,
112 Double_t nchLow)
113 : fEtaMin(etaMin), fEtaMax(etaMax), fLow(nchLow), fN(0), fD(0), fAxis(1,0,1)
114{
115}
116//____________________________________________________________________
117void
118AliForwardMultDists::BinSpec::Push(UShort_t n, Double_t d)
119{
120 Int_t l = fN.GetSize();
121 fN.Set(l+1); // Not terribly efficient, but that's ok because we do
122 fD.Set(l+1); // this infrequently
123 fN[l] = n;
124 fD[l] = d;
125}
126//____________________________________________________________________
127const TAxis&
128AliForwardMultDists::BinSpec::Axis() const
129{
130 if (fAxis.GetNbins() > 1) return fAxis;
131 if (fN.GetSize() <= 0) return fAxis;
132 if (fN.GetSize() == 1) {
133 // Exactly one spec,
134 fAxis.Set(fN[0], fLow, fN[0]*fD[0]+fLow);
135 }
136 else {
137 Int_t n = fN.GetSum();
138 Int_t l = fN.GetSize();
139 Int_t i = 0;
140 TArrayD b(n+1);
141 b[0] = fLow;
142 for (Int_t j = 0; j < l; j++) {
143 for (Int_t k = 0; k < fN[j]; k++) {
144 b[i+1] = b[i] + fD[j];
145 i++;
146 }
147 }
148 if (i != n) {
149 ::Warning("Axis", "Assigned bins, and number of bins not equal");
150 n = TMath::Min(i, n);
151 }
152 fAxis.Set(n, b.GetArray());
153 }
154 return fAxis;
155}
52b36573 156
c8b1a7db 157//====================================================================
bfab35d9 158AliForwardMultDists::AliForwardMultDists()
c8b1a7db 159 : AliBaseAODTask(),
bfab35d9 160 fBins(),
161 fSymmetric(0),
162 fNegative(0),
163 fPositive(0),
281a2bf8 164 fMCVertex(0),
165 fDiag(0),
bfab35d9 166 fForwardCache(0),
167 fCentralCache(0),
168 fMCCache(0),
c8b1a7db 169 fUsePhiAcc(true),
170 fIsSelected(false)
bfab35d9 171{}
172
173//____________________________________________________________________
174AliForwardMultDists::AliForwardMultDists(const char* name)
c8b1a7db 175 : AliBaseAODTask(name),
bfab35d9 176 fBins(),
177 fSymmetric(0),
178 fNegative(0),
179 fPositive(0),
281a2bf8 180 fMCVertex(0),
181 fDiag(0),
bfab35d9 182 fForwardCache(0),
183 fCentralCache(0),
184 fMCCache(0),
c8b1a7db 185 fUsePhiAcc(true),
186 fIsSelected(false)
bfab35d9 187{
188 /**
189 * User constructor
190 *
191 * @param name Name of the task
192 */
bfab35d9 193}
194
c8b1a7db 195
bfab35d9 196//____________________________________________________________________
c8b1a7db 197Bool_t
198AliForwardMultDists::Book()
bfab35d9 199{
c8b1a7db 200 /**
201 * Create output objects - called at start of job in slave
202 *
203 */
204 return true;
bfab35d9 205}
bfab35d9 206//____________________________________________________________________
c8b1a7db 207Bool_t
208AliForwardMultDists::PreData()
bfab35d9 209{
210 /**
c8b1a7db 211 * Set-up internal structures on first seen event
bfab35d9 212 *
c8b1a7db 213 * @param hist Basic histogram template from AOD object
bfab35d9 214 */
c8b1a7db 215 AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
216 AliAODForwardMult* forward = GetForward(*aod);
217 const TH2D& hist = forward->GetHistogram();
218 Bool_t useMC = GetPrimary(*aod) != 0;
bfab35d9 219
c8b1a7db 220 fSums->Add(AliForwardUtil::MakeParameter("minIpZ", fMinIpZ));
221 fSums->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
bfab35d9 222
c8b1a7db 223 TAxis* xaxis = hist.GetXaxis();
224 if (!xaxis->GetXbins() || xaxis->GetXbins()->GetSize() <= 0) {
225 fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
226 xaxis->GetNbins(), xaxis->GetXmin(),
227 xaxis->GetXmax());
228 fCentralCache = new TH1D("centralCache", "Projection of Central data",
229 xaxis->GetNbins(), xaxis->GetXmin(),
230 xaxis->GetXmax());
231 if (useMC)
232 fMCCache = new TH1D("mcCache", "Projection of Mc data",
233 xaxis->GetNbins(), xaxis->GetXmin(),
234 xaxis->GetXmax());
235 }
236 else {
237 fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
238 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
239 fCentralCache = new TH1D("centralCache", "Projection of Central data",
240 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
241 if (useMC)
242 fMCCache = new TH1D("mcCache", "Projection of Mc data",
243 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
244 }
245 fForwardCache->SetDirectory(0);
246 fForwardCache->GetXaxis()->SetTitle("#eta");
247 fForwardCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
248 fForwardCache->Sumw2();
249 fCentralCache->SetDirectory(0);
250 fCentralCache->GetXaxis()->SetTitle("#eta");
251 fCentralCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
252 fCentralCache->Sumw2();
bfab35d9 253
c8b1a7db 254 if (useMC) {
255 fMCCache->SetDirectory(0);
256 fMCCache->GetXaxis()->SetTitle("#eta");
257 fMCCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
258 fMCCache->Sumw2();
bfab35d9 259
c8b1a7db 260 fMCVertex = static_cast<TH1*>(fVertex->Clone("mcVertex"));
261 fMCVertex->SetTitle("Vertex distribution from MC");
262 fMCVertex->SetDirectory(0);
263 fMCVertex->SetFillColor(kBlue+2);
264 fSums->Add(fMCVertex);
265 }
bfab35d9 266
c8b1a7db 267 UShort_t xBase = kTrigger-1;
268 UShort_t yBase = kAnalysis-1;
269 fDiag = new TH2I("diagnostics", "Event selection",
270 kTriggerVertex-xBase, kTrigger-.5, kTriggerVertex+.5,
271 kTriggerVertex-yBase, kAnalysis-.5, kTriggerVertex+.5);
272 fDiag->SetDirectory(0);
273 fDiag->GetXaxis()->SetTitle("Selection");
274 fDiag->GetXaxis()->SetBinLabel(kTrigger -xBase, "Trigger");
275 fDiag->GetXaxis()->SetBinLabel(kVertex -xBase, "Vertex");
276 fDiag->GetXaxis()->SetBinLabel(kTriggerVertex-xBase, "Trigger+Vertex");
277 fDiag->GetYaxis()->SetTitle("Type/MC selection");
278 fDiag->GetYaxis()->SetBinLabel(kAnalysis -yBase, "Analysis");
279 fDiag->GetYaxis()->SetBinLabel(kMC -yBase, "MC");
280 fDiag->GetYaxis()->SetBinLabel(kTrigger -yBase, "MC Trigger");
281 fDiag->GetYaxis()->SetBinLabel(kVertex -yBase, "MC Vertex");
282 fDiag->GetYaxis()->SetBinLabel(kTriggerVertex-yBase, "MC Trigger+Vertex");
283 fDiag->SetMarkerSize(3);
284 fDiag->SetMarkerColor(kWhite);
285 fSums->Add(fDiag);
286
287 TIter next(&fBins);
288 EtaBin* bin = 0;
289 while ((bin = static_cast<EtaBin*>(next()))) {
290 bin->SetupForData(fSums, hist, useMC);
291 }
292 return true;
bfab35d9 293}
294//____________________________________________________________________
c8b1a7db 295Bool_t
296AliForwardMultDists::Event(AliAODEvent& aod)
bfab35d9 297{
298 /**
299 * Analyse a single event
300 *
301 * @param option Not used
302 */
c8b1a7db 303 AliAODForwardMult* forward = GetForward(aod, false, true);
304 if (!forward) return false;
bfab35d9 305
c8b1a7db 306 AliAODCentralMult* central = GetCentral(aod, false, true);
307 if (!central) return false;
bfab35d9 308
c8b1a7db 309 TH2* primary = GetPrimary(aod);
bfab35d9 310
311 const TH2& forwardData = forward->GetHistogram();
312 const TH2& centralData = central->GetHistogram();
313
c8b1a7db 314 Double_t vz = forward->GetIpZ();
315 Bool_t trg = forward->IsTriggerBits(fTriggerMask);
316 Bool_t vtx = (vz <= fMaxIpZ && vz >= fMinIpZ);
317 Bool_t ok = true; // Should bins process this event?
318 Bool_t mcTrg = false;
319 Bool_t mcVtx = false;
320 // If we have MC inpunt
281a2bf8 321 if (primary) {
322 // For MC, we need to check if we should process the event for
323 // MC information or not.
324 if ((fTriggerMask & (AliAODForwardMult::kV0AND|AliAODForwardMult::kNSD))
325 && !forward->IsTriggerBits(AliAODForwardMult::kMCNSD))
326 // Bail out if this event is not MC NSD event
327 ok = false;
c8b1a7db 328 else
329 mcTrg = true;
281a2bf8 330 Double_t mcVz = primary->GetBinContent(0,0);
331 fMCVertex->Fill(mcVz);
332 if (mcVz > fMaxIpZ || mcVz < fMinIpZ)
333 // Bail out if this event was not in the valid range
334 ok = false;
335 else
336 mcVtx = true;
337 }
c8b1a7db 338 // Fill diagnostics
281a2bf8 339 if (trg) {
340 fDiag->Fill(kTrigger, kAnalysis);
341 if (mcTrg) fDiag->Fill(kTrigger, kTrigger);
342 if (mcVtx) fDiag->Fill(kTrigger, kVertex);
343 if (mcTrg && mcVtx) fDiag->Fill(kTrigger, kTriggerVertex);
344 }
345 if (vtx) {
346 fDiag->Fill(kVertex, kAnalysis);
347 if (mcTrg) fDiag->Fill(kVertex, kTrigger);
348 if (mcVtx) fDiag->Fill(kVertex, kVertex);
349 if (mcTrg && mcVtx) fDiag->Fill(kVertex, kTriggerVertex);
350 }
351 if (trg && vtx) {
352 fDiag->Fill(kTriggerVertex, kAnalysis);
353 if (mcTrg) fDiag->Fill(kTriggerVertex, kTrigger);
354 if (mcVtx) fDiag->Fill(kTriggerVertex, kVertex);
355 if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kTriggerVertex);
356 }
357 if (primary) {
358 if (mcTrg) fDiag->Fill(kTrigger, kMC);
359 if (mcVtx) fDiag->Fill(kVertex, kMC);
360 if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kMC);
361 }
c8b1a7db 362 if (!fIsSelected && !ok) {
281a2bf8 363 // Printf("Event is neither accepted by analysis nor by MC");
c8b1a7db 364 return false;
281a2bf8 365 }
366
bfab35d9 367 // forward->Print();
368 // Info("", "Event vz=%f", forward->GetIpZ());
369
370 ProjectX(forwardData, *fForwardCache, fUsePhiAcc);
371 ProjectX(centralData, *fCentralCache, fUsePhiAcc);
372 ProjectX(primary, fMCCache);
373
374 TIter next(&fBins);
375 EtaBin* bin = 0;
376 while ((bin = static_cast<EtaBin*>(next()))) {
377 bin->Process(*fForwardCache, *fCentralCache,
281a2bf8 378 forwardData, centralData,
c8b1a7db 379 fIsSelected, fMCCache);
bfab35d9 380 }
381
c8b1a7db 382 return true;
bfab35d9 383}
bfab35d9 384
c8b1a7db 385//_____________________________________________________________________
386Bool_t
387AliForwardMultDists::CheckEvent(const AliAODForwardMult& fwd)
388{
389 fIsSelected = AliBaseAODTask::CheckEvent(fwd);
390 // We always return true, so that we can process MC truth;
391 return true;
bfab35d9 392}
c8b1a7db 393
bfab35d9 394//____________________________________________________________________
c8b1a7db 395Bool_t
396AliForwardMultDists::Finalize()
bfab35d9 397{
398 /**
399 * Called at the end of the final processing of the job on the
400 * full data set (merged data)
401 *
402 *
403 * @param option Not used
404 */
c8b1a7db 405 fResults->Add(fSums->FindObject("triggers")->Clone());
406 fResults->Add(fSums->FindObject("status")->Clone());
407 fResults->Add(fSums->FindObject("diagnostics")->Clone());
bfab35d9 408
409 THStack* sym = new THStack("all", "All distributions");
410 THStack* pos = new THStack("all", "All distributions");
411 THStack* neg = new THStack("all", "All distributions");
412 THStack* oth = new THStack("all", "All distributions");
413 THStack* sta = 0;
414 EtaBin* bin = 0;
415 TIter next(&fBins);
416 while ((bin = static_cast<EtaBin*>(next()))) {
c8b1a7db 417 bin->Terminate(fSums, fResults);
bfab35d9 418
419 sta = oth;
420 if (bin->IsSymmetric()) sta = sym;
421 else if (bin->IsNegative()) sta = neg;
422 else if (bin->IsPositive()) sta = pos;
423
281a2bf8 424 TH1* hh[] = { bin->fSum, bin->fTruth, bin->fTruthAccepted, 0 };
bfab35d9 425 TH1** ph = hh;
426
427 Int_t idx = sta->GetHists() ? sta->GetHists()->GetEntries() : 0;
428 if (bin->fTruth) idx /= 2;
429
430 Int_t mkrBits = GetIndexMarker(idx);
431 while (*ph) {
432 Int_t factor = TMath::Power(10, idx);
433 TH1* h = static_cast<TH1*>((*ph)->Clone());
434 h->SetDirectory(0);
435 h->Scale(factor);
436 h->SetTitle(Form("%s (#times%d)", h->GetTitle(), Int_t(factor)));
437 h->SetMarkerStyle(GetMarkerStyle(mkrBits));
438 mkrBits ^= kHollow;
439
440 sta->Add(h, "p");
441 ph++;
442 }
443 }
c8b1a7db 444 TList* lsym = static_cast<TList*>(fResults->FindObject("symmetric"));
445 TList* lneg = static_cast<TList*>(fResults->FindObject("negative"));
446 TList* lpos = static_cast<TList*>(fResults->FindObject("positive"));
447 TList* loth = static_cast<TList*>(fResults->FindObject("other"));
bfab35d9 448
449 if (lsym) lsym->Add(sym);
450 if (lneg) lneg->Add(neg);
451 if (lpos) lpos->Add(pos);
452 if (loth) loth->Add(oth);
453
c8b1a7db 454 return true;
bfab35d9 455}
456
bfab35d9 457//____________________________________________________________________
458void
459AliForwardMultDists::ProjectX(const TH2& input, TH1& cache, Bool_t usePhiAcc)
460{
461 /**
462 * Project a 2D histogram into a 1D histogram taking care to use
c8b1a7db 463 * either the @f$\phi@f$ acceptance stored in the overflow bins, or
bfab35d9 464 * the @f$\eta@f$ coverage stored in the underflow bins.
465 *
466 * @param input 2D histogram to project
467 * @param cache 1D histogram to project into
c8b1a7db 468 * @param usePhiAcc If true, use the @f$\phi@f$ acceptance stored in
bfab35d9 469 * the overflow bins, or if false the @f$\eta@f$ coverage stored in
470 * the underflow bins.
471 */
472 cache.Reset();
473
474 Int_t nPhi = input.GetNbinsY();
475 Int_t nEta = input.GetNbinsX();
476
477 for (Int_t iEta = 1; iEta <= nEta; iEta++) {
478 Double_t phiAcc = input.GetBinContent(iEta, nPhi+1);
479 Double_t etaCov = input.GetBinContent(iEta, 0);
480 Double_t factor = usePhiAcc ? phiAcc : etaCov;
481
482 if (factor < 1e-3) continue;
483 Double_t sum = 0;
484 Double_t e2sum = 0;
485 for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
486 Double_t c = input.GetBinContent(iEta, iPhi);
487 Double_t e = input.GetBinError(iEta, iPhi);
488
489 sum += c;
490 e2sum += e * e;
491 }
492 cache.SetBinContent(iEta, factor * sum);
493 cache.SetBinError(iEta, factor * TMath::Sqrt(e2sum));
494 }
495}
496//____________________________________________________________________
497void
498AliForwardMultDists::ProjectX(const TH2* input, TH1* cache)
499{
500 /**
501 * Project on @f$\eta@f$ axis. If any of the pointers passed is
502 * zero, do nothing.
503 *
504 * @param input
505 * @param cache
506 */
507 if (!input || !cache) return;
508 cache->Reset();
509
510 Int_t nPhi = input->GetNbinsY();
511 Int_t nEta = input->GetNbinsX();
512
513 for (Int_t iEta = 1; iEta <= nEta; iEta++) {
514 Double_t sum = 0;
515 Double_t e2sum = 0;
516 for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
517 Double_t c = input->GetBinContent(iEta, iPhi);
518 Double_t e = input->GetBinError(iEta, iPhi);
519
520 sum += c;
521 e2sum += e * e;
522 }
523 cache->SetBinContent(iEta, sum);
524 cache->SetBinError(iEta, TMath::Sqrt(e2sum));
525 }
526}
527//____________________________________________________________________
528void
52b36573 529AliForwardMultDists::AddBin(const BinSpec& spec)
530{
531 AddBin(spec.fEtaMin, spec.fEtaMax, spec.Axis());
532}
533
534//____________________________________________________________________
535void
536AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax,
537 const TAxis& mAxis)
bfab35d9 538{
539 /**
540 * Add an @f$\eta@f$ bin
541 *
542 * @param etaLow Low cut on @f$\eta@f$
543 * @param etaMax High cut on @f$\eta@f$
544 */
545 // Symmetric bin
546 if (etaMax >= kInvalidEta) {
547 etaLow = -TMath::Abs(etaLow);
548 etaMax = +TMath::Abs(etaLow);
549 }
52b36573 550 EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
551 // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
552 fBins.Add(bin);
553}
554//____________________________________________________________________
555void
556AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax,
557 UShort_t nMax, UShort_t nDiv)
558{
559 /**
560 * Add an @f$\eta@f$ bin
561 *
562 * @param etaLow Low cut on @f$\eta@f$
563 * @param etaMax High cut on @f$\eta@f$
564 */
565 // Symmetric bin
566 if (etaMax >= kInvalidEta) {
567 etaLow = -TMath::Abs(etaLow);
568 etaMax = +TMath::Abs(etaLow);
569 }
570 TAxis mAxis((nMax+1)/nDiv, -.5, nMax+.5);
571 EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
bfab35d9 572 // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
573 fBins.Add(bin);
574}
c8b1a7db 575#define PFV(N,VALUE) \
576 do { \
577 AliForwardUtil::PrintName(N); \
578 std::cout << (VALUE) << std::endl; } while(false)
579#define PFB(N,FLAG) \
580 do { \
581 AliForwardUtil::PrintName(N); \
582 std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
583 } while(false)
584
bfab35d9 585//____________________________________________________________________
586void
c8b1a7db 587AliForwardMultDists::Print(Option_t* option) const
bfab35d9 588{
589 /**
590 * Print this task
591 *
592 * @param option Not used
593 */
c8b1a7db 594 AliBaseAODTask::Print(option);
595 gROOT->IncreaseDirLevel();
596 PFB("Use phi acceptance", fUsePhiAcc);
597 // AliForwardUtil::PrintName("Bins");
598
599 Int_t first = true;
bfab35d9 600 TIter next(&fBins);
601 EtaBin* bin = 0;
602 while ((bin = static_cast<EtaBin*>(next()))) {
c8b1a7db 603 PFV(first ? "Bins" : "", bin->GetName());
604 first = false;
bfab35d9 605 }
c8b1a7db 606 gROOT->DecreaseDirLevel();
bfab35d9 607}
608
609//====================================================================
610AliForwardMultDists::EtaBin::EtaBin()
611 : fName(""),
52b36573 612 fMAxis(1,0,0),
613 fTAxis(1,0,0),
bfab35d9 614 fMinEta(0),
615 fMaxEta(0),
616 fMinBin(0),
617 fMaxBin(0),
618 fSum(0),
619 fCorr(0),
620 fResponse(0),
621 fTruth(0),
281a2bf8 622 fTruthAccepted(0),
bfab35d9 623 fCoverage(0)
624{
625 /**
626 * I/O constructor
627 */
628}
629
630//____________________________________________________________________
52b36573 631AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta,
632 const TAxis& mAxis)
bfab35d9 633 : fName(""),
52b36573 634 fMAxis(1,0,0),
635 fTAxis(1,0,0),
bfab35d9 636 fMinEta(minEta),
637 fMaxEta(maxEta),
638 fMinBin(0),
639 fMaxBin(0),
640 fSum(0),
641 fCorr(0),
642 fResponse(0),
643 fTruth(0),
281a2bf8 644 fTruthAccepted(0),
bfab35d9 645 fCoverage(0)
646{
647 /**
648 * User constructor
649 *
650 * @param minEta Least @f$\eta@f$ to consider
651 * @param maxEta Largest @f$\eta@f$ to consider
652 */
653 fName = TString::Format("%+05.2f_%+05.2f", fMinEta, fMaxEta);
654 fName.ReplaceAll("-", "m");
655 fName.ReplaceAll("+", "p");
656 fName.ReplaceAll(".", "d");
52b36573 657
658 // Copy to other our object
659 mAxis.Copy(fMAxis);
660
661 if (mAxis.GetXbins() && mAxis.GetXbins()->GetArray()) {
662 const TArrayD& mA = *(mAxis.GetXbins());
663 TArrayD tA(mA.GetSize());
664 Int_t j = 0;
665 Double_t min = mA[0];
666 tA[0] = min;
667 for (Int_t i = 1; i < mA.GetSize(); i++) {
668 Double_t d = mA[i] - min;
669 if (d < 1)
670 // Not full integer bin
671 continue;
672 tA[j+1] = tA[j] + Int_t(d);
673 min = tA[j+1];
674 j++;
675 }
676 fTAxis.Set(j, tA.GetArray());
677 }
678 else {
679 // Rounded down maximum and minimum
680 Int_t max = mAxis.GetXmax();
681 Int_t min = mAxis.GetXmin();
682 fTAxis.Set((max-min)+1, min-.5, max+.5);
683 }
bfab35d9 684}
685//____________________________________________________________________
686AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o)
687 : TObject(o),
688 fName(o.fName),
52b36573 689 fMAxis(o.fMAxis),
690 fTAxis(o.fTAxis),
bfab35d9 691 fMinEta(o.fMinEta),
692 fMaxEta(o.fMaxEta),
693 fMinBin(o.fMinBin),
694 fMaxBin(o.fMaxBin),
695 fSum(o.fSum),
696 fCorr(o.fCorr),
697 fResponse(o.fResponse),
698 fTruth(o.fTruth),
281a2bf8 699 fTruthAccepted(o.fTruthAccepted),
bfab35d9 700 fCoverage(o.fCoverage)
701{}
702//____________________________________________________________________
703AliForwardMultDists::EtaBin&
704AliForwardMultDists::EtaBin::operator=(const EtaBin& o)
705{
706 if (&o == this) return *this;
707
281a2bf8 708 fName = o.fName;
52b36573 709 fMAxis = o.fMAxis;
710 fTAxis = o.fTAxis;
281a2bf8 711 fMinEta = o.fMinEta;
712 fMaxEta = o.fMaxEta;
713 fMinBin = o.fMinBin;
714 fMaxBin = o.fMaxBin;
715 fSum = o.fSum;
716 fCorr = o.fCorr;
717 fResponse = o.fResponse;
718 fTruth = o.fTruth;
719 fTruthAccepted = o.fTruthAccepted;
720 fCoverage = o.fCoverage;
bfab35d9 721
722 return *this;
723}
724//____________________________________________________________________
725Bool_t
726AliForwardMultDists::EtaBin::IsSymmetric() const
727{
728 return TMath::Abs(fMaxEta + fMinEta) < 1e-6;
729}
730//____________________________________________________________________
731Bool_t
732AliForwardMultDists::EtaBin::IsNegative() const
733{
734 return TMath::Abs(fMaxEta) < 1e-6 && fMinEta < 0;
735}
736//____________________________________________________________________
737Bool_t
738AliForwardMultDists::EtaBin::IsPositive() const
739{
740 return TMath::Abs(fMinEta) < 1e-6 && fMaxEta > 0;
741}
742//____________________________________________________________________
743const char*
744AliForwardMultDists::EtaBin::ParentName() const
745{
746 if (IsSymmetric()) return "symmetric";
747 else if (IsNegative()) return "negative";
748 else if (IsPositive()) return "positive";
749 return "other";
750}
751//____________________________________________________________________
752TList*
753AliForwardMultDists::EtaBin::FindParent(TList* l, Bool_t create) const
754{
755 const char* parent = ParentName();
756 TObject* op = l->FindObject(parent);
757
758 if (op) return static_cast<TList*>(op);
759 if (!create) return 0;
760
761 // Info("FindParent", "Parent %s not found in %s, creating and adding",
762 // parent, l->GetName());
763 TList* p = new TList;
764 p->SetName(parent);
765 p->SetOwner();
766 l->Add(p);
767
768 TList* lowEdges = new TList;
769 lowEdges->SetName("lowEdges");
770 lowEdges->SetOwner();
771 p->Add(lowEdges);
772
773 TList* highEdges = new TList;
774 highEdges->SetName("highEdges");
775 highEdges->SetOwner();
776 p->Add(highEdges);
777
778 return p;
779}
780
52b36573 781//____________________________________________________________________
782TH1*
783AliForwardMultDists::EtaBin::CreateH1(const char* name,
784 const char* title,
785 const TAxis& xAxis)
786{
787 TH1* ret = 0;
788
789 if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray())
790 ret = new TH1D(name,title,xAxis.GetNbins(), xAxis.GetXbins()->GetArray());
791 else
792 ret = new TH1D(name,title,xAxis.GetNbins(),xAxis.GetXmin(),xAxis.GetXmax());
793 return ret;
794}
795//____________________________________________________________________
796TH2*
797AliForwardMultDists::EtaBin::CreateH2(const char* name,
798 const char* title,
799 const TAxis& xAxis,
800 const TAxis& yAxis)
801{
802 TH2* ret = 0;
803
804 if (xAxis.GetXbins() && xAxis.GetXbins()->GetArray()) {
805 if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray())
806 // Variable variable
807 ret = new TH2D(name,title,
808 xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
809 yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
810 else
811 // variable fixed
812 ret = new TH2D(name,title,
813 xAxis.GetNbins(), xAxis.GetXbins()->GetArray(),
814 yAxis.GetNbins(),yAxis.GetXmin(),yAxis.GetXmax());
815 }
816 else {
817 if (yAxis.GetXbins() && yAxis.GetXbins()->GetArray())
818 // fixed variable
819 ret = new TH2D(name,title,
820 xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
821 yAxis.GetNbins(), yAxis.GetXbins()->GetArray());
822 else
823 // fixed fixed
824 ret = new TH2D(name,title,
825 xAxis.GetNbins(), xAxis.GetXmin(), xAxis.GetXmax(),
826 yAxis.GetNbins(), yAxis.GetXmin(), yAxis.GetXmax());
827 }
828 return ret;
829}
830
bfab35d9 831//____________________________________________________________________
832void
833AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist,
281a2bf8 834 Bool_t useMC)
bfab35d9 835{
836 /**
837 * Set-up internal structures on first event.
838 *
839 * @param list List to add information to
840 * @param hist Template histogram
841 * @param max Maximum number of particles
842 * @param useMC Whether to set-up for MC input
843 */
844 TList* p = FindParent(list, true);
845 TList* l = new TList;
846 l->SetName(GetName());
847 l->SetOwner();
848 p->Add(l);
849 TList* le = static_cast<TList*>(p->FindObject("lowEdges"));
850 TList* he = static_cast<TList*>(p->FindObject("highEdges"));
851 if (!le || !he) {
852 AliError("Failed to get bin edge lists from parent");
853 return;
854 }
855 else {
856 Int_t n = le->GetEntries();
857 TParameter<double>* lp =
858 new TParameter<double>(Form("minEta%02d", n), fMinEta);
859 TParameter<double>* hp =
860 new TParameter<double>(Form("maxEta%02d", n), fMaxEta);
861 lp->SetMergeMode('f');
862 hp->SetMergeMode('f');
863 le->Add(lp);
864 he->Add(hp);
865 }
866
867 fMinBin = hist.GetXaxis()->FindBin(fMinEta);
868 fMaxBin = hist.GetXaxis()->FindBin(fMaxEta-.00001);
281a2bf8 869
52b36573 870 TString t(Form("%+5.2f<#eta<%+5.2f", fMinEta, fMaxEta));
871 fSum = CreateH1("rawDist",Form("Raw P(#it{N}_{ch}) in %s",t.Data()),fMAxis);
bfab35d9 872 fSum->SetDirectory(0);
873 fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
874 fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
875 fSum->SetFillColor(kRed+1);
876 fSum->SetFillStyle(0);
877 // fSum->SetOption("hist e");
878 fSum->SetMarkerStyle(20);
879 fSum->SetMarkerColor(kRed+1);
880 fSum->SetLineColor(kBlack);
881 fSum->Sumw2();
882 l->Add(fSum);
883
52b36573 884 fCorr = CreateH2("corr",Form("C_{SPD,FMD} in %s", t.Data()),fTAxis,fTAxis);
bfab35d9 885 fCorr->SetDirectory(0);
886 fCorr->GetXaxis()->SetTitle("Forward");
887 fCorr->GetYaxis()->SetTitle("Central");
888 fCorr->SetOption("colz");
889 l->Add(fCorr);
890
891 fCoverage = new TH1D("coverage", "Fraction of bins with coverage",
892 101, -.5, 100.5);
893 fCoverage->SetDirectory(0);
894 fCoverage->SetXTitle("Fraction of bins [%]");
895 fCoverage->SetFillStyle(3001);
896 fCoverage->SetFillColor(kGreen+1);
897 l->Add(fCoverage);
898
899 if (!useMC) return;
52b36573 900 fResponse = CreateH2("response", Form("Reponse matrix in %s", t.Data()),
901 fMAxis, fTAxis);
bfab35d9 902 fResponse->SetDirectory(0);
281a2bf8 903 fResponse->SetYTitle("MC");
904 fResponse->SetXTitle("Signal");
bfab35d9 905 fResponse->SetOption("colz");
906 l->Add(fResponse);
907
52b36573 908 fTruth = CreateH1("truth",Form("True P(#it{N}_{ch}) in %s",t.Data()),fTAxis);
281a2bf8 909 fTruth->SetXTitle(fSum->GetXaxis()->GetTitle());
910 fTruth->SetYTitle(fSum->GetYaxis()->GetTitle());
bfab35d9 911 fTruth->SetLineColor(kBlack);
912 fTruth->SetFillColor(kBlue+1);
913 fTruth->SetFillStyle(0);
914 fTruth->SetDirectory(0);
915 /// fTruth->SetOption("");
916 fTruth->SetMarkerColor(kBlue+1);
917 fTruth->SetMarkerStyle(24);
281a2bf8 918 fTruth->Sumw2();
bfab35d9 919 l->Add(fTruth);
281a2bf8 920
921 fTruthAccepted = static_cast<TH1D*>(fTruth->Clone("truthAccepted"));
922 fTruthAccepted->SetTitle(Form("True (accepted) P(#it{N}_{ch}) in %s",
52b36573 923 t.Data()));
281a2bf8 924 fTruthAccepted->SetLineColor(kGray+2);
925 fTruthAccepted->SetFillColor(kOrange+2);
926 fTruthAccepted->SetDirectory(0);
927 /// fTruth->SetOption("");
928 fTruthAccepted->SetMarkerColor(kOrange+2);
929 l->Add(fTruthAccepted);
bfab35d9 930}
931//____________________________________________________________________
932void
933AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
934 const TH1& sumCentral,
935 const TH2& forward,
936 const TH2& central,
281a2bf8 937 Bool_t accepted,
bfab35d9 938 const TH1* mc)
939{
940 /**
941 * Process a single event
942 *
943 * @param sumForward Projection of forward data
944 * @param sumCentral Projection of the central data
945 * @param forward The original forward data
946 * @param central The original central data
947 */
281a2bf8 948 if (!mc && !accepted)
949 // If we're not looking at MC data, and the event wasn't selected,
950 // we bail out - this check is already done, but we do it again
951 // for safety.
952 return;
953
bfab35d9 954 Double_t sum = 0;
955 Double_t e2sum = 0;
956 Int_t covered = 0;
957 Double_t fsum = -1;
958 Double_t csum = -1;
959 Double_t mcsum = 0;
960 Double_t mce2sum = 0;
961 for (Int_t iEta = fMinBin; iEta <= fMaxBin; iEta++) {
281a2bf8 962 if (mc) {
963 Double_t cM = mc->GetBinContent(iEta);
964 Double_t eM = mc->GetBinError(iEta);
965 mcsum += cM;
966 mce2sum += eM * eM;
967 }
968 if (!accepted)
969 // Event wasn't selected, but we still need to get the rest of
970 // the MC data.
971 continue;
972
bfab35d9 973 Double_t cF = sumForward.GetBinContent(iEta);
974 Double_t eF = sumForward.GetBinError(iEta);
975 Bool_t bF = forward.GetBinContent(iEta, 0) > 0;
976 Double_t cC = sumCentral.GetBinContent(iEta);
977 Double_t eC = sumCentral.GetBinError(iEta);
978 Bool_t bC = central.GetBinContent(iEta, 0) > 0;
bfab35d9 979 Double_t c = 0;
980 Double_t e = 0;
281a2bf8 981
bfab35d9 982 // If we have an overlap - as given by the eta-coverage,
983 // calculate the mean
984 if (bF & bC) {
985 c = (cF + cC) / 2;
986 e = TMath::Sqrt(eF*eF + eC*eC);
987 covered++;
988 }
989 // Else, if we have coverage from forward, use that
990 else if (bF) {
991 c = cF;
992 e = eF;
993 covered++;
994 }
995 // Else, if we have covrage from central, use that
996 else if (bC) {
997 c = cC;
998 e = eC;
999 covered++;
1000 }
1001 // Else, we have incomplete coverage
1002
1003 if (bF) {
1004 if (fsum < 0) fsum = 0;
1005 fsum += cF;
1006 }
1007 if (bC) {
1008 if (csum < 0) csum = 0;
1009 csum += cC;
1010 }
1011
1012 sum += c;
1013 e2sum += e*e;
1014 }
1015
281a2bf8 1016 if (accepted) {
1017 // Only update the histograms if the event was triggered.
1018 // Fill with weight
1019 Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
1020 fSum->Fill(sum, w);
52b36573 1021 fCorr->Fill((fsum<0?0:fsum), (csum<0?0:csum));
281a2bf8 1022
1023 Int_t nTotal = fMaxBin - fMinBin + 1;
1024 fCoverage->Fill(100*float(covered) / nTotal);
1025 }
bfab35d9 1026
1027 if (mc) {
281a2bf8 1028 Double_t w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
bfab35d9 1029 if (fTruth) {
bfab35d9 1030 fTruth->Fill(mcsum, w);
1031 }
281a2bf8 1032 if (accepted) {
1033 if (fResponse)
1034 // Only fill response matrix for accepted events
1035 fResponse->Fill(sum, mcsum);
1036 if (fTruthAccepted)
1037 fTruthAccepted->Fill(mcsum, w);
1038 }
bfab35d9 1039 }
1040}
1041//____________________________________________________________________
1042void
52b36573 1043AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out)
bfab35d9 1044{
1045 /**
1046 * Called at the end of the final processing of the job on the
1047 * full data set (merged data)
1048 *
1049 * @param in Input list
1050 * @param out Output list
1051 * @param maxN Maximum number of @f$N_{ch}@f$ to consider
1052 */
1053 TList* ip = FindParent(in, false);
1054 if (!ip) {
1055 AliErrorF("Parent folder %s not found in input", ParentName());
1056 return;
1057 }
1058
1059 TList* i = dynamic_cast<TList*>(ip->FindObject(GetName()));
1060 if (!i) {
1061 AliErrorF("List %s not found in input", GetName());
1062 return;
1063 }
1064
1065 TList* op = FindParent(out, true);
1066 TList* o = static_cast<TList*>(i->Clone());
1067 o->SetOwner();
1068 op->Add(o);
1069
281a2bf8 1070 fSum = static_cast<TH1*>(o->FindObject("rawDist"));
1071 fTruth = static_cast<TH1*>(o->FindObject("truth"));
1072 fTruthAccepted = static_cast<TH1*>(o->FindObject("truthAccepted"));
1073
1074 TH1* hists[] = { fSum, fTruth, fTruthAccepted, 0 };
bfab35d9 1075 TH1** phist = hists;
1076 while (*phist) {
1077 TH1* h = *phist;
1078 if (h) {
52b36573 1079 Int_t maxN = h->GetNbinsX();
bfab35d9 1080 Double_t intg = h->Integral(1, maxN);
281a2bf8 1081 h->Scale(1. / intg, "width");
bfab35d9 1082 }
1083 phist++;
1084 }
281a2bf8 1085
1086 if (fTruth && fTruthAccepted) {
1087 TH1* trgVtx = static_cast<TH1*>(fTruthAccepted->Clone("triggerVertex"));
1088 TString tit(fTruth->GetTitle());
1089 tit.ReplaceAll("True P(#it{N}_{ch})", "C_{trigger,vertex}");
1090 trgVtx->SetTitle(tit);
1091 trgVtx->SetYTitle("P_{MC}(#it{N}_{ch})/P_{MC,accepted}(#it{N}_{ch})");
1092 trgVtx->Divide(fTruth);
1093 trgVtx->SetDirectory(0);
1094 o->Add(trgVtx);
1095 }
1096
bfab35d9 1097}
1098//====================================================================
1099//
1100// EOF
1101//