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