Compilation with Root6: TH1::GetXaxis returns now const TAxis*
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardMultDists.cxx
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 {
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 #if 0
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   }
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   /** 
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   }
110 #endif
111 }
112
113 //====================================================================
114 AliForwardMultDists::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 //____________________________________________________________________
121 void
122 AliForwardMultDists::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 //____________________________________________________________________
131 const TAxis&
132 AliForwardMultDists::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 { 
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 //====================================================================
162 AliForwardMultDists::AliForwardMultDists()
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 //____________________________________________________________________
178 AliForwardMultDists::AliForwardMultDists(const char* name)
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
202 AliForwardMultDists::Book()
203 {
204   /* 
205    * Create output objects - called at start of job in slave 
206    * 
207    */
208   return true;
209 }
210 //____________________________________________________________________
211 Bool_t
212 AliForwardMultDists::PreData()
213 {
214   /* 
215    * Set-up internal structures on first seen event 
216    * 
217    */
218   AliAODEvent*       aod     = AliForwardUtil::GetAODEvent(this);
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", fMinIpZ));
224   fSums->Add(AliForwardUtil::MakeParameter("maxIpZ", fMaxIpZ));
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
299 AliForwardMultDists::Event(AliAODEvent& aod)
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 <= 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
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;
331     else 
332       mcTrg = true;
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   }
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) {
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   }
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()))) {
380     bin->Process(*fForwardCache, *fCentralCache,
381                  forwardData,    centralData, 
382                  fIsSelected,    fMCCache);
383   }
384
385   return true;
386 }
387
388 //_____________________________________________________________________
389 Bool_t 
390 AliForwardMultDists::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;
395 }
396
397 //____________________________________________________________________
398 Bool_t
399 AliForwardMultDists::Finalize()
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
529 AliForwardMultDists::AddBin(const BinSpec& spec)
530 {
531   AddBin(spec.fEtaMin, spec.fEtaMax, spec.Axis());
532 }
533
534 //____________________________________________________________________
535 void
536 AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax, 
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
556 AliForwardMultDists::AddBin(Double_t etaLow, Double_t etaMax, 
557                             UShort_t nMax, UShort_t nDiv)
558 {
559   /* 
560    * Add an @f$\eta@f$ bin
561    * 
562    * @param etaLow Low cut on @f$\eta@f$
563    * @param etaMax High cut on @f$\eta@f$
564    */
565   // Symmetric bin
566   if (etaMax >= kInvalidEta) { 
567     etaLow = -TMath::Abs(etaLow);
568     etaMax = +TMath::Abs(etaLow);
569   }
570   TAxis mAxis((nMax+1)/nDiv, -.5, nMax+.5);
571   EtaBin* bin = new EtaBin(etaLow, etaMax, mAxis);
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
587 AliForwardMultDists::Print(Option_t* option) const 
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 //====================================================================
610 AliForwardMultDists::EtaBin::EtaBin() 
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 //____________________________________________________________________
631 AliForwardMultDists::EtaBin::EtaBin(Double_t minEta, Double_t maxEta,
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 //____________________________________________________________________
686 AliForwardMultDists::EtaBin::EtaBin(const EtaBin& o) 
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),
699     fTruthAccepted(o.fTruthAccepted),
700     fCoverage(o.fCoverage)
701 {}
702 //____________________________________________________________________
703 AliForwardMultDists::EtaBin&
704 AliForwardMultDists::EtaBin::operator=(const EtaBin& o) 
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;
719   fTruthAccepted = o.fTruthAccepted;
720   fCoverage      = o.fCoverage;
721
722   return *this;
723 }
724 //____________________________________________________________________
725 Bool_t
726 AliForwardMultDists::EtaBin::IsSymmetric() const
727 {
728   return TMath::Abs(fMaxEta + fMinEta) < 1e-6;
729 }
730 //____________________________________________________________________
731 Bool_t
732 AliForwardMultDists::EtaBin::IsNegative() const
733 {
734   return TMath::Abs(fMaxEta) < 1e-6 && fMinEta < 0;
735 }
736 //____________________________________________________________________
737 Bool_t
738 AliForwardMultDists::EtaBin::IsPositive() const
739 {
740   return TMath::Abs(fMinEta) < 1e-6 && fMaxEta > 0;
741 }
742 //____________________________________________________________________
743 const char*
744 AliForwardMultDists::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 //____________________________________________________________________
752 TList*
753 AliForwardMultDists::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
781 //____________________________________________________________________
782 TH1*
783 AliForwardMultDists::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 //____________________________________________________________________
796 TH2*
797 AliForwardMultDists::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
831 //____________________________________________________________________
832 void
833 AliForwardMultDists::EtaBin::SetupForData(TList* list, const TH2& hist, 
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);
916   /// fTruth->SetOption("");
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);
928   /// fTruth->SetOption("");
929   fTruthAccepted->SetMarkerColor(kOrange+2);
930   l->Add(fTruthAccepted);
931 }
932 //____________________________________________________________________
933 void
934 AliForwardMultDists::EtaBin::Process(const TH1& sumForward, 
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
1044 AliForwardMultDists::EtaBin::Terminate(TList* in, TList* out)
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 //