]> 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
Updates
[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 }
329a685b 54#if 0
c8b1a7db 55 /**
56 * Get the marker option bits from a ROOT style
57 *
58 * @param style ROOT style
59 *
60 * @return Pattern of marker options
61 */
62 UShort_t GetMarkerBits(Int_t style)
63 {
64 UShort_t bits = 0;
65 switch (style) {
66 case 24: case 25: case 26: case 27: case 28: case 30: case 32:
67 bits |= kHollow; break;
68 }
69 switch (style) {
70 case 20: case 24: bits |= kCircle; break;
71 case 21: case 25: bits |= kSquare; break;
72 case 22: case 26: bits |= kUpTriangle; break;
73 case 23: case 32: bits |= kDownTriangle; break;
74 case 27: case 33: bits |= kDiamond; break;
75 case 28: case 34: bits |= kCross; break;
76 case 29: case 30: bits |= kStar; break;
77 }
78 return bits;
79 }
329a685b 80#endif
c8b1a7db 81 static Int_t GetIndexMarker(Int_t idx)
82 {
83 const UShort_t nMarkers = 7;
84 UShort_t markers[] = {
85 kCircle,
86 kSquare,
87 kUpTriangle,
88 kDownTriangle,
89 kDiamond,
90 kCross,
91 kStar
92 };
93
94 return markers[idx % nMarkers];
95 }
329a685b 96#if 0
c8b1a7db 97 /**
98 * Flip the 'hollow-ness' of a marker
99 *
100 * @param style ROOT Style
101 *
102 * @return ROOT style
103 */
104 Int_t FlipHollowStyle(Int_t style)
105 {
106 UShort_t bits = GetMarkerBits(style);
107 Int_t ret = GetMarkerStyle(bits ^ kHollow);
108 return ret;
109 }
329a685b 110#endif
c8b1a7db 111}
112
113//====================================================================
52b36573 114AliForwardMultDists::BinSpec::BinSpec(Double_t etaMin,
115 Double_t etaMax,
116 Double_t nchLow)
117 : fEtaMin(etaMin), fEtaMax(etaMax), fLow(nchLow), fN(0), fD(0), fAxis(1,0,1)
118{
119}
120//____________________________________________________________________
121void
122AliForwardMultDists::BinSpec::Push(UShort_t n, Double_t d)
123{
124 Int_t l = fN.GetSize();
125 fN.Set(l+1); // Not terribly efficient, but that's ok because we do
126 fD.Set(l+1); // this infrequently
127 fN[l] = n;
128 fD[l] = d;
129}
130//____________________________________________________________________
131const TAxis&
132AliForwardMultDists::BinSpec::Axis() const
133{
134 if (fAxis.GetNbins() > 1) return fAxis;
135 if (fN.GetSize() <= 0) return fAxis;
136 if (fN.GetSize() == 1) {
137 // Exactly one spec,
138 fAxis.Set(fN[0], fLow, fN[0]*fD[0]+fLow);
139 }
140 else {
77f97e3f 141 Int_t n = Int_t(fN.GetSum());
52b36573 142 Int_t l = fN.GetSize();
143 Int_t i = 0;
144 TArrayD b(n+1);
145 b[0] = fLow;
146 for (Int_t j = 0; j < l; j++) {
147 for (Int_t k = 0; k < fN[j]; k++) {
148 b[i+1] = b[i] + fD[j];
149 i++;
150 }
151 }
152 if (i != n) {
153 ::Warning("Axis", "Assigned bins, and number of bins not equal");
154 n = TMath::Min(i, n);
155 }
156 fAxis.Set(n, b.GetArray());
157 }
158 return fAxis;
159}
52b36573 160
c8b1a7db 161//====================================================================
bfab35d9 162AliForwardMultDists::AliForwardMultDists()
c8b1a7db 163 : AliBaseAODTask(),
bfab35d9 164 fBins(),
165 fSymmetric(0),
166 fNegative(0),
167 fPositive(0),
281a2bf8 168 fMCVertex(0),
169 fDiag(0),
bfab35d9 170 fForwardCache(0),
171 fCentralCache(0),
172 fMCCache(0),
c8b1a7db 173 fUsePhiAcc(true),
174 fIsSelected(false)
bfab35d9 175{}
176
177//____________________________________________________________________
178AliForwardMultDists::AliForwardMultDists(const char* name)
bb64e168 179 : AliBaseAODTask(name,"AliForwardMultDists"),
bfab35d9 180 fBins(),
181 fSymmetric(0),
182 fNegative(0),
183 fPositive(0),
281a2bf8 184 fMCVertex(0),
185 fDiag(0),
bfab35d9 186 fForwardCache(0),
187 fCentralCache(0),
188 fMCCache(0),
c8b1a7db 189 fUsePhiAcc(true),
190 fIsSelected(false)
bfab35d9 191{
67a4bb96 192 /*
bfab35d9 193 * User constructor
194 *
195 * @param name Name of the task
196 */
bfab35d9 197}
198
c8b1a7db 199
bfab35d9 200//____________________________________________________________________
c8b1a7db 201Bool_t
202AliForwardMultDists::Book()
bfab35d9 203{
67a4bb96 204 /*
c8b1a7db 205 * Create output objects - called at start of job in slave
206 *
207 */
208 return true;
bfab35d9 209}
bfab35d9 210//____________________________________________________________________
c8b1a7db 211Bool_t
212AliForwardMultDists::PreData()
bfab35d9 213{
67a4bb96 214 /*
c8b1a7db 215 * Set-up internal structures on first seen event
bfab35d9 216 *
217 */
c8b1a7db 218 AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
219 AliAODForwardMult* forward = GetForward(*aod);
220 const TH2D& hist = forward->GetHistogram();
221 Bool_t useMC = GetPrimary(*aod) != 0;
bfab35d9 222
c8b1a7db 223 fSums->Add(AliForwardUtil::MakeParameter("minIpZ", fMinIpZ));
224 fSums->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
bfab35d9 225
700a5a26 226 const TAxis* xaxis = hist.GetXaxis();
c8b1a7db 227 if (!xaxis->GetXbins() || xaxis->GetXbins()->GetSize() <= 0) {
228 fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
229 xaxis->GetNbins(), xaxis->GetXmin(),
230 xaxis->GetXmax());
231 fCentralCache = new TH1D("centralCache", "Projection of Central data",
232 xaxis->GetNbins(), xaxis->GetXmin(),
233 xaxis->GetXmax());
234 if (useMC)
235 fMCCache = new TH1D("mcCache", "Projection of Mc data",
236 xaxis->GetNbins(), xaxis->GetXmin(),
237 xaxis->GetXmax());
238 }
239 else {
240 fForwardCache = new TH1D("forwardCache", "Projection of Forward data",
241 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
242 fCentralCache = new TH1D("centralCache", "Projection of Central data",
243 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
244 if (useMC)
245 fMCCache = new TH1D("mcCache", "Projection of Mc data",
246 xaxis->GetNbins(),xaxis->GetXbins()->GetArray());
247 }
248 fForwardCache->SetDirectory(0);
249 fForwardCache->GetXaxis()->SetTitle("#eta");
250 fForwardCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
251 fForwardCache->Sumw2();
252 fCentralCache->SetDirectory(0);
253 fCentralCache->GetXaxis()->SetTitle("#eta");
254 fCentralCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
255 fCentralCache->Sumw2();
bfab35d9 256
c8b1a7db 257 if (useMC) {
258 fMCCache->SetDirectory(0);
259 fMCCache->GetXaxis()->SetTitle("#eta");
260 fMCCache->GetYaxis()->SetTitle("#sum#frac{d^{2}N}{d#etadphi}");
261 fMCCache->Sumw2();
bfab35d9 262
c8b1a7db 263 fMCVertex = static_cast<TH1*>(fVertex->Clone("mcVertex"));
264 fMCVertex->SetTitle("Vertex distribution from MC");
265 fMCVertex->SetDirectory(0);
266 fMCVertex->SetFillColor(kBlue+2);
267 fSums->Add(fMCVertex);
268 }
bfab35d9 269
c8b1a7db 270 UShort_t xBase = kTrigger-1;
271 UShort_t yBase = kAnalysis-1;
272 fDiag = new TH2I("diagnostics", "Event selection",
273 kTriggerVertex-xBase, kTrigger-.5, kTriggerVertex+.5,
274 kTriggerVertex-yBase, kAnalysis-.5, kTriggerVertex+.5);
275 fDiag->SetDirectory(0);
276 fDiag->GetXaxis()->SetTitle("Selection");
277 fDiag->GetXaxis()->SetBinLabel(kTrigger -xBase, "Trigger");
278 fDiag->GetXaxis()->SetBinLabel(kVertex -xBase, "Vertex");
279 fDiag->GetXaxis()->SetBinLabel(kTriggerVertex-xBase, "Trigger+Vertex");
280 fDiag->GetYaxis()->SetTitle("Type/MC selection");
281 fDiag->GetYaxis()->SetBinLabel(kAnalysis -yBase, "Analysis");
282 fDiag->GetYaxis()->SetBinLabel(kMC -yBase, "MC");
283 fDiag->GetYaxis()->SetBinLabel(kTrigger -yBase, "MC Trigger");
284 fDiag->GetYaxis()->SetBinLabel(kVertex -yBase, "MC Vertex");
285 fDiag->GetYaxis()->SetBinLabel(kTriggerVertex-yBase, "MC Trigger+Vertex");
286 fDiag->SetMarkerSize(3);
287 fDiag->SetMarkerColor(kWhite);
288 fSums->Add(fDiag);
289
290 TIter next(&fBins);
291 EtaBin* bin = 0;
292 while ((bin = static_cast<EtaBin*>(next()))) {
293 bin->SetupForData(fSums, hist, useMC);
294 }
295 return true;
bfab35d9 296}
297//____________________________________________________________________
c8b1a7db 298Bool_t
299AliForwardMultDists::Event(AliAODEvent& aod)
bfab35d9 300{
67a4bb96 301 /*
bfab35d9 302 * Analyse a single event
303 *
67a4bb96 304 * @param aod Input event
bfab35d9 305 */
c8b1a7db 306 AliAODForwardMult* forward = GetForward(aod, false, true);
307 if (!forward) return false;
bfab35d9 308
c8b1a7db 309 AliAODCentralMult* central = GetCentral(aod, false, true);
310 if (!central) return false;
bfab35d9 311
c8b1a7db 312 TH2* primary = GetPrimary(aod);
bfab35d9 313
314 const TH2& forwardData = forward->GetHistogram();
315 const TH2& centralData = central->GetHistogram();
316
c8b1a7db 317 Double_t vz = forward->GetIpZ();
318 Bool_t trg = forward->IsTriggerBits(fTriggerMask);
319 Bool_t vtx = (vz <= fMaxIpZ && vz >= fMinIpZ);
320 Bool_t ok = true; // Should bins process this event?
321 Bool_t mcTrg = false;
322 Bool_t mcVtx = false;
323 // If we have MC inpunt
281a2bf8 324 if (primary) {
325 // For MC, we need to check if we should process the event for
326 // MC information or not.
327 if ((fTriggerMask & (AliAODForwardMult::kV0AND|AliAODForwardMult::kNSD))
328 && !forward->IsTriggerBits(AliAODForwardMult::kMCNSD))
329 // Bail out if this event is not MC NSD event
330 ok = false;
c8b1a7db 331 else
332 mcTrg = true;
281a2bf8 333 Double_t mcVz = primary->GetBinContent(0,0);
334 fMCVertex->Fill(mcVz);
335 if (mcVz > fMaxIpZ || mcVz < fMinIpZ)
336 // Bail out if this event was not in the valid range
337 ok = false;
338 else
339 mcVtx = true;
340 }
c8b1a7db 341 // Fill diagnostics
281a2bf8 342 if (trg) {
343 fDiag->Fill(kTrigger, kAnalysis);
344 if (mcTrg) fDiag->Fill(kTrigger, kTrigger);
345 if (mcVtx) fDiag->Fill(kTrigger, kVertex);
346 if (mcTrg && mcVtx) fDiag->Fill(kTrigger, kTriggerVertex);
347 }
348 if (vtx) {
349 fDiag->Fill(kVertex, kAnalysis);
350 if (mcTrg) fDiag->Fill(kVertex, kTrigger);
351 if (mcVtx) fDiag->Fill(kVertex, kVertex);
352 if (mcTrg && mcVtx) fDiag->Fill(kVertex, kTriggerVertex);
353 }
354 if (trg && vtx) {
355 fDiag->Fill(kTriggerVertex, kAnalysis);
356 if (mcTrg) fDiag->Fill(kTriggerVertex, kTrigger);
357 if (mcVtx) fDiag->Fill(kTriggerVertex, kVertex);
358 if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kTriggerVertex);
359 }
360 if (primary) {
361 if (mcTrg) fDiag->Fill(kTrigger, kMC);
362 if (mcVtx) fDiag->Fill(kVertex, kMC);
363 if (mcTrg && mcVtx) fDiag->Fill(kTriggerVertex, kMC);
364 }
c8b1a7db 365 if (!fIsSelected && !ok) {
281a2bf8 366 // Printf("Event is neither accepted by analysis nor by MC");
c8b1a7db 367 return false;
281a2bf8 368 }
369
bfab35d9 370 // forward->Print();
371 // Info("", "Event vz=%f", forward->GetIpZ());
372
373 ProjectX(forwardData, *fForwardCache, fUsePhiAcc);
374 ProjectX(centralData, *fCentralCache, fUsePhiAcc);
375 ProjectX(primary, fMCCache);
376
377 TIter next(&fBins);
378 EtaBin* bin = 0;
379 while ((bin = static_cast<EtaBin*>(next()))) {
380 bin->Process(*fForwardCache, *fCentralCache,
281a2bf8 381 forwardData, centralData,
c8b1a7db 382 fIsSelected, fMCCache);
bfab35d9 383 }
384
c8b1a7db 385 return true;
bfab35d9 386}
bfab35d9 387
c8b1a7db 388//_____________________________________________________________________
389Bool_t
390AliForwardMultDists::CheckEvent(const AliAODForwardMult& fwd)
391{
392 fIsSelected = AliBaseAODTask::CheckEvent(fwd);
393 // We always return true, so that we can process MC truth;
394 return true;
bfab35d9 395}
c8b1a7db 396
bfab35d9 397//____________________________________________________________________
c8b1a7db 398Bool_t
399AliForwardMultDists::Finalize()
bfab35d9 400{
67a4bb96 401 /*
bfab35d9 402 * Called at the end of the final processing of the job on the
403 * full data set (merged data)
bfab35d9 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) {
77f97e3f 432 Int_t factor = Int_t(TMath::Power(10, idx));
bfab35d9 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{
67a4bb96 461 /*
bfab35d9 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{
67a4bb96 500 /*
bfab35d9 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{
67a4bb96 539 /*
bfab35d9 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{
67a4bb96 559 /*
52b36573 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{
67a4bb96 589 /*
bfab35d9 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{
67a4bb96 625 /*
bfab35d9 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{
67a4bb96 647 /*
bfab35d9 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
77f97e3f
CHC		 680 Int_t max = Int_t(mAxis.GetXmax());
681 Int_t min = Int_t(mAxis.GetXmin());
52b36573 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{
67a4bb96 836 /*
bfab35d9 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 }
700a5a26 866
867 const TAxis * xax = hist.GetXaxis();
868 fMinBin = xax->FindFixBin(fMinEta);
869 fMaxBin = xax->FindFixBin(fMaxEta-.00001);
281a2bf8 870
52b36573 871 TString t(Form("%+5.2f<#eta<%+5.2f", fMinEta, fMaxEta));
872 fSum = CreateH1("rawDist",Form("Raw P(#it{N}_{ch}) in %s",t.Data()),fMAxis);
bfab35d9 873 fSum->SetDirectory(0);
874 fSum->GetXaxis()->SetTitle("#it{N}_{ch}");
875 fSum->GetYaxis()->SetTitle("Raw P(#it{N}_{ch})");
876 fSum->SetFillColor(kRed+1);
877 fSum->SetFillStyle(0);
878 // fSum->SetOption("hist e");
879 fSum->SetMarkerStyle(20);
880 fSum->SetMarkerColor(kRed+1);
881 fSum->SetLineColor(kBlack);
882 fSum->Sumw2();
883 l->Add(fSum);
884
52b36573 885 fCorr = CreateH2("corr",Form("C_{SPD,FMD} in %s", t.Data()),fTAxis,fTAxis);
bfab35d9 886 fCorr->SetDirectory(0);
887 fCorr->GetXaxis()->SetTitle("Forward");
888 fCorr->GetYaxis()->SetTitle("Central");
889 fCorr->SetOption("colz");
890 l->Add(fCorr);
891
892 fCoverage = new TH1D("coverage", "Fraction of bins with coverage",
893 101, -.5, 100.5);
894 fCoverage->SetDirectory(0);
895 fCoverage->SetXTitle("Fraction of bins [%]");
896 fCoverage->SetFillStyle(3001);
897 fCoverage->SetFillColor(kGreen+1);
898 l->Add(fCoverage);
899
900 if (!useMC) return;
52b36573 901 fResponse = CreateH2("response", Form("Reponse matrix in %s", t.Data()),
902 fMAxis, fTAxis);
bfab35d9 903 fResponse->SetDirectory(0);
281a2bf8 904 fResponse->SetYTitle("MC");
905 fResponse->SetXTitle("Signal");
bfab35d9 906 fResponse->SetOption("colz");
907 l->Add(fResponse);
908
52b36573 909 fTruth = CreateH1("truth",Form("True P(#it{N}_{ch}) in %s",t.Data()),fTAxis);
281a2bf8 910 fTruth->SetXTitle(fSum->GetXaxis()->GetTitle());
911 fTruth->SetYTitle(fSum->GetYaxis()->GetTitle());
bfab35d9 912 fTruth->SetLineColor(kBlack);
913 fTruth->SetFillColor(kBlue+1);
914 fTruth->SetFillStyle(0);
915 fTruth->SetDirectory(0);
916 /// fTruth->SetOption("");
917 fTruth->SetMarkerColor(kBlue+1);
918 fTruth->SetMarkerStyle(24);
281a2bf8 919 fTruth->Sumw2();
bfab35d9 920 l->Add(fTruth);
281a2bf8 921
922 fTruthAccepted = static_cast<TH1D*>(fTruth->Clone("truthAccepted"));
923 fTruthAccepted->SetTitle(Form("True (accepted) P(#it{N}_{ch}) in %s",
52b36573 924 t.Data()));
281a2bf8 925 fTruthAccepted->SetLineColor(kGray+2);
926 fTruthAccepted->SetFillColor(kOrange+2);
927 fTruthAccepted->SetDirectory(0);
928 /// fTruth->SetOption("");
929 fTruthAccepted->SetMarkerColor(kOrange+2);
930 l->Add(fTruthAccepted);
bfab35d9 931}
932//____________________________________________________________________
933void
934AliForwardMultDists::EtaBin::Process(const TH1& sumForward,
935 const TH1& sumCentral,
936 const TH2& forward,
937 const TH2& central,
281a2bf8 938 Bool_t accepted,
bfab35d9 939 const TH1* mc)
940{
67a4bb96 941 /*
bfab35d9 942 * Process a single event
943 *
944 * @param sumForward Projection of forward data
945 * @param sumCentral Projection of the central data
946 * @param forward The original forward data
947 * @param central The original central data
948 */
281a2bf8 949 if (!mc && !accepted)
950 // If we're not looking at MC data, and the event wasn't selected,
951 // we bail out - this check is already done, but we do it again
952 // for safety.
953 return;
954
bfab35d9 955 Double_t sum = 0;
956 Double_t e2sum = 0;
957 Int_t covered = 0;
958 Double_t fsum = -1;
959 Double_t csum = -1;
960 Double_t mcsum = 0;
961 Double_t mce2sum = 0;
962 for (Int_t iEta = fMinBin; iEta <= fMaxBin; iEta++) {
281a2bf8 963 if (mc) {
964 Double_t cM = mc->GetBinContent(iEta);
965 Double_t eM = mc->GetBinError(iEta);
966 mcsum += cM;
967 mce2sum += eM * eM;
968 }
969 if (!accepted)
970 // Event wasn't selected, but we still need to get the rest of
971 // the MC data.
972 continue;
973
bfab35d9 974 Double_t cF = sumForward.GetBinContent(iEta);
975 Double_t eF = sumForward.GetBinError(iEta);
976 Bool_t bF = forward.GetBinContent(iEta, 0) > 0;
977 Double_t cC = sumCentral.GetBinContent(iEta);
978 Double_t eC = sumCentral.GetBinError(iEta);
979 Bool_t bC = central.GetBinContent(iEta, 0) > 0;
bfab35d9 980 Double_t c = 0;
981 Double_t e = 0;
281a2bf8 982
bfab35d9 983 // If we have an overlap - as given by the eta-coverage,
984 // calculate the mean
985 if (bF & bC) {
986 c = (cF + cC) / 2;
987 e = TMath::Sqrt(eF*eF + eC*eC);
988 covered++;
989 }
990 // Else, if we have coverage from forward, use that
991 else if (bF) {
992 c = cF;
993 e = eF;
994 covered++;
995 }
996 // Else, if we have covrage from central, use that
997 else if (bC) {
998 c = cC;
999 e = eC;
1000 covered++;
1001 }
1002 // Else, we have incomplete coverage
1003
1004 if (bF) {
1005 if (fsum < 0) fsum = 0;
1006 fsum += cF;
1007 }
1008 if (bC) {
1009 if (csum < 0) csum = 0;
1010 csum += cC;
1011 }
1012
1013 sum += c;
1014 e2sum += e*e;
1015 }
1016
281a2bf8 1017 if (accepted) {
1018 // Only update the histograms if the event was triggered.
1019 // Fill with weight
1020 Double_t w = 1; // e2sum > 0 ? 1/TMath::Sqrt(e2sum) : 1
1021 fSum->Fill(sum, w);
52b36573 1022 fCorr->Fill((fsum<0?0:fsum), (csum<0?0:csum));
281a2bf8 1023
1024 Int_t nTotal = fMaxBin - fMinBin + 1;
1025 fCoverage->Fill(100*float(covered) / nTotal);
1026 }
bfab35d9 1027
1028 if (mc) {
281a2bf8 1029 Double_t w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
bfab35d9 1030 if (fTruth) {
bfab35d9 1031 fTruth->Fill(mcsum, w);
1032 }
281a2bf8 1033 if (accepted) {
1034 if (fResponse)
1035 // Only fill response matrix for accepted events
1036 fResponse->Fill(sum, mcsum);
1037 if (fTruthAccepted)
1038 fTruthAccepted->Fill(mcsum, w);
1039 }
bfab35d9 1040 }
1041}
1042//____________________________________________________________________
1043void
52b36573 1044AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out)
bfab35d9 1045{
67a4bb96 1046 /*
bfab35d9 1047 * Called at the end of the final processing of the job on the
1048 * full data set (merged data)
1049 *
1050 * @param in Input list
1051 * @param out Output list
1052 * @param maxN Maximum number of @f$N_{ch}@f$ to consider
1053 */
1054 TList* ip = FindParent(in, false);
1055 if (!ip) {
1056 AliErrorF("Parent folder %s not found in input", ParentName());
1057 return;
1058 }
1059
1060 TList* i = dynamic_cast<TList*>(ip->FindObject(GetName()));
1061 if (!i) {
1062 AliErrorF("List %s not found in input", GetName());
1063 return;
1064 }
1065
1066 TList* op = FindParent(out, true);
1067 TList* o = static_cast<TList*>(i->Clone());
1068 o->SetOwner();
1069 op->Add(o);
1070
281a2bf8 1071 fSum = static_cast<TH1*>(o->FindObject("rawDist"));
1072 fTruth = static_cast<TH1*>(o->FindObject("truth"));
1073 fTruthAccepted = static_cast<TH1*>(o->FindObject("truthAccepted"));
1074
1075 TH1* hists[] = { fSum, fTruth, fTruthAccepted, 0 };
bfab35d9 1076 TH1** phist = hists;
1077 while (*phist) {
1078 TH1* h = *phist;
1079 if (h) {
52b36573 1080 Int_t maxN = h->GetNbinsX();
bfab35d9 1081 Double_t intg = h->Integral(1, maxN);
281a2bf8 1082 h->Scale(1. / intg, "width");
bfab35d9 1083 }
1084 phist++;
1085 }
281a2bf8 1086
1087 if (fTruth && fTruthAccepted) {
1088 TH1* trgVtx = static_cast<TH1*>(fTruthAccepted->Clone("triggerVertex"));
1089 TString tit(fTruth->GetTitle());
1090 tit.ReplaceAll("True P(#it{N}_{ch})", "C_{trigger,vertex}");
1091 trgVtx->SetTitle(tit);
1092 trgVtx->SetYTitle("P_{MC}(#it{N}_{ch})/P_{MC,accepted}(#it{N}_{ch})");
1093 trgVtx->Divide(fTruth);
1094 trgVtx->SetDirectory(0);
1095 o->Add(trgVtx);
1096 }
1097
bfab35d9 1098}
1099//====================================================================
1100//
1101// EOF
1102//