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