Various fixes, and possibility to set trigger efficiency
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliBasedNdetaTask.cxx
1 //====================================================================
2 #include "AliBasedNdetaTask.h"
3 #include <TMath.h>
4 #include <TH2D.h>
5 #include <TH1D.h>
6 #include <THStack.h>
7 #include <TList.h>
8 #include <AliAnalysisManager.h>
9 #include <AliAODEvent.h>
10 #include <AliAODHandler.h>
11 #include <AliAODInputHandler.h>
12 #include "AliForwardUtil.h"
13 #include "AliAODForwardMult.h"
14 #include <TFile.h>
15 #include <TStyle.h>
16
17 //____________________________________________________________________
18 AliBasedNdetaTask::AliBasedNdetaTask()
19   : AliAnalysisTaskSE(), 
20     fSums(0),           // Container of sums 
21     fOutput(0),         // Container of output 
22     fVtxMin(0),         // Minimum v_z
23     fVtxMax(0),         // Maximum v_z
24     fTriggerMask(0),    // Trigger mask 
25     fRebin(0),          // Rebinning factor 
26     fCutEdges(false), 
27     fSymmetrice(true),
28     fCorrEmpty(true), 
29     fUseROOTProj(false),
30     fTriggerEff(1),
31   fTriggerEff0(1),
32     fShapeCorr(0),
33     fListOfCentralities(0),
34     fSNNString(0),
35     fSysString(0),
36     fCent(0),
37     fCentAxis(0),
38     fNormalizationScheme(kFull), 
39     fSchemeString(0), 
40     fTriggerString(0),
41     fFinalMCCorrFile("")
42 {
43   // 
44   // Constructor
45   // 
46   DGUARD(0,0,"Default construction of AliBasedNdetaTask");
47 }
48
49 //____________________________________________________________________
50 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
51   : AliAnalysisTaskSE(name), 
52     fSums(0),           // Container of sums 
53     fOutput(0),         // Container of output 
54     fVtxMin(-10),       // Minimum v_z
55     fVtxMax(10),        // Maximum v_z
56     fTriggerMask(AliAODForwardMult::kInel), 
57     fRebin(5),          // Rebinning factor 
58     fCutEdges(false), 
59     fSymmetrice(true),
60     fCorrEmpty(true), 
61     fUseROOTProj(false),
62     fTriggerEff(1),
63     fTriggerEff0(1),
64     fShapeCorr(0),
65     fListOfCentralities(0),
66     fSNNString(0),
67     fSysString(0),
68     fCent(0),
69     fCentAxis(0),
70     fNormalizationScheme(kFull), 
71     fSchemeString(0),
72     fTriggerString(0),
73     fFinalMCCorrFile("")
74 {
75   // 
76   // Constructor
77   // 
78   DGUARD(0,0,"Named construction of AliBasedNdetaTask: %s", name);
79   fListOfCentralities = new TObjArray(1);
80   
81   // Set the normalisation scheme 
82   SetNormalizationScheme(kFull);
83
84   // Set the trigger mask
85   SetTriggerMask(AliAODForwardMult::kInel);
86
87   // Output slot #1 writes into a TH1 container
88   DefineOutput(1, TList::Class()); 
89   DefineOutput(2, TList::Class()); 
90 }
91
92 //____________________________________________________________________
93 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
94   : AliAnalysisTaskSE(o), 
95     fSums(o.fSums),             // TList* - Container of sums 
96     fOutput(o.fOutput),         // Container of output 
97     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
98     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
99     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask 
100     fRebin(o.fRebin),           // Int_t - Rebinning factor 
101     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
102     fSymmetrice(o.fSymmetrice),
103     fCorrEmpty(o.fCorrEmpty), 
104     fUseROOTProj(o.fUseROOTProj),
105     fTriggerEff(o.fTriggerEff),
106     fTriggerEff0(o.fTriggerEff0),
107     fShapeCorr(o.fShapeCorr),
108     fListOfCentralities(o.fListOfCentralities),
109     fSNNString(o.fSNNString),
110     fSysString(o.fSysString),
111     fCent(o.fCent),
112     fCentAxis(o.fCentAxis),
113     fNormalizationScheme(o.fNormalizationScheme), 
114     fSchemeString(o.fSchemeString), 
115     fTriggerString(o.fTriggerString),
116     fFinalMCCorrFile(o.fFinalMCCorrFile)
117 {
118   DGUARD(0,0,"Copy construction of AliBasedNdetaTask");
119 }
120
121 //____________________________________________________________________
122 AliBasedNdetaTask::~AliBasedNdetaTask()
123 {
124   // 
125   // Destructor
126   // 
127   DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask");
128   if (fSums) { 
129     fSums->Delete();
130     delete fSums;
131     fSums = 0;
132   }
133   if (fOutput) { 
134     fOutput->Delete();
135     delete fOutput;
136     fOutput = 0;
137   }
138 }
139
140 //________________________________________________________________________
141 void 
142 AliBasedNdetaTask::SetDebugLevel(Int_t lvl)
143 {
144   AliAnalysisTaskSE::SetDebugLevel(lvl);
145   for (Int_t i = 0; i < fListOfCentralities->GetEntries(); i++) { 
146     CentralityBin* bin = 
147       static_cast<CentralityBin*>(fListOfCentralities->At(i));
148     bin->SetDebugLevel(lvl);
149   }
150 }
151
152 //________________________________________________________________________
153 void 
154 AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
155 {
156   DGUARD(fDebug,3,"Set centrality axis, %d bins", n);
157   if (!fCentAxis) { 
158     fCentAxis = new TAxis();
159     fCentAxis->SetName("centAxis");
160     fCentAxis->SetTitle("Centrality [%]");
161   }
162   TArrayD dbins(n+1);
163   for (UShort_t i = 0; i <= n; i++) 
164     dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
165   fCentAxis->Set(n, dbins.GetArray());
166 }
167     
168 //________________________________________________________________________
169 void 
170 AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
171 {
172   // 
173   // Add a centrality bin 
174   // 
175   // Parameters:
176   //    low  Low cut
177   //    high High cut
178   //
179   DGUARD(fDebug,3,"Add a centrality bin [%f,%f] @ %d", low, high, at);
180   CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
181   bin->SetDebugLevel(fDebug);
182   fListOfCentralities->AddAtAndExpand(bin, at);
183 }
184
185 //________________________________________________________________________
186 AliBasedNdetaTask::CentralityBin*
187 AliBasedNdetaTask::MakeCentralityBin(const char* name, 
188                                      Short_t low, Short_t high) const
189 {
190   // 
191   // Make a centrality bin 
192   // 
193   // Parameters:
194   //    name  Name used for histograms
195   //    low   Low cut in percent
196   //    high  High cut in percent
197   // 
198   // Return:
199   //    A newly created centrality bin 
200   //
201   DGUARD(fDebug,3,"Make a centrality bin [%f,%f]: %s", low, high, name);
202   return new CentralityBin(name, low, high);
203 }
204
205 #define TESTAPPEND(SCHEME,BIT,STRING) \
206   do { if (!(SCHEME & BIT)) break;                                      \
207     if (!s.IsNull()) s.Append(","); s.Append(STRING); } while(false) 
208   
209 //________________________________________________________________________
210 const Char_t*
211 AliBasedNdetaTask::NormalizationSchemeString(UShort_t scheme)
212 {
213   // Create a string from normalization scheme bits 
214   static TString s;
215   s = "";
216
217   if (scheme == kNone) 
218     return s.Data();
219   if (scheme == kFull) { 
220     s = "FULL";
221     return s.Data();
222   }
223   TESTAPPEND(scheme, kEventLevel,        "EVENT");
224   TESTAPPEND(scheme, kShape,             "SHAPE");
225   TESTAPPEND(scheme, kBackground,        "BACKGROUND");
226   TESTAPPEND(scheme, kTriggerEfficiency, "TRIGGER");
227   TESTAPPEND(scheme, kZeroBin,           "ZEROBIN");
228
229   return s.Data();
230 }
231 //________________________________________________________________________
232 UShort_t
233 AliBasedNdetaTask::ParseNormalizationScheme(const char* what)
234 {
235   UShort_t    scheme = 0;
236   TString     twhat(what);
237   twhat.ToUpper();
238   TObjString* opt;
239   TIter       next(twhat.Tokenize(" ,|"));
240   while ((opt = static_cast<TObjString*>(next()))) { 
241     TString s(opt->GetString());
242     if      (s.IsNull()) continue;
243     Bool_t add = true;
244     switch (s[0]) { 
245     case '-': add = false; // Fall through 
246     case '+': s.Remove(0,1);  // Remove character 
247     }
248     UShort_t bit = 0;
249     if      (s.CompareTo("EVENT")     == 0) bit = kEventLevel;
250     else if (s.CompareTo("SHAPE")     == 0) bit = kShape;
251     else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
252     else if (s.CompareTo("TRIGGER")   == 0) bit = kTriggerEfficiency;
253     else if (s.CompareTo("FULL")      == 0) bit = kFull;
254     else if (s.CompareTo("NONE")      == 0) bit = kNone;
255     else if (s.CompareTo("ZEROBIN")   == 0) bit = kZeroBin;
256     else 
257       ::Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
258     if (add) scheme |= bit;
259     else     scheme ^= bit;
260   }
261   return scheme;
262 }  
263 //________________________________________________________________________
264 void 
265 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
266 {
267   // 
268   // Set normalisation scheme 
269   // 
270   DGUARD(fDebug,3,"Set the normalization scheme: %s", what);
271   SetNormalizationScheme(ParseNormalizationScheme(what));
272 }
273 //________________________________________________________________________
274 void 
275 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme) 
276 {
277   DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme);
278   fNormalizationScheme = scheme; 
279   if (fSchemeString) delete fSchemeString;
280   fSchemeString = AliForwardUtil::MakeParameter("scheme", scheme);
281 }
282 //________________________________________________________________________
283 void 
284 AliBasedNdetaTask::SetTriggerMask(const char* mask)
285 {
286   // 
287   // Set the trigger maskl 
288   // 
289   // Parameters:
290   //    mask Trigger mask
291   //
292   DGUARD(fDebug,3,"Set the trigger mask: %s", mask);
293   SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
294 }
295 //________________________________________________________________________
296 void 
297 AliBasedNdetaTask::SetTriggerMask(UShort_t mask) 
298
299   DGUARD(fDebug,3,"Set the trigger mask: 0x%0x", mask);
300   fTriggerMask = mask; 
301   if (fTriggerString) delete fTriggerString;
302   fTriggerString = AliForwardUtil::MakeParameter("trigger", fTriggerMask);
303 }
304
305 //________________________________________________________________________
306 void 
307 AliBasedNdetaTask::SetShapeCorrection(const TH2F* c)
308 {
309   // 
310   // Set the shape correction (a.k.a., track correction) for selected
311   // trigger(s)
312   // 
313   // Parameters:
314   //    h Correction
315   //
316   DGUARD(fDebug,3,"Set the shape correction: %p", c);
317   if (!c) return;
318   fShapeCorr = static_cast<TH2F*>(c->Clone());
319   fShapeCorr->SetDirectory(0);
320 }
321
322 //________________________________________________________________________
323 void 
324 AliBasedNdetaTask::InitializeCentBins()
325 {
326   if (fListOfCentralities->GetEntries() > 0) return;
327
328   // Automatically add 'all' centrality bin if nothing has been defined. 
329   AddCentralityBin(0, 0, 0);
330   if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) { 
331     const TArrayD* bins = fCentAxis->GetXbins();
332     Int_t          nbin = fCentAxis->GetNbins(); 
333     for (Int_t i = 0; i < nbin; i++) 
334       AddCentralityBin(i+1,  Short_t((*bins)[i]), Short_t((*bins)[i+1]));
335   }
336 }
337
338 //________________________________________________________________________
339 void 
340 AliBasedNdetaTask::UserCreateOutputObjects()
341 {
342   // 
343   // Create output objects.  
344   //
345   // This is called once per slave process 
346   //
347   DGUARD(fDebug,1,"Create user ouput object");
348   fSums = new TList;
349   fSums->SetName(Form("%s_sums", GetName()));
350   fSums->SetOwner();
351
352   InitializeCentBins();
353   if (fCentAxis) fSums->Add(fCentAxis);
354
355
356   // Centrality histogram 
357   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
358   fCent->SetDirectory(0);
359   fCent->SetXTitle(0);
360   fSums->Add(fCent);
361
362   // Loop over centrality bins 
363   TIter next(fListOfCentralities);
364   CentralityBin* bin = 0;
365   while ((bin = static_cast<CentralityBin*>(next()))) 
366     bin->CreateOutputObjects(fSums);
367   
368
369   // Post data for ALL output slots >0 here, to get at least an empty
370   // histogram
371   PostData(1, fSums); 
372 }
373   
374 //____________________________________________________________________
375 void 
376 AliBasedNdetaTask::UserExec(Option_t *) 
377 {
378   // 
379   // Process a single event 
380   // 
381   // Parameters:
382   //    option Not used
383   //
384   // Main loop
385   DGUARD(fDebug,1,"Analyse the AOD event");
386   AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
387   if (!aod) return;
388
389   
390   TObject* obj = aod->FindListObject("Forward");
391   if (!obj) { 
392     AliWarning("No forward object found");
393     return;
394   }
395   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
396   
397   // Fill centrality histogram 
398   Float_t cent = forward->GetCentrality();
399   fCent->Fill(cent);
400
401   // Get the histogram(s) 
402   TH2D* data   = GetHistogram(aod, false);
403   TH2D* dataMC = GetHistogram(aod, true);
404
405   Bool_t isZero = ((fNormalizationScheme & kZeroBin) &&
406                    !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0));
407
408
409   // Loop over centrality bins 
410   CentralityBin* allBin = 
411     static_cast<CentralityBin*>(fListOfCentralities->At(0));
412   allBin->ProcessEvent(forward, fTriggerMask, isZero, 
413                        fVtxMin, fVtxMax, data, dataMC);
414   
415   // Find this centrality bin 
416   if (fCentAxis && fCentAxis->GetNbins() > 0) {
417     Int_t          icent   = fCentAxis->FindBin(cent);
418     CentralityBin* thisBin = 0;
419     if (icent >= 1 && icent <= fCentAxis->GetNbins()) 
420       thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
421     if (thisBin)
422       thisBin->ProcessEvent(forward, fTriggerMask, isZero, fVtxMin, fVtxMax, 
423                             data, dataMC);
424   }
425
426   // Here, we get the update 
427   if (!fSNNString) { 
428     fSNNString = AliForwardUtil::MakeParameter("sNN", forward->GetSNN());
429     fSysString = AliForwardUtil::MakeParameter("sys", forward->GetSystem());
430
431     fSums->Add(fSNNString);
432     fSums->Add(fSysString);
433     fSums->Add(fSchemeString);
434     fSums->Add(fTriggerString);
435
436     // Print();
437   }
438   PostData(1, fSums);
439 }
440
441 //________________________________________________________________________
442 void 
443 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
444                                           const char* title, const char* ytitle)
445 {
446   // 
447   // Set histogram graphical options, etc. 
448   // 
449   // Parameters:
450   //    h       Histogram to modify
451   //    colour  Marker color 
452   //    marker  Marker style
453   //    title   Title of histogram
454   //    ytitle  Title on y-axis. 
455   //
456   h->SetTitle(title);
457   h->SetMarkerColor(colour);
458   h->SetMarkerStyle(marker);
459   h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1);
460   h->SetFillStyle(0);
461   h->SetYTitle(ytitle);
462   h->GetXaxis()->SetTitleFont(132);
463   h->GetXaxis()->SetLabelFont(132);
464   h->GetXaxis()->SetNdivisions(10);
465   h->GetYaxis()->SetTitleFont(132);
466   h->GetYaxis()->SetLabelFont(132);
467   h->GetYaxis()->SetNdivisions(10);
468   h->GetYaxis()->SetDecimals();
469   h->SetStats(0);
470 }
471
472 //________________________________________________________________________
473 void
474 AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm) 
475 {
476   // Normalize to the acceptance -
477   // dndeta->Divide(accNorm);
478   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { 
479     Double_t a = norm->GetBinContent(i);
480     for (Int_t j = 1; j <= copy->GetNbinsY(); j++) { 
481       if (a <= 0) { 
482         copy->SetBinContent(i,j,0);
483         copy->SetBinError(i,j,0);
484         continue;
485       }
486       Double_t c = copy->GetBinContent(i, j);
487       Double_t e = copy->GetBinError(i, j);
488       copy->SetBinContent(i, j, c / a);
489       copy->SetBinError(i, j, e / a);
490     }
491   }
492 }
493
494 //________________________________________________________________________
495 TH1D*
496 AliBasedNdetaTask::ProjectX(const TH2D* h, 
497                             const char* name,
498                             Int_t firstbin, 
499                             Int_t lastbin, 
500                             bool  useRoot,
501                             bool  corr,
502                             bool  error)
503 {
504   // 
505   // Project onto the X axis 
506   // 
507   // Parameters:
508   //    h         2D histogram 
509   //    name      New name 
510   //    firstbin  First bin to use 
511   //    lastbin   Last bin to use
512   //    error     Whether to calculate errors
513   // 
514   // Return:
515   //    Newly created histogram or null
516   //
517   if (!h) return 0;
518   if (useRoot) 
519     return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
520   
521   TAxis* xaxis = h->GetXaxis();
522   TAxis* yaxis = h->GetYaxis();
523   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(), 
524                           xaxis->GetXmin(), xaxis->GetXmax());
525   static_cast<const TAttLine*>(h)->Copy(*ret);
526   static_cast<const TAttFill*>(h)->Copy(*ret);
527   static_cast<const TAttMarker*>(h)->Copy(*ret);
528   ret->GetXaxis()->ImportAttributes(xaxis);
529
530   Int_t first = firstbin; 
531   Int_t last  = lastbin;
532   if (first < 0)                         first = 0;
533   else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
534   if      (last  < 0)                    last  = yaxis->GetNbins();
535   else if (last  >  yaxis->GetNbins()+1) last  = yaxis->GetNbins();
536   if (last-first < 0) { 
537     AliWarningGeneral("AliBasedNdetaTask", 
538                       Form("Nothing to project [%d,%d]", first, last));
539     return 0;
540     
541   }
542
543   // Loop over X bins 
544   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
545   Int_t ybins = (last-first+1);
546   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) { 
547     Double_t content = 0;
548     Double_t error2  = 0;
549     Int_t    nbins   = 0;
550     
551     
552     for (Int_t ybin = first; ybin <= last; ybin++) { 
553       Double_t c1 = h->GetCellContent(xbin, ybin);
554       Double_t e1 = h->GetCellError(xbin, ybin);
555
556       // Ignore empty bins 
557       if (c1 < 1e-12) continue;
558       if (e1 < 1e-12) {
559         if (error) continue; 
560         e1 = 1;
561       }
562
563       content    += c1;
564       error2     += e1*e1;
565       nbins++;
566     } // for (ybin)
567     if(content > 0 && nbins > 0) {
568       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
569 #if 0
570       AliWarningGeneral(ret->GetName(), 
571                         Form("factor @ %d is %d/%d -> %f", 
572                              xbin, ybins, nbins, factor));
573 #endif
574       if (error) { 
575         // Calculate weighted average
576         ret->SetBinContent(xbin, content * factor);
577         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
578       }
579       else 
580         ret->SetBinContent(xbin, factor * content);
581     }
582   } // for (xbin)
583   
584   return ret;
585 }
586   
587 //________________________________________________________________________
588 void 
589 AliBasedNdetaTask::Terminate(Option_t *) 
590 {
591   // 
592   // Called at end of event processing.. 
593   //
594   // This is called once in the master 
595   // 
596   // Parameters:
597   //    option Not used 
598   //
599   // Draw result to screen, or perform fitting, normalizations Called
600   // once at the end of the query
601   DGUARD(fDebug,1,"Process final merged results");
602
603   fSums = dynamic_cast<TList*> (GetOutputData(1));
604   if(!fSums) {
605     AliError("Could not retrieve TList fSums"); 
606     return; 
607   }
608   
609   fOutput = new TList;
610   fOutput->SetName(Form("%s_result", GetName()));
611   fOutput->SetOwner();
612   
613   fSNNString     = fSums->FindObject("sNN");
614   fSysString     = fSums->FindObject("sys");
615   fCentAxis      = static_cast<TAxis*>(fSums->FindObject("centAxis"));
616   fSchemeString  = fSums->FindObject("scheme");
617   fTriggerString = fSums->FindObject("trigger");
618
619   if(fSysString && fSNNString && 
620      fSysString->GetUniqueID() == AliForwardUtil::kPP)
621     LoadNormalizationData(fSysString->GetUniqueID(),
622                           fSNNString->GetUniqueID());
623
624   InitializeCentBins();
625
626   // Print before we loop
627   Print();
628
629   // Loop over centrality bins 
630   TIter next(fListOfCentralities);
631   CentralityBin* bin = 0;
632   gStyle->SetPalette(1);
633   THStack* dndetaStack        = new THStack("dndeta", "dN/d#eta");
634   THStack* dndetaStackRebin   = new THStack(Form("dndeta_rebin%02d", fRebin), 
635                                             "dN_{ch}/d#eta");
636   THStack* dndetaMCStack      = new THStack("dndetaMC", "dN_{ch}/d#eta");
637   THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin), 
638                                             "dN_{ch}/d#eta");
639   
640   TList* mclist = 0;
641   TList* truthlist = 0;
642   
643   if (fFinalMCCorrFile.Contains(".root")) {
644     TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
645     if(ftest) {
646       mclist    = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
647       truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
648     }
649     else 
650       AliWarning("MC analysis file invalid - no final MC correction possible");
651   }
652   Int_t style = GetMarker();
653   Int_t color = GetColor();
654   
655   AliInfo(Form("Marker style=%d, color=%d", style, color));
656   while ((bin = static_cast<CentralityBin*>(next()))) {
657     
658     bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, 
659              fTriggerEff, fTriggerEff0, 
660              fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges, 
661              fTriggerMask, style, color, mclist, truthlist);
662     if (fCentAxis && bin->IsAllBin()) continue;
663     TH1* dndeta      = bin->GetResult(0, false, "");
664     TH1* dndetaSym   = bin->GetResult(0, true,  "");
665     TH1* dndetaMC    = bin->GetResult(0, false, "MC");
666     TH1* dndetaMCSym = bin->GetResult(0, true,  "MC");
667     if (dndeta)      dndetaStack->Add(dndeta);
668     if (dndetaSym)   dndetaStack->Add(dndetaSym);
669     if (dndetaMC)    dndetaMCStack->Add(dndetaMC);
670     if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
671     if (fRebin > 1) { 
672       dndeta      = bin->GetResult(fRebin, false, "");
673       dndetaSym   = bin->GetResult(fRebin, true,  "");
674       dndetaMC    = bin->GetResult(fRebin, false, "MC");
675       dndetaMCSym = bin->GetResult(fRebin, true,  "MC");
676       if (dndeta)      dndetaStackRebin->Add(dndeta);
677       if (dndetaSym)   dndetaStackRebin->Add(dndetaSym);
678       if (dndetaMC)    dndetaMCStackRebin->Add(dndetaMC);
679       if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
680     }
681   }
682   // Output the stack
683   fOutput->Add(dndetaStack);
684
685   // If available output rebinned stack 
686   if (!dndetaStackRebin->GetHists() || 
687       dndetaStackRebin->GetHists()->GetEntries() <= 0) {
688     AliWarning("No rebinned histograms found");
689     delete dndetaStackRebin;
690     dndetaStackRebin = 0;
691   }
692   if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
693
694   // If available, output track-ref stack
695   if (!dndetaMCStack->GetHists() || 
696       dndetaMCStack->GetHists()->GetEntries() <= 0) {
697     AliWarning("No MC histograms found");
698     delete dndetaMCStack;
699     dndetaMCStack = 0;
700   }
701   if (dndetaMCStack) fOutput->Add(dndetaMCStack);
702
703   // If available, output rebinned track-ref stack
704   if (!dndetaMCStackRebin->GetHists() || 
705       dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
706     AliWarning("No rebinned MC histograms found");
707     delete dndetaMCStackRebin;
708     dndetaMCStackRebin = 0;
709   }
710   if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
711
712   // Output collision energy string 
713   if (fSNNString) { 
714     UShort_t sNN = fSNNString->GetUniqueID();
715     TNamed* sNNObj = new TNamed(fSNNString->GetName(),
716                                 AliForwardUtil::CenterOfMassEnergyString(sNN));
717     sNNObj->SetUniqueID(sNN);
718     fOutput->Add(sNNObj); // fSNNString->Clone());
719   }
720
721   // Output collision system string 
722   if (fSysString) { 
723     UShort_t sys = fSysString->GetUniqueID();
724     TNamed* sysObj = new TNamed(fSysString->GetName(),
725                                 AliForwardUtil::CollisionSystemString(sys));
726     sysObj->SetUniqueID(sys);
727     fOutput->Add(sysObj); // fSysString->Clone());
728   }
729
730   // Output centrality axis 
731   if (fCentAxis) fOutput->Add(fCentAxis);
732
733   // Output trigger string 
734   if (fTriggerString) { 
735     UShort_t mask = fTriggerString->GetUniqueID();
736     TNamed* maskObj = new TNamed(fTriggerString->GetName(), 
737                                  AliAODForwardMult::GetTriggerString(mask));
738     maskObj->SetUniqueID(mask);
739     fOutput->Add(maskObj); // fTriggerString->Clone());
740   }
741   
742   // Normalization string 
743   if (fSchemeString) {
744     UShort_t scheme = fSchemeString->GetUniqueID();
745     TNamed* schemeObj = new TNamed(fSchemeString->GetName(), 
746                                    NormalizationSchemeString(scheme));
747     schemeObj->SetUniqueID(scheme);
748     fOutput->Add(schemeObj); // fSchemeString->Clone());
749   }
750
751   // Output vertex axis 
752   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
753   vtxAxis->SetName("vtxAxis");
754   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
755   fOutput->Add(vtxAxis);
756
757   // Output trigger efficiency and shape correction 
758   fOutput->Add(AliForwardUtil::MakeParameter("triggerEff", fTriggerEff));
759   fOutput->Add(AliForwardUtil::MakeParameter("triggerEff0", fTriggerEff0));
760   if (fShapeCorr) fOutput->Add(fShapeCorr);
761
762   TNamed* options = new TNamed("options","");
763   TString str;
764   str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
765   str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
766   str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
767   options->SetTitle(str);
768   fOutput->Add(options);
769
770   PostData(2, fOutput);
771 }
772 //________________________________________________________________________
773 void
774 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
775 {
776   // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
777   DGUARD(fDebug,1,"Load the normalization data for sys=%d, energy=%d",
778          sys, energy);
779   TString type("pp");
780   TString snn("900");
781   if(energy == 7000) snn.Form("7000");
782   if(energy == 2750) snn.Form("2750"); 
783   
784   if(fShapeCorr &&  (fTriggerEff != 1)) {
785     AliInfo("Objects already set for normalization - no action taken"); 
786     return; 
787   }
788   
789   TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/"
790                                 "Normalization/normalizationHists_%s_%s.root",
791                                 type.Data(),snn.Data()));
792   if(!fin) {
793     AliWarning(Form("no file for normalization of %d/%d", sys, energy));
794     return;
795   }
796
797   // Shape correction
798   if ((fNormalizationScheme & kShape) && !fShapeCorr) {
799     TString trigName("All");
800     if (fTriggerMask == AliAODForwardMult::kInel || 
801         fTriggerMask == AliAODForwardMult::kNClusterGt0) 
802       trigName = "Inel";
803     else if (fTriggerMask == AliAODForwardMult::kNSD)
804       trigName = "NSD";
805     else if (fTriggerMask == AliAODForwardMult::kInelGt0)
806       trigName = "InelGt0";
807     else {
808       AliWarning(Form("Normalization for trigger %s not known, using all",
809                       AliAODForwardMult::GetTriggerString(fTriggerMask)));
810     }
811       
812     TString shapeCorName(Form("h%sNormalization", trigName.Data()));
813     TH2F*   shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
814     if (shapeCor) SetShapeCorrection(shapeCor);
815     else { 
816       AliWarning(Form("No shape correction found for %s", trigName.Data()));
817     }
818   }
819
820   // Trigger efficiency
821   TString effName(Form("%sTriggerEff", 
822                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
823                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
824                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
825                        "inelgt0" : "all"));
826
827   Double_t trigEff = 1; 
828   if (fNormalizationScheme & kTriggerEfficiency) { 
829     TObject* eff = fin->Get(effName);
830     if (eff) AliForwardUtil::GetParameter(eff, trigEff);
831   }
832   if (fTriggerEff != 1) SetTriggerEff(trigEff);
833   if (fTriggerEff < 0)  fTriggerEff = 1;
834
835   // Trigger efficiency
836   TString eff0Name(Form("%sTriggerEff0", 
837                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
838                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
839                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
840                        "inelgt0" : "all"));
841
842   Double_t trigEff0 = 1; 
843   if (fNormalizationScheme & kTriggerEfficiency) { 
844     TObject* eff = fin->Get(eff0Name);
845     if (eff) AliForwardUtil::GetParameter(eff, trigEff0);
846   }
847   if (fTriggerEff0 != 1) SetTriggerEff0(trigEff0);
848   if (fTriggerEff0 < 0)  fTriggerEff0 = 1;
849
850   // TEMPORARY FIX
851   // Rescale the shape correction by the trigger efficiency 
852   if (fShapeCorr) {
853     AliWarning(Form("Rescaling shape correction by trigger efficiency: "
854                     "1/E_X=1/%f", trigEff));
855     fShapeCorr->Scale(1. / trigEff);
856   }
857
858   // Print - out
859   if (fShapeCorr && fTriggerEff) AliInfo("Loaded objects for normalization.");
860 }
861
862
863 //________________________________________________________________________
864 void 
865 AliBasedNdetaTask::Print(Option_t*) const
866 {
867   // 
868   // Print information 
869   // 
870   std::cout << this->ClassName() << ": " << this->GetName() << "\n"
871             << std::boolalpha 
872             << " Trigger:                    " << (fTriggerString ? 
873                                                    fTriggerString->GetTitle() :
874                                                    "none") << "\n"
875             << " Vertex range:               [" << fVtxMin << ":" 
876             << fVtxMax << "]\n"
877             << " Rebin factor:               " << fRebin << "\n" 
878             << " Cut edges:                  " << fCutEdges << "\n"
879             << " Symmertrice:                " << fSymmetrice << "\n"
880             << " Use TH2::ProjectionX:       " << fUseROOTProj << "\n"
881             << " Correct for empty:          " << fCorrEmpty << "\n"
882             << " Normalization scheme:       " << (fSchemeString ? 
883                                                    fSchemeString->GetTitle() : 
884                                              "none") <<"\n"
885             << " Trigger efficiency:         " << fTriggerEff << "\n" 
886             << " Bin-0 Trigger efficiency:   " << fTriggerEff0 << "\n" 
887             << " Shape correction:           " << (fShapeCorr ? 
888                                                    fShapeCorr->GetName() : 
889                                                    "none") << "\n"
890             << " sqrt(s_NN):                 " << (fSNNString ? 
891                                                    fSNNString->GetTitle() : 
892                                                    "unknown") << "\n"
893             << " Collision system:           " << (fSysString ? 
894                                                    fSysString->GetTitle() : 
895                                                    "unknown") << "\n"
896             << " Centrality bins:            " << (fCentAxis ? "" : "none");
897   if (fCentAxis) { 
898     Int_t           nBins = fCentAxis->GetNbins();
899     const Double_t* bins  = fCentAxis->GetXbins()->GetArray();
900     for (Int_t i = 0; i <= nBins; i++) 
901       std::cout << (i==0 ? "" : "-") << bins[i];
902   }
903   std::cout << std::noboolalpha << std::endl;
904   
905 }
906
907 //________________________________________________________________________
908 TH1D*
909 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
910 {
911   // 
912   // Make a copy of the input histogram and rebin that histogram
913   // 
914   // Parameters:
915   //    h  Histogram to rebin
916   // 
917   // Return:
918   //    New (rebinned) histogram
919   //
920   if (rebin <= 1) return 0;
921
922   Int_t nBins = h->GetNbinsX();
923   if(nBins % rebin != 0) {
924     AliWarningGeneral("AliBasedNdetaTask", 
925                       Form("Rebin factor %d is not a devisor of current number "
926                            "of bins %d in the histogram %s", 
927                            rebin, nBins, h->GetName()));
928     return 0;
929   }
930     
931   // Make a copy 
932   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d", 
933                                                h->GetName(), rebin)));
934   tmp->Rebin(rebin);
935   tmp->SetDirectory(0);
936
937   // The new number of bins 
938   Int_t nBinsNew = nBins / rebin;
939   for(Int_t i = 1;i<= nBinsNew; i++) {
940     Double_t content = 0;
941     Double_t sumw    = 0;
942     Double_t wsum    = 0;
943     Int_t    nbins   = 0;
944     for(Int_t j = 1; j<=rebin;j++) {
945       Int_t    bin = (i-1)*rebin + j;
946       Double_t c   =  h->GetBinContent(bin);
947       if (c <= 0) continue;
948       
949       if (cutEdges) {
950         if (h->GetBinContent(bin+1)<=0 || 
951             h->GetBinContent(bin-1)<=0) {
952 #if 0
953           AliWarningGeneral("AliBasedNdetaTask", 
954                             Form("removing bin %d=%f of %s (%d=%f,%d=%f)", 
955                                  bin, c, h->GetName(), 
956                                  bin+1, h->GetBinContent(bin+1), 
957                                  bin-1, h->GetBinContent(bin-1)));
958 #endif
959           continue;
960         }       
961       }
962       Double_t e =  h->GetBinError(bin);
963       Double_t w =  1 / (e*e); // 1/c/c
964       content    += c;
965       sumw       += w;
966       wsum       += w * c;
967       nbins++;
968     }
969       
970     if(content > 0 && nbins > 0) {
971       tmp->SetBinContent(i, wsum / sumw);
972       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
973     }
974   }
975   
976   return tmp;
977 }
978
979 //__________________________________________________________________
980 TH1* 
981 AliBasedNdetaTask::Symmetrice(const TH1* h)
982 {
983   // 
984   // Make an extension of @a h to make it symmetric about 0 
985   // 
986   // Parameters:
987   //    h Histogram to symmertrice 
988   // 
989   // Return:
990   //    Symmetric extension of @a h 
991   //
992   Int_t nBins = h->GetNbinsX();
993   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
994   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
995   s->Reset();
996   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
997   s->SetMarkerColor(h->GetMarkerColor());
998   s->SetMarkerSize(h->GetMarkerSize());
999   s->SetMarkerStyle(FlipHollowStyle(h->GetMarkerStyle()));
1000   s->SetFillColor(h->GetFillColor());
1001   s->SetFillStyle(h->GetFillStyle());
1002   s->SetDirectory(0);
1003
1004   // Find the first and last bin with data 
1005   Int_t first = nBins+1;
1006   Int_t last  = 0;
1007   for (Int_t i = 1; i <= nBins; i++) { 
1008     if (h->GetBinContent(i) <= 0) continue; 
1009     first = TMath::Min(first, i);
1010     last  = TMath::Max(last,  i);
1011   }
1012     
1013   Double_t xfirst = h->GetBinCenter(first-1);
1014   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
1015   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
1016   for (Int_t i = f1, j=l2; i <= last; i++,j--) { 
1017     s->SetBinContent(j, h->GetBinContent(i));
1018     s->SetBinError(j, h->GetBinError(i));
1019   }
1020   // Fill in overlap bin 
1021   s->SetBinContent(l2+1, h->GetBinContent(first));
1022   s->SetBinError(l2+1, h->GetBinError(first));
1023   return s;
1024 }
1025
1026 //__________________________________________________________________
1027 Int_t
1028 AliBasedNdetaTask::GetMarkerStyle(UShort_t bits)
1029 {
1030   Int_t  base   = bits & (0xFE);
1031   Bool_t hollow = bits & kHollow;
1032   switch (base) { 
1033   case kCircle:       return (hollow ? 24 : 20);
1034   case kSquare:       return (hollow ? 25 : 21);
1035   case kUpTriangle:   return (hollow ? 26 : 22);
1036   case kDownTriangle: return (hollow ? 32 : 23);
1037   case kDiamond:      return (hollow ? 27 : 33); 
1038   case kCross:        return (hollow ? 28 : 34); 
1039   case kStar:         return (hollow ? 30 : 29); 
1040   }
1041   return 1;
1042 }
1043 //__________________________________________________________________
1044 UShort_t
1045 AliBasedNdetaTask::GetMarkerBits(Int_t style) 
1046
1047   UShort_t bits = 0;
1048   switch (style) { 
1049   case 24: case 25: case 26: case 27: case 28: case 30: case 32: 
1050     bits |= kHollow; break;
1051   }
1052   switch (style) { 
1053   case 20: case 24: bits |= kCircle;       break;
1054   case 21: case 25: bits |= kSquare;       break;
1055   case 22: case 26: bits |= kUpTriangle;   break;
1056   case 23: case 32: bits |= kDownTriangle; break;
1057   case 27: case 33: bits |= kDiamond;      break;
1058   case 28: case 34: bits |= kCross;        break;
1059   case 29: case 30: bits |= kStar;         break;
1060   }
1061   return bits;
1062 }
1063 //__________________________________________________________________
1064 Int_t
1065 AliBasedNdetaTask::FlipHollowStyle(Int_t style) 
1066 {
1067   UShort_t bits = GetMarkerBits(style);
1068   Int_t    ret  = GetMarkerStyle(bits ^ kHollow);
1069   return ret;
1070 }
1071
1072 //====================================================================
1073 void
1074 AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col)
1075 {
1076   DGUARD(fDebug,1,"Initializing sums with %s", data->GetName());
1077   TString n(GetHistName(0));
1078   TString n0(GetHistName(1));
1079   const char* postfix = GetTitle();
1080
1081   fSum = static_cast<TH2D*>(data->Clone(n));
1082   if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix));
1083   fSum->SetDirectory(0);
1084   fSum->SetMarkerColor(col);
1085   fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid));
1086   fSum->Reset();
1087   list->Add(fSum);
1088
1089   fSum0 = static_cast<TH2D*>(data->Clone(n0));
1090   if (postfix) 
1091     fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix));
1092   else   
1093     fSum0->SetTitle(Form("%s 0-bin", data->GetTitle()));
1094   fSum0->SetDirectory(0);
1095   fSum0->SetMarkerColor(col);
1096   fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow));
1097   fSum0->Reset();
1098   list->Add(fSum0);
1099
1100   fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5);
1101   fEvents->SetDirectory(0);
1102   fEvents->GetXaxis()->SetBinLabel(1, "Non-zero");
1103   fEvents->GetXaxis()->SetBinLabel(2, "Zero");
1104   list->Add(fEvents);
1105 }
1106
1107 //____________________________________________________________________
1108 TString
1109 AliBasedNdetaTask::Sum::GetHistName(const char* name, Int_t what, const char* post)
1110 {
1111   TString n(name);
1112   if      (what == 1) n.Append("0");
1113   else if (what == 2) n.Append("Events");
1114   if (post && post[0] != '\0')  n.Append(post);
1115   return n;
1116 }
1117
1118 //____________________________________________________________________
1119 TString
1120 AliBasedNdetaTask::Sum::GetHistName(Int_t what) const
1121 {
1122   return GetHistName(GetName(), what, GetTitle());
1123 }
1124
1125 //____________________________________________________________________
1126 void
1127 AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero)
1128 {
1129   DGUARD(fDebug,2,"Adding %s to sums", data->GetName());
1130   if (isZero) fSum0->Add(data);
1131   else        fSum->Add(data);
1132   fEvents->Fill(isZero ? 1 : 0);
1133 }
1134
1135 //____________________________________________________________________
1136 TH2D*
1137 AliBasedNdetaTask::Sum::CalcSum(TList*       output, 
1138                                 Double_t&    ntotal,
1139                                 Double_t     epsilon0, 
1140                                 Double_t     epsilon,
1141                                 Int_t        marker,
1142                                 Bool_t       rootProj, 
1143                                 Bool_t       corrEmpty) const
1144 {
1145   DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName());
1146
1147   TH2D* ret      = static_cast<TH2D*>(fSum->Clone(fSum->GetName()));
1148   ret->SetDirectory(0);
1149   ret->Reset();
1150   Int_t n        = Int_t(fEvents->GetBinContent(1));
1151   Int_t n0       = Int_t(fEvents->GetBinContent(2));
1152
1153   AliInfoF("Adding histograms %s and %s with weights %f and %f resp.",
1154            fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1155   DMSG(fDebug,2, "Adding histograms %s and %s with weights %f and %f resp.",
1156        fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1157   // Generate merged histogram 
1158   ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon); 
1159   ntotal = n / epsilon + n0 / epsilon0;
1160
1161   TList* out = new TList;
1162   out->SetOwner();
1163   const char* postfix = GetTitle();
1164   if (!postfix) postfix = "";
1165   out->SetName(Form("partial%s", postfix));
1166   output->Add(out);
1167
1168   // Now make copies, normalize them, and store in output list 
1169   TH2D* sumCopy  = static_cast<TH2D*>(fSum->Clone("sum"));
1170   TH2D* sum0Copy = static_cast<TH2D*>(fSum0->Clone("sum0"));
1171   TH2D* retCopy  = static_cast<TH2D*>(ret->Clone("sumAll"));
1172   sumCopy->SetMarkerStyle(FlipHollowStyle(marker));
1173   sumCopy->SetDirectory(0);
1174   sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4));
1175   sum0Copy->SetDirectory(0);
1176   retCopy->SetMarkerStyle(marker);
1177   retCopy->SetDirectory(0);
1178
1179   TH1D* norm    = ProjectX(fSum,  "norm",    0, 0, rootProj, corrEmpty, false);
1180   TH1D* norm0   = ProjectX(fSum0, "norm0",   0, 0, rootProj, corrEmpty, false);
1181   TH1D* normAll = ProjectX(ret,   "normAll", 0, 0, rootProj, corrEmpty, false);
1182   norm->SetDirectory(0);
1183   norm0->SetDirectory(0);
1184   normAll->SetDirectory(0);
1185   
1186   ScaleToCoverage(sumCopy, norm);
1187   ScaleToCoverage(sum0Copy, norm0);
1188   ScaleToCoverage(retCopy, normAll);
1189
1190   Int_t nY = fSum->GetNbinsY();
1191   TH1D* sumCopyPx  = ProjectX(sumCopy,  "average",    1, nY,rootProj,corrEmpty);
1192   TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0",   1, nY,rootProj,corrEmpty);
1193   TH1D* retCopyPx  = ProjectX(retCopy,  "averageAll", 1, nY,rootProj,corrEmpty);
1194   sumCopyPx->SetDirectory(0);
1195   sum0CopyPx->SetDirectory(0);
1196   retCopyPx->SetDirectory(0);
1197
1198   // Scale our 1D histograms
1199   sumCopyPx->Scale(1., "width");
1200   sum0CopyPx->Scale(1., "width");  
1201   retCopyPx->Scale(1., "width");  
1202
1203   AliInfo(Form("Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(), 
1204                sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()));
1205
1206   // Scale the normalization - they should be 1 at the maximum
1207   norm->Scale(n > 0   ? 1. / n  : 1);
1208   norm0->Scale(n0 > 0 ? 1. / n0 : 1);
1209   normAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1210
1211   out->Add(sumCopy);
1212   out->Add(sum0Copy);
1213   out->Add(retCopy);
1214   out->Add(sumCopyPx);
1215   out->Add(sum0CopyPx);
1216   out->Add(retCopyPx);
1217   out->Add(norm);
1218   out->Add(norm0);
1219   out->Add(normAll);
1220
1221   AliInfoF("Returning  (1/%f * %s + 1/%f * %s), "
1222            "1/%f * %d + 1/%f * %d = %d", 
1223            epsilon0, fSum0->GetName(), epsilon, fSum->GetName(), 
1224            epsilon0, n0, epsilon, n, int(ntotal));
1225 #if 0
1226   for (Int_t i = 1; i <= ret->GetNbinsX(); i++) { 
1227     Double_t nc  = sum->GetBinContent(i, 0);
1228     Double_t nc0 = sum0->GetBinContent(i, 0);
1229     ret->SetBinContent(i, 0, nc + nc0); // Just count events 
1230   }
1231 #endif
1232  
1233   return ret;
1234 }
1235
1236 //====================================================================
1237 AliBasedNdetaTask::CentralityBin::CentralityBin()
1238   : TNamed("", ""), 
1239     fSums(0), 
1240     fOutput(0),
1241     fSum(0), 
1242     fSumMC(0), 
1243     fTriggers(0), 
1244     fLow(0), 
1245     fHigh(0),
1246     fDoFinalMCCorrection(false), 
1247     fDebug(0)
1248 {
1249   // 
1250   // Constructor 
1251   //
1252   DGUARD(0,0,"Centrality bin default construction");
1253 }
1254 //____________________________________________________________________
1255 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name, 
1256                                                 Short_t low, Short_t high)
1257   : TNamed(name, ""), 
1258     fSums(0), 
1259     fOutput(0),
1260     fSum(0), 
1261     fSumMC(0), 
1262     fTriggers(0),
1263     fLow(low), 
1264     fHigh(high),
1265     fDoFinalMCCorrection(false), 
1266     fDebug(0)
1267 {
1268   // 
1269   // Constructor 
1270   // 
1271   // Parameters:
1272   //    name Name used for histograms (e.g., Forward)
1273   //    low  Lower centrality cut in percent 
1274   //    high Upper centrality cut in percent 
1275   //
1276   DGUARD(0,0,"Named Centrality bin construction: %s [%d,%d]", name);
1277   if (low <= 0 && high <= 0) { 
1278     fLow  = 0; 
1279     fHigh = 0;
1280     SetTitle("All centralities");
1281   }
1282   else {
1283     fLow  = low;
1284     fHigh = high;
1285     SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
1286   }
1287 }
1288 //____________________________________________________________________
1289 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
1290   : TNamed(o), 
1291     fSums(o.fSums), 
1292     fOutput(o.fOutput),
1293     fSum(o.fSum), 
1294     fSumMC(o.fSumMC), 
1295     fTriggers(o.fTriggers), 
1296     fLow(o.fLow), 
1297     fHigh(o.fHigh),
1298     fDoFinalMCCorrection(o.fDoFinalMCCorrection), 
1299     fDebug(o.fDebug)
1300 {
1301   // 
1302   // Copy constructor 
1303   // 
1304   // Parameters:
1305   //    other Object to copy from 
1306   //
1307   DGUARD(0,0,"Copy Centrality bin construction");
1308 }
1309 //____________________________________________________________________
1310 AliBasedNdetaTask::CentralityBin::~CentralityBin()
1311 {
1312   // 
1313   // Destructor 
1314   //
1315   DGUARD(fDebug,3,"Centrality bin desctruction");
1316   if (fSums) fSums->Delete();
1317   if (fOutput) fOutput->Delete();
1318 }
1319
1320 //____________________________________________________________________
1321 AliBasedNdetaTask::CentralityBin&
1322 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
1323 {
1324   // 
1325   // Assignment operator 
1326   // 
1327   // Parameters:
1328   //    other Object to assign from 
1329   // 
1330   // Return:
1331   //    Reference to this 
1332   //
1333   DGUARD(fDebug,3,"Centrality bin assignment");
1334   if (&o == this) return *this; 
1335   SetName(o.GetName());
1336   SetTitle(o.GetTitle());
1337   fSums      = o.fSums;
1338   fOutput    = o.fOutput;
1339   fSum       = o.fSum;
1340   fSumMC     = o.fSumMC;
1341   fTriggers  = o.fTriggers;
1342   fLow       = o.fLow;
1343   fHigh      = o.fHigh;
1344   fDoFinalMCCorrection = o.fDoFinalMCCorrection;
1345
1346   return *this;
1347 }
1348 //____________________________________________________________________
1349 Int_t
1350 AliBasedNdetaTask::CentralityBin::GetColor(Int_t fallback) const 
1351 {
1352   if (IsAllBin()) return fallback;
1353   Float_t  fc       = (fLow+double(fHigh-fLow)/2) / 100;
1354   Int_t    nCol     = gStyle->GetNumberOfColors();
1355   Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
1356   Int_t    col      = gStyle->GetColorPalette(icol);
1357   return col;
1358 }
1359 //____________________________________________________________________
1360 const char* 
1361 AliBasedNdetaTask::CentralityBin::GetListName() const
1362 {
1363   // 
1364   // Get the list name 
1365   // 
1366   // Return:
1367   //    List Name 
1368   //
1369   if (IsAllBin()) return "all"; 
1370   return Form("cent%03d_%03d", fLow, fHigh);
1371 }
1372 //____________________________________________________________________
1373 void
1374 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
1375 {
1376   // 
1377   // Create output objects 
1378   // 
1379   // Parameters:
1380   //    dir   Parent list
1381   //
1382   DGUARD(fDebug,1,"Create centrality bin output objects");
1383   fSums = new TList;
1384   fSums->SetName(GetListName());
1385   fSums->SetOwner();
1386   dir->Add(fSums);
1387
1388   fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers");
1389   fTriggers->SetDirectory(0);
1390   fSums->Add(fTriggers);
1391 }
1392 //____________________________________________________________________
1393 void
1394 AliBasedNdetaTask::CentralityBin::SetDebugLevel(Int_t lvl)
1395 {
1396   fDebug = lvl;
1397   if (fSum) fSum->fDebug = lvl;
1398   if (fSumMC) fSumMC->fDebug = lvl;
1399 }
1400
1401 //____________________________________________________________________
1402 Bool_t
1403 AliBasedNdetaTask::CentralityBin::ReadSum(TList* list, bool mc)
1404 {
1405   const char* post = (mc ? "MC" : "");
1406   TString     sn   = Sum::GetHistName(GetName(),0,post);
1407   TString     sn0  = Sum::GetHistName(GetName(),1,post);
1408   TString     ev   = Sum::GetHistName(GetName(),2,post);
1409   TH2D* sum      = static_cast<TH2D*>(list->FindObject(sn));
1410   TH2D* sum0     = static_cast<TH2D*>(list->FindObject(sn0));
1411   TH1I* events   = static_cast<TH1I*>(list->FindObject(ev));
1412   if (!sum || !sum0 || !events) {
1413     AliWarningF("Failed to find one or more histograms: "
1414                 "%s (%p) %s (%p) %s (%p)", 
1415                 sn.Data(), sum, 
1416                 sn0.Data(), sum0, 
1417                 ev.Data(), events); 
1418     return false;
1419   }
1420   Sum* ret     = new Sum(GetName(), post);
1421   ret->fSum    = sum;
1422   ret->fSum0   = sum0;
1423   ret->fEvents = events;
1424   ret->fDebug  = fDebug;
1425   if (mc) fSumMC = ret;
1426   else    fSum   = ret;
1427
1428   return true;
1429 }
1430
1431 //____________________________________________________________________
1432 void
1433 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
1434 {
1435   // 
1436   // Create sum histogram 
1437   // 
1438   // Parameters:
1439   //    data  Data histogram to clone 
1440   //    mc    (optional) MC histogram to clone 
1441   //
1442   DGUARD(fDebug,1,"Create centrality bin sums from %s", 
1443          data ? data->GetName() : "(null)");
1444   if (data) {
1445     fSum = new Sum(GetName(),"");
1446     fSum->Init(fSums, data, GetColor());
1447     fSum->fDebug = fDebug;
1448   }
1449
1450   // If no MC data is given, then do not create MC sum histogram 
1451   if (!mc) return;
1452
1453   fSumMC = new Sum(GetName(), "MC");
1454   fSumMC->Init(fSums, mc, GetColor());
1455   fSumMC->fDebug = fDebug;
1456 }
1457
1458 //____________________________________________________________________
1459 Bool_t
1460 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
1461                                              Int_t triggerMask,
1462                                              Double_t vzMin, Double_t vzMax)
1463 {
1464   // 
1465   // Check the trigger, vertex, and centrality
1466   // 
1467   // Parameters:
1468   //    aod Event input 
1469   // 
1470   // Return:
1471   //    true if the event is to be used 
1472   //
1473   if (!forward) return false;
1474
1475   DGUARD(fDebug,2,"Check the event");
1476   // We do not check for centrality here - it's already done 
1477   return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, fTriggers);
1478 }
1479   
1480   
1481 //____________________________________________________________________
1482 void
1483 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
1484                                                Int_t triggerMask, Bool_t isZero,
1485                                                Double_t vzMin, Double_t vzMax,
1486                                                const TH2D* data, const TH2D* mc)
1487 {
1488   // 
1489   // Process an event
1490   // 
1491   // Parameters:
1492   //    forward     Forward data (for trigger, vertex, & centrality)
1493   //    triggerMask Trigger mask 
1494   //    vzMin       Minimum IP z coordinate
1495   //    vzMax       Maximum IP z coordinate
1496   //    data        Data histogram 
1497   //    mc          MC histogram
1498   //
1499   DGUARD(fDebug,1,"Process one event for %s a given centrality bin", 
1500          data ? data->GetName() : "(null)");
1501   if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
1502   if (!data) return;
1503   if (!fSum) CreateSums(data, mc);
1504
1505   fSum->Add(data, isZero);
1506   if (mc) fSumMC->Add(mc, isZero);
1507 }
1508
1509 //________________________________________________________________________
1510 Double_t 
1511 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1512                                                 UShort_t    scheme,
1513                                                 Double_t    trigEff,
1514                                                 Double_t&   ntotal) const
1515 {
1516   // 
1517   // Calculate normalization 
1518   // 
1519   // Parameters: 
1520   //    t       Trigger histogram
1521   //    scheme  Normaliztion scheme 
1522   //    trigEff From MC
1523   //    ntotal  On return, contains the number of events. 
1524   //
1525   DGUARD(fDebug,1,"Normalize centrality bin with %s", t.GetName());
1526   Double_t nAll        = t.GetBinContent(AliAODForwardMult::kBinAll);
1527   Double_t nB          = t.GetBinContent(AliAODForwardMult::kBinB);
1528   Double_t nA          = t.GetBinContent(AliAODForwardMult::kBinA);
1529   Double_t nC          = t.GetBinContent(AliAODForwardMult::kBinC);
1530   Double_t nE          = t.GetBinContent(AliAODForwardMult::kBinE);
1531   Double_t nOffline    = t.GetBinContent(AliAODForwardMult::kBinOffline);
1532   Double_t nTriggered  = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1533   Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1534   Double_t nAccepted   = ntotal; 
1535   ntotal               = 0;
1536   
1537   if (nTriggered <= 0.1) { 
1538     AliError("Number of triggered events <= 0");
1539     return -1;
1540   }
1541   if (nWithVertex <= 0.1) { 
1542     AliError("Number of events with vertex <= 0");
1543     return -1;
1544   }
1545   ntotal          = nAccepted;
1546   Double_t vtxEff = nWithVertex / nTriggered;
1547   Double_t scaler = 1;
1548   Double_t beta   = nA + nC - 2*nE;
1549
1550   if (scheme & kEventLevel && !(scheme & kZeroBin)) {
1551     ntotal = nAccepted / vtxEff;
1552     scaler = vtxEff;
1553     AliInfo(Form("Calculating event normalisation as\n"
1554                  " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1555                  Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex), 
1556                  ntotal, scaler));
1557             
1558     if (scheme & kBackground) {
1559       //          1            E_V             E_V
1560       //   s = --------- = ------------- = ------------ 
1561       //        1 - beta   1 - beta E_V    1 - beta N_V 
1562       //       ---  ----       --------        ---- ---
1563       //       E_V  N_V          N_V           N_V  N_T 
1564       // 
1565       //          E_V
1566       //     = ------------
1567       //        1 - beta 
1568       //            ----
1569       //             N_T 
1570       // 
1571       ntotal -= nAccepted * beta / nWithVertex;
1572       // This one is direct and correct. 
1573       // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1574       // A simpler expresion
1575       scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1576       AliInfo(Form("Calculating event normalisation as\n"
1577                    " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1578                    " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1579                    Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1580                    nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta), 
1581                    Int_t(nWithVertex), ntotal, scaler));
1582     }
1583   }
1584   if (scheme & kZeroBin) {
1585     // Calculate as 
1586     // 
1587     //  N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta)
1588     //    = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1589     // 
1590     //  s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1591     //    = N_V / (N_V + 1/E_X (N_T - N_V - beta)) 
1592     // 
1593     if (!(scheme & kBackground)) beta = 0;
1594     ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1 
1595                                          - beta / nWithVertex));
1596     scaler = nWithVertex / (nWithVertex + 
1597                             1/trigEff * (nTriggered-nWithVertex-beta));
1598     AliInfo(Form("Calculating event normalisation as\n"
1599                  "  beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1600                  "  N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = "
1601                  "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)",
1602                  Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1603                  Int_t(nAccepted), trigEff, Int_t(nTriggered), 
1604                  Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex), 
1605                  ntotal, scaler));
1606   }
1607   if (scheme & kTriggerEfficiency && !(scheme & kZeroBin)) {
1608     ntotal /= trigEff;
1609     scaler *= trigEff;
1610     AliInfo(Form("Correcting for trigger efficiency:\n"
1611                  " N = 1 / E_X * N = 1 / %f * %d = %f (%f)", 
1612                  trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1613   }
1614
1615   AliInfo(Form("\n"
1616                " Total of        %9d events for %s\n"
1617                "  of these       %9d have an offline trigger\n"
1618                "  of these N_T = %9d has the selected trigger\n" 
1619                "  of these N_V = %9d has a vertex\n" 
1620                "  of these N_A = %9d were in the selected range\n"
1621                "  Triggers by hardware type:\n"
1622                "    N_b   =  %9d\n"
1623                "    N_ac  =  %9d (%d+%d)\n"
1624                "    N_e   =  %9d\n"
1625                "  Vertex efficiency:          %f\n"
1626                "  Trigger efficiency:         %f\n"
1627                "  Total number of events: N = %f\n"
1628                "  Scaler (N_A/N):             %f",
1629                Int_t(nAll), GetTitle(), Int_t(nOffline), 
1630                Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1631                Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE), 
1632                vtxEff, trigEff, ntotal, scaler));
1633   return scaler;
1634 }
1635
1636 //________________________________________________________________________
1637 const char* 
1638 AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
1639                                                 Bool_t sym, 
1640                                                 const char* postfix) const
1641 {
1642   static TString n;
1643   n = Form("dndeta%s%s",GetName(), postfix);
1644   if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
1645   if (sym)       n.Append("_mirror");
1646   return n.Data();
1647 }
1648 //________________________________________________________________________
1649 TH1* 
1650 AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
1651                                             Bool_t sym, 
1652                                             const char* postfix) const
1653 {
1654   if (!fOutput) { 
1655     AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(), 
1656                     fLow, fHigh));
1657     return 0;
1658   }
1659   TString  n = GetResultName(rebin, sym, postfix);
1660   TObject* o = fOutput->FindObject(n.Data());
1661   if (!o) { 
1662     // AliWarning(Form("Object %s not found in output list", n.Data()));
1663     return 0;
1664   }
1665   return static_cast<TH1*>(o);
1666 }
1667
1668 //________________________________________________________________________
1669 void 
1670 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,  
1671                                              const char* postfix, 
1672                                              bool        rootProj, 
1673                                              bool        corrEmpty,
1674                                              const TH2F* shapeCorr,
1675                                              Double_t    scaler,
1676                                              bool        symmetrice, 
1677                                              Int_t       rebin, 
1678                                              bool        cutEdges, 
1679                                              Int_t       marker,
1680                                              Int_t       color,
1681                                              TList*      mclist, 
1682                                              TList*      truthlist)
1683 {
1684   // 
1685   // Generate the dN/deta result from input 
1686   // 
1687   // Parameters: 
1688   //     sum        Sum of 2D hists 
1689   //     postfix    Post fix on names
1690   //     rootProj   Whether to use ROOT TH2::ProjectionX
1691   //     corrEmpty  Correct for empty bins 
1692   //     shapeCorr  Shape correction to use 
1693   //     scaler     Event-level normalization scaler  
1694   //     symmetrice Whether to make symmetric extensions 
1695   //     rebin      Whether to rebin
1696   //     cutEdges   Whether to cut edges when rebinning 
1697   //
1698   DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
1699   TH2D* copy    = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s", 
1700                                                      GetName(), postfix)));
1701   TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0, 
1702                            rootProj, corrEmpty, false);
1703   accNorm->SetDirectory(0);
1704
1705   // ---- Scale by shape correction ----------------------------------
1706   if (shapeCorr) copy->Divide(shapeCorr);
1707   else AliInfo("No shape correction specified, or disabled");
1708   
1709   // --- Normalize to the coverage -----------------------------------
1710   ScaleToCoverage(copy, accNorm);
1711
1712   // --- Event-level normalization -----------------------------------
1713   copy->Scale(scaler);
1714
1715   // --- Project on X axis -------------------------------------------
1716   TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix), 
1717                           1, sum->GetNbinsY(), rootProj, corrEmpty);
1718   dndeta->SetDirectory(0);
1719   // Event-level normalization 
1720   dndeta->Scale(1., "width");
1721   copy->Scale(1., "width");
1722   
1723   TH1D* dndetaMCCorrection = 0;
1724   TList* centlist          = 0;
1725   TH1D* dndetaMCtruth      = 0;
1726   TList* truthcentlist     = 0;
1727   
1728   // Possible final correction to <MC analysis> / <MC truth>
1729   if(mclist) 
1730     centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
1731   if(centlist)
1732     dndetaMCCorrection = 
1733       static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s",
1734                                                    GetName(), postfix)));
1735   if(truthlist) 
1736     truthcentlist =static_cast<TList*>(truthlist->FindObject(GetListName()));
1737   if(truthcentlist)
1738     dndetaMCtruth =static_cast<TH1D*>(truthcentlist->FindObject("dndetaTruth"));
1739
1740   if(dndetaMCCorrection && dndetaMCtruth) {
1741     AliInfo("Correcting with final MC correction");
1742     dndetaMCCorrection->Divide(dndetaMCtruth);
1743     dndetaMCCorrection->SetTitle("Final MC correction");
1744     dndetaMCCorrection->SetName("finalMCCorr");
1745     dndeta->Divide(dndetaMCCorrection);    
1746   }
1747   else 
1748     AliInfo("No final MC correction applied");
1749   
1750   // --- Set some histogram attributes -------------------------------
1751   TString post;
1752   Int_t rColor = GetColor(color);
1753   if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
1754   SetHistogramAttributes(dndeta,  rColor, marker, 
1755                          Form("ALICE %s%s", GetName(), post.Data()));
1756   SetHistogramAttributes(accNorm, rColor, marker, 
1757                          Form("ALICE %s normalisation%s", 
1758                               GetName(), post.Data())); 
1759
1760   // --- Make symmetric extensions and rebinnings --------------------
1761   if (symmetrice)   fOutput->Add(Symmetrice(dndeta));
1762   fOutput->Add(dndeta);
1763   fOutput->Add(accNorm);
1764   fOutput->Add(copy);
1765   fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1766   if (symmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1767   if (dndetaMCCorrection) fOutput->Add(dndetaMCCorrection);
1768 }  
1769
1770 //________________________________________________________________________
1771 void 
1772 AliBasedNdetaTask::CentralityBin::End(TList*      sums, 
1773                                       TList*      results,
1774                                       UShort_t    scheme,
1775                                       const TH2F* shapeCorr, 
1776                                       Double_t    trigEff,
1777                                       Double_t    trigEff0,
1778                                       Bool_t      symmetrice,
1779                                       Int_t       rebin, 
1780                                       Bool_t      rootProj,
1781                                       Bool_t      corrEmpty, 
1782                                       Bool_t      cutEdges, 
1783                                       Int_t       triggerMask,
1784                                       Int_t       marker,
1785                                       Int_t       color, 
1786                                       TList*      mclist,
1787                                       TList*      truthlist) 
1788 {
1789   // 
1790   // End of processing 
1791   // 
1792   // Parameters:
1793   //    sums        List of sums
1794   //    results     Output list of results
1795   //    shapeCorr   Shape correction or nil
1796   //    trigEff     Trigger efficiency 
1797   //    symmetrice  Whether to symmetrice the results
1798   //    rebin       Whether to rebin the results
1799   //    corrEmpty   Whether to correct for empty bins
1800   //    cutEdges    Whether to cut edges when rebinning
1801   //    triggerMask Trigger mask 
1802   //
1803   DGUARD(fDebug,1,"End centrality bin procesing");
1804
1805   fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
1806   if(!fSums) {
1807     AliError("Could not retrieve TList fSums"); 
1808     return; 
1809   }
1810   
1811   fOutput = new TList;
1812   fOutput->SetName(GetListName());
1813   fOutput->SetOwner();
1814   results->Add(fOutput);
1815
1816   if (!fSum) { 
1817     if (!ReadSum(fSums, false)) {
1818       AliInfo("This task did not produce any output");
1819       return;
1820     }
1821   }
1822   if (!fSumMC) ReadSum(fSums, true);
1823
1824   fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
1825
1826   if (!fTriggers) { 
1827     AliError("Couldn't find histogram 'triggers' in list");
1828     return;
1829   }
1830
1831   // --- Get normalization scaler ------------------------------------
1832   Double_t epsilonT  = trigEff;
1833   Double_t epsilonT0 = trigEff0;
1834   AliInfoF("Using epsilonT=%f, epsilonT0=%f for %d", 
1835            epsilonT, epsilonT0, triggerMask);
1836 #if 0
1837   // TEMPORARY FIX
1838   if (triggerMask == AliAODForwardMult::kNSD) {
1839     // This is a local change 
1840     epsilonT = 0.96; 
1841     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
1842   }
1843   else if (triggerMask == AliAODForwardMult::kInel) {
1844     // This is a local change 
1845     epsilonT = 0.934; 
1846     AliWarning(Form("Using hard-coded Inel trigger efficiency of %f",epsilonT));
1847   }
1848   if (scheme & kZeroBin) { 
1849     if (triggerMask==AliAODForwardMult::kInel)
1850       epsilonT0 = 0.785021; // 0.100240;
1851     else if (triggerMask==AliAODForwardMult::kInelGt0)
1852       epsilonT0 = 0;
1853     else if (triggerMask==AliAODForwardMult::kNSD)
1854       epsilonT0 = .706587;
1855     epsilonT = 1;
1856     AliWarning(Form("Using hard-coded NCluster>0 trigger efficiency of %f",
1857                     epsilonT0));
1858   }
1859 #endif
1860
1861   // Get our histograms 
1862   Double_t nSum   = 0;
1863   TH2D*    sum    = fSum->CalcSum(fOutput, nSum, epsilonT0, epsilonT, 
1864                                  marker, rootProj, corrEmpty);
1865   Double_t nSumMC = 0;
1866   TH2D*    sumMC  = 0;
1867   if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC, 
1868                                       epsilonT0, epsilonT, marker,
1869                                       rootProj, corrEmpty);
1870   if (!sum) { 
1871     AliError("Failed to get sum from summer - bailing out");
1872     return;
1873   }
1874     
1875   Double_t ntotal = nSum;
1876   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
1877   if (scaler < 0) { 
1878     AliError("Failed to calculate normalization - bailing out");
1879     return;
1880   }
1881   fOutput->Add(fTriggers->Clone());
1882
1883   // --- Make result and store ---------------------------------------
1884   MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
1885              scaler, symmetrice, rebin, cutEdges, marker, color, 
1886              mclist, truthlist);
1887
1888   // --- Process result from TrackRefs -------------------------------
1889   if (sumMC) 
1890     MakeResult(sumMC, "MC", rootProj, corrEmpty, 
1891                (scheme & kShape) ? shapeCorr : 0,
1892                scaler, symmetrice, rebin, cutEdges, 
1893                GetMarkerStyle(GetMarkerBits(marker)+4), color, 
1894                mclist, truthlist);
1895   
1896   // Temporary stuff 
1897   // if (!IsAllBin()) return;
1898
1899 }
1900
1901 //
1902 // EOF
1903 //