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