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