AliPhysics  5403132 (5403132)
AliForwardMultDists.cxx
Go to the documentation of this file.
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>
12 #include <TError.h>
13 #include <iostream>
14 #include <iomanip>
15 
16 //====================================================================
17 namespace {
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  };
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 #if 0
55 
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  }
80 #endif
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  }
96 #if 0
97 
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  }
110 #endif
111 }
112 
113 //====================================================================
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 //____________________________________________________________________
121 void
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 //____________________________________________________________________
131 const TAxis&
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 {
141  Int_t n = Int_t(fN.GetSum());
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 }
160 
161 //====================================================================
163  : AliBaseAODTask(),
164  fBins(),
165  fSymmetric(0),
166  fNegative(0),
167  fPositive(0),
168  fMCVertex(0),
169  fDiag(0),
170  fForwardCache(0),
171  fCentralCache(0),
172  fMCCache(0),
173  fUsePhiAcc(true),
174  fIsSelected(false)
175 {}
176 
177 //____________________________________________________________________
179  : AliBaseAODTask(name,"AliForwardMultDists"),
180  fBins(),
181  fSymmetric(0),
182  fNegative(0),
183  fPositive(0),
184  fMCVertex(0),
185  fDiag(0),
186  fForwardCache(0),
187  fCentralCache(0),
188  fMCCache(0),
189  fUsePhiAcc(true),
190  fIsSelected(false)
191 {
192  /*
193  * User constructor
194  *
195  * @param name Name of the task
196  */
197 }
198 
199 
200 //____________________________________________________________________
201 Bool_t
203 {
204  /*
205  * Create output objects - called at start of job in slave
206  *
207  */
208  return true;
209 }
210 //____________________________________________________________________
211 Bool_t
213 {
214  /*
215  * Set-up internal structures on first seen event
216  *
217  */
219  AliAODForwardMult* forward = GetForward(*aod);
220  const TH2D& hist = forward->GetHistogram();
221  Bool_t useMC = GetPrimary(*aod) != 0;
222 
223  fSums->Add(AliForwardUtil::MakeParameter("minIpZ", fIPzAxis.GetXmin()));
224  fSums->Add(AliForwardUtil::MakeParameter("maxIpZ", fIPzAxis.GetXmax()));
225 
226  const TAxis* xaxis = hist.GetXaxis();
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();
256 
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();
262 
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  }
269 
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;
296 }
297 //____________________________________________________________________
298 Bool_t
300 {
301  /*
302  * Analyse a single event
303  *
304  * @param aod Input event
305  */
306  AliAODForwardMult* forward = GetForward(aod, false, true);
307  if (!forward) return false;
308 
309  AliAODCentralMult* central = GetCentral(aod, false, true);
310  if (!central) return false;
311 
312  TH2* primary = GetPrimary(aod);
313 
314  const TH2& forwardData = forward->GetHistogram();
315  const TH2& centralData = central->GetHistogram();
316 
317  Double_t vz = forward->GetIpZ();
318  Bool_t trg = forward->IsTriggerBits(fTriggerMask);
319  Bool_t vtx = (vz <= fIPzAxis.GetXmax() && vz >= fIPzAxis.GetXmin());
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
324  if (primary) {
325  // For MC, we need to check if we should process the event for
326  // MC information or not.
329  // Bail out if this event is not MC NSD event
330  ok = false;
331  else
332  mcTrg = true;
333  Double_t mcVz = primary->GetBinContent(0,0);
334  fMCVertex->Fill(mcVz);
335  if (mcVz > fIPzAxis.GetXmax() || mcVz < fIPzAxis.GetXmin())
336  // Bail out if this event was not in the valid range
337  ok = false;
338  else
339  mcVtx = true;
340  }
341  // Fill diagnostics
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) {
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  }
365  if (!fIsSelected && !ok) {
366  // Printf("Event is neither accepted by analysis nor by MC");
367  return false;
368  }
369 
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()))) {
381  forwardData, centralData,
383  }
384 
385  return true;
386 }
387 
388 //_____________________________________________________________________
389 Bool_t
391 {
393  // We always return true, so that we can process MC truth;
394  return true;
395 }
396 
397 //____________________________________________________________________
398 Bool_t
400 {
401  /*
402  * Called at the end of the final processing of the job on the
403  * full data set (merged data)
404  */
405  fResults->Add(fSums->FindObject("triggers")->Clone());
406  fResults->Add(fSums->FindObject("status")->Clone());
407  fResults->Add(fSums->FindObject("diagnostics")->Clone());
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()))) {
417  bin->Terminate(fSums, fResults);
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 
424  TH1* hh[] = { bin->fSum, bin->fTruth, bin->fTruthAccepted, 0 };
425  TH1** ph = hh;
426 
427  Int_t idx = sta->GetHists() ? sta->GetHists()->GetEntries() : 0;
428  if (bin->fTruth) idx /= 2;
429 
430  Int_t mkrBits = GetIndexMarker(idx);
431  while (*ph) {
432  Int_t factor = Int_t(TMath::Power(10, idx));
433  TH1* h = static_cast<TH1*>((*ph)->Clone());
434  h->SetDirectory(0);
435  h->Scale(factor);
436  h->SetTitle(Form("%s (#times%d)", h->GetTitle(), Int_t(factor)));
437  h->SetMarkerStyle(GetMarkerStyle(mkrBits));
438  mkrBits ^= kHollow;
439 
440  sta->Add(h, "p");
441  ph++;
442  }
443  }
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"));
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 
454  return true;
455 }
456 
457 //____________________________________________________________________
458 void
459 AliForwardMultDists::ProjectX(const TH2& input, TH1& cache, Bool_t usePhiAcc)
460 {
461  /*
462  * Project a 2D histogram into a 1D histogram taking care to use
463  * either the @f$\phi@f$ acceptance stored in the overflow bins, or
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
468  * @param usePhiAcc If true, use the @f$\phi@f$ acceptance stored in
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 //____________________________________________________________________
497 void
498 AliForwardMultDists::ProjectX(const TH2* input, TH1* cache)
499 {
500  /*
501  * Project on @f$\eta@f$ axis. If any of the pointers passed is
502  * zero, do nothing.
503  *
504  * @param input
505  * @param cache
506  */
507  if (!input || !cache) return;
508  cache->Reset();
509 
510  Int_t nPhi = input->GetNbinsY();
511  Int_t nEta = input->GetNbinsX();
512 
513  for (Int_t iEta = 1; iEta <= nEta; iEta++) {
514  Double_t sum = 0;
515  Double_t e2sum = 0;
516  for (Int_t iPhi = 1; iPhi <= nPhi; iPhi++) {
517  Double_t c = input->GetBinContent(iEta, iPhi);
518  Double_t e = input->GetBinError(iEta, iPhi);
519 
520  sum += c;
521  e2sum += e * e;
522  }
523  cache->SetBinContent(iEta, sum);
524  cache->SetBinError(iEta, TMath::Sqrt(e2sum));
525  }
526 }
527 //____________________________________________________________________
528 void
530 {
531  AddBin(spec.fEtaMin, spec.fEtaMax, spec.Axis());
532 }
533 
534 //____________________________________________________________________
535 void
537  const TAxis& mAxis)
538 {
539  /*
540  * Add an @f$\eta@f$ bin
541  *
542  * @param etaLow Low cut on @f$\eta@f$
543  * @param etaMax High cut on @f$\eta@f$
544  */
545  // Symmetric bin
546  if (etaMax >= kInvalidEta) {
547  etaLow = -TMath::Abs(etaLow);
548  etaMax = +TMath::Abs(etaLow);
549  }
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 //____________________________________________________________________
555 void
557  UShort_t nMax, UShort_t nDiv)
558 {
559  /*
560  * Add an @f$\eta@f$ bin
561  *
562  * @param etaLow Low cut on @f$\eta@f$
563  * @param etaMax High cut on @f$\eta@f$
564  */
565  // Symmetric bin
566  if (etaMax >= kInvalidEta) {
567  etaLow = -TMath::Abs(etaLow);
568  etaMax = +TMath::Abs(etaLow);
569  }
570  TAxis mAxis((nMax+1)/nDiv, -.5, nMax+.5);
571  EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
572  // AliInfoF("Adding bin %f, %f: %s", etaLow, etaMax, bin->GetName());
573  fBins.Add(bin);
574 }
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 
585 //____________________________________________________________________
586 void
588 {
589  /*
590  * Print this task
591  *
592  * @param option Not used
593  */
594  AliBaseAODTask::Print(option);
595  gROOT->IncreaseDirLevel();
596  PFB("Use phi acceptance", fUsePhiAcc);
597  // AliForwardUtil::PrintName("Bins");
598 
599  Int_t first = true;
600  TIter next(&fBins);
601  EtaBin* bin = 0;
602  while ((bin = static_cast<EtaBin*>(next()))) {
603  PFV(first ? "Bins" : "", bin->GetName());
604  first = false;
605  }
606  gROOT->DecreaseDirLevel();
607 }
608 
609 //====================================================================
611  : fName(""),
612  fMAxis(1,0,0),
613  fTAxis(1,0,0),
614  fMinEta(0),
615  fMaxEta(0),
616  fMinBin(0),
617  fMaxBin(0),
618  fSum(0),
619  fCorr(0),
620  fResponse(0),
621  fTruth(0),
622  fTruthAccepted(0),
623  fCoverage(0)
624 {
625  /*
626  * I/O constructor
627  */
628 }
629 
630 //____________________________________________________________________
632  const TAxis& mAxis)
633  : fName(""),
634  fMAxis(1,0,0),
635  fTAxis(1,0,0),
636  fMinEta(minEta),
637  fMaxEta(maxEta),
638  fMinBin(0),
639  fMaxBin(0),
640  fSum(0),
641  fCorr(0),
642  fResponse(0),
643  fTruth(0),
644  fTruthAccepted(0),
645  fCoverage(0)
646 {
647  /*
648  * User constructor
649  *
650  * @param minEta Least @f$\eta@f$ to consider
651  * @param maxEta Largest @f$\eta@f$ to consider
652  */
653  fName = TString::Format("%+05.2f_%+05.2f", fMinEta, fMaxEta);
654  fName.ReplaceAll("-", "m");
655  fName.ReplaceAll("+", "p");
656  fName.ReplaceAll(".", "d");
657 
658  // Copy to other our object
659  mAxis.Copy(fMAxis);
660 
661  if (mAxis.GetXbins() && mAxis.GetXbins()->GetArray()) {
662  const TArrayD& mA = *(mAxis.GetXbins());
663  TArrayD tA(mA.GetSize());
664  Int_t j = 0;
665  Double_t min = mA[0];
666  tA[0] = min;
667  for (Int_t i = 1; i < mA.GetSize(); i++) {
668  Double_t d = mA[i] - min;
669  if (d < 1)
670  // Not full integer bin
671  continue;
672  tA[j+1] = tA[j] + Int_t(d);
673  min = tA[j+1];
674  j++;
675  }
676  fTAxis.Set(j, tA.GetArray());
677  }
678  else {
679  // Rounded down maximum and minimum
680  Int_t max = Int_t(mAxis.GetXmax());
681  Int_t min = Int_t(mAxis.GetXmin());
682  fTAxis.Set((max-min)+1, min-.5, max+.5);
683  }
684 }
685 //____________________________________________________________________
687  : TObject(o),
688  fName(o.fName),
689  fMAxis(o.fMAxis),
690  fTAxis(o.fTAxis),
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),
701 {}
702 //____________________________________________________________________
705 {
706  if (&o == this) return *this;
707 
708  fName = o.fName;
709  fMAxis = o.fMAxis;
710  fTAxis = o.fTAxis;
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;
720  fCoverage = o.fCoverage;
721 
722  return *this;
723 }
724 //____________________________________________________________________
725 Bool_t
727 {
728  return TMath::Abs(fMaxEta + fMinEta) < 1e-6;
729 }
730 //____________________________________________________________________
731 Bool_t
733 {
734  return TMath::Abs(fMaxEta) < 1e-6 && fMinEta < 0;
735 }
736 //____________________________________________________________________
737 Bool_t
739 {
740  return TMath::Abs(fMinEta) < 1e-6 && fMaxEta > 0;
741 }
742 //____________________________________________________________________
743 const char*
745 {
746  if (IsSymmetric()) return "symmetric";
747  else if (IsNegative()) return "negative";
748  else if (IsPositive()) return "positive";
749  return "other";
750 }
751 //____________________________________________________________________
752 TList*
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 
781 //____________________________________________________________________
782 TH1*
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 //____________________________________________________________________
796 TH2*
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 
831 //____________________________________________________________________
832 void
834  Bool_t useMC)
835 {
836  /*
837  * Set-up internal structures on first event.
838  *
839  * @param list List to add information to
840  * @param hist Template histogram
841  * @param max Maximum number of particles
842  * @param useMC Whether to set-up for MC input
843  */
844  TList* p = FindParent(list, true);
845  TList* l = new TList;
846  l->SetName(GetName());
847  l->SetOwner();
848  p->Add(l);
849  TList* le = static_cast<TList*>(p->FindObject("lowEdges"));
850  TList* he = static_cast<TList*>(p->FindObject("highEdges"));
851  if (!le || !he) {
852  AliError("Failed to get bin edge lists from parent");
853  return;
854  }
855  else {
856  Int_t n = le->GetEntries();
857  TParameter<double>* lp =
858  new TParameter<double>(Form("minEta%02d", n), fMinEta);
859  TParameter<double>* hp =
860  new TParameter<double>(Form("maxEta%02d", n), fMaxEta);
861  lp->SetMergeMode('f');
862  hp->SetMergeMode('f');
863  le->Add(lp);
864  he->Add(hp);
865  }
866 
867  const TAxis * xax = hist.GetXaxis();
868  fMinBin = xax->FindFixBin(fMinEta);
869  fMaxBin = xax->FindFixBin(fMaxEta-.00001);
870 
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);
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 
885  fCorr = CreateH2("corr",Form("C_{SPD,FMD} in %s", t.Data()),fTAxis,fTAxis);
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;
901  fResponse = CreateH2("response", Form("Reponse matrix in %s", t.Data()),
902  fMAxis, fTAxis);
903  fResponse->SetDirectory(0);
904  fResponse->SetYTitle("MC");
905  fResponse->SetXTitle("Signal");
906  fResponse->SetOption("colz");
907  l->Add(fResponse);
908 
909  fTruth = CreateH1("truth",Form("True P(#it{N}_{ch}) in %s",t.Data()),fTAxis);
910  fTruth->SetXTitle(fSum->GetXaxis()->GetTitle());
911  fTruth->SetYTitle(fSum->GetYaxis()->GetTitle());
912  fTruth->SetLineColor(kBlack);
913  fTruth->SetFillColor(kBlue+1);
914  fTruth->SetFillStyle(0);
915  fTruth->SetDirectory(0);
917  fTruth->SetMarkerColor(kBlue+1);
918  fTruth->SetMarkerStyle(24);
919  fTruth->Sumw2();
920  l->Add(fTruth);
921 
922  fTruthAccepted = static_cast<TH1D*>(fTruth->Clone("truthAccepted"));
923  fTruthAccepted->SetTitle(Form("True (accepted) P(#it{N}_{ch}) in %s",
924  t.Data()));
925  fTruthAccepted->SetLineColor(kGray+2);
926  fTruthAccepted->SetFillColor(kOrange+2);
927  fTruthAccepted->SetDirectory(0);
929  fTruthAccepted->SetMarkerColor(kOrange+2);
930  l->Add(fTruthAccepted);
931 }
932 //____________________________________________________________________
933 void
935  const TH1& sumCentral,
936  const TH2& forward,
937  const TH2& central,
938  Bool_t accepted,
939  const TH1* mc)
940 {
941  /*
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  */
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 
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++) {
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 
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;
980  Double_t c = 0;
981  Double_t e = 0;
982 
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 
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);
1022  fCorr->Fill((fsum<0?0:fsum), (csum<0?0:csum));
1023 
1024  Int_t nTotal = fMaxBin - fMinBin + 1;
1025  fCoverage->Fill(100*float(covered) / nTotal);
1026  }
1027 
1028  if (mc) {
1029  Double_t w = 1; // mce2sum > 0 ? 1/TMath::Sqrt(mce2sum) : 1
1030  if (fTruth) {
1031  fTruth->Fill(mcsum, w);
1032  }
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  }
1040  }
1041 }
1042 //____________________________________________________________________
1043 void
1045 {
1046  /*
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 
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 };
1076  TH1** phist = hists;
1077  while (*phist) {
1078  TH1* h = *phist;
1079  if (h) {
1080  Int_t maxN = h->GetNbinsX();
1081  Double_t intg = h->Integral(1, maxN);
1082  h->Scale(1. / intg, "width");
1083  }
1084  phist++;
1085  }
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 
1098 }
1099 //====================================================================
1100 //
1101 // EOF
1102 //
static TH2 * CreateH2(const char *name, const char *title, const TAxis &xAxis, const TAxis &yAxis)
virtual Bool_t CheckEvent(const AliAODForwardMult &forward)
const char * ParentName() const
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
TCanvas * c
Definition: TestFitELoss.C:172
void AddBin(const BinSpec &spec)
static AliAODEvent * GetAODEvent(AliAnalysisTaskSE *task)
void SetupForData(TList *list, const TH2 &hist, Bool_t useMC)
Float_t GetIpZ() const
TList * FindParent(TList *l, Bool_t create=true) const
#define PFV(N, VALUE)
static void ProjectX(const TH2 &input, TH1 &cache, Bool_t usePhiAcc=true)
Per-event per bin.
int Int_t
Definition: External.C:63
static TH1 * CreateH1(const char *name, const char *title, const TAxis &xAxis)
static Bool_t IsTriggerBits(UInt_t bits, UInt_t trg)
Various utilities used in PWGLF/FORWARD.
void Push(UShort_t n, Double_t d)
Definition: External.C:228
Definition: External.C:212
BinSpec(Double_t etaMin, Double_t etaMax, Double_t nchLow)
void Terminate(TList *in, TList *out)
void Print(Option_t *option="") const
void Print(Option_t *option="") const
EtaBin & operator=(const EtaBin &o)
static TObject * MakeParameter(const char *name, UShort_t value)
const TH2D & GetHistogram() const
const char * fwd
Class to make raw distributions.
const char * GetName() const
Definition: External.C:220
AliAODForwardMult * GetForward(const AliAODEvent &aod, Bool_t mc=false, Bool_t verb=true)
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
TH2D * GetPrimary(const AliAODEvent &aod)
AliAODCentralMult * GetCentral(const AliAODEvent &aod, Bool_t mc=false, Bool_t verb=true)
bool Bool_t
Definition: External.C:53
Bool_t CheckEvent(const AliAODForwardMult &fwd)
#define PFB(N, FLAG)
Definition: External.C:196
Bool_t Event(AliAODEvent &aod)
const TH2D & GetHistogram() const
void Process(const TH1 &sumForward, const TH1 &sumCentral, const TH2 &forward, const TH2 &central, Bool_t accepted, const TH1 *mc)