Fixes
[u/mrichter/AliRoot.git] / PWG2 / 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
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     fShapeCorr(0),
32     fListOfCentralities(0),
33     fSNNString(0),
34     fSysString(0),
35     fCent(0),
36     fCentAxis(0),
37     fNormalizationScheme(kFull), 
38     fSchemeString(0), 
39     fTriggerString(0)
40 {
41   // 
42   // Constructor
43   // 
44 }
45
46 //____________________________________________________________________
47 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
48   : AliAnalysisTaskSE(name), 
49     fSums(0),           // Container of sums 
50     fOutput(0),         // Container of output 
51     fVtxMin(-10),       // Minimum v_z
52     fVtxMax(10),        // Maximum v_z
53     fTriggerMask(AliAODForwardMult::kInel), 
54     fRebin(5),          // Rebinning factor 
55     fCutEdges(false), 
56     fSymmetrice(true),
57     fCorrEmpty(true), 
58     fUseROOTProj(false),
59     fTriggerEff(1),
60     fShapeCorr(0),
61     fListOfCentralities(0),
62     fSNNString(0),
63     fSysString(0),
64     fCent(0),
65     fCentAxis(0),
66     fNormalizationScheme(kFull), 
67     fSchemeString(0),
68     fTriggerString(0)
69 {
70   // 
71   // Constructor
72   // 
73   fListOfCentralities = new TObjArray(1);
74   
75   // Set the normalisation scheme 
76   SetNormalizationScheme(kFull);
77
78   // Set the trigger mask
79   SetTriggerMask(AliAODForwardMult::kInel);
80
81   // Output slot #1 writes into a TH1 container
82   DefineOutput(1, TList::Class()); 
83   DefineOutput(2, TList::Class()); 
84 }
85
86 //____________________________________________________________________
87 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
88   : AliAnalysisTaskSE(o), 
89     fSums(o.fSums),             // TList* - Container of sums 
90     fOutput(o.fOutput),         // Container of output 
91     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
92     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
93     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask 
94     fRebin(o.fRebin),           // Int_t - Rebinning factor 
95     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
96     fSymmetrice(o.fSymmetrice),
97     fCorrEmpty(o.fCorrEmpty), 
98     fUseROOTProj(o.fUseROOTProj),
99     fTriggerEff(o.fTriggerEff),
100     fShapeCorr(o.fShapeCorr),
101     fListOfCentralities(o.fListOfCentralities),
102     fSNNString(o.fSNNString),
103     fSysString(o.fSysString),
104     fCent(o.fCent),
105     fCentAxis(o.fCentAxis),
106     fNormalizationScheme(o.fNormalizationScheme), 
107     fSchemeString(o.fSchemeString), 
108     fTriggerString(o.fTriggerString)
109 {}
110
111 //____________________________________________________________________
112 AliBasedNdetaTask::~AliBasedNdetaTask()
113 {
114   // 
115   // Destructor
116   // 
117   if (fSums) { 
118     fSums->Delete();
119     delete fSums;
120     fSums = 0;
121   }
122   if (fOutput) { 
123     fOutput->Delete();
124     delete fOutput;
125     fOutput = 0;
126   }
127 }
128
129 //________________________________________________________________________
130 void 
131 AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
132 {
133   if (!fCentAxis) { 
134     fCentAxis = new TAxis();
135     fCentAxis->SetName("centAxis");
136     fCentAxis->SetTitle("Centrality [%]");
137   }
138   TArrayD dbins(n+1);
139   for (UShort_t i = 0; i <= n; i++) 
140     dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
141   fCentAxis->Set(n, dbins.GetArray());
142 }
143     
144 //________________________________________________________________________
145 void 
146 AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
147 {
148   // 
149   // Add a centrality bin 
150   // 
151   // Parameters:
152   //    low  Low cut
153   //    high High cut
154   //
155   CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
156   AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at));
157   fListOfCentralities->AddAtAndExpand(bin, at);
158 }
159
160 //________________________________________________________________________
161 AliBasedNdetaTask::CentralityBin*
162 AliBasedNdetaTask::MakeCentralityBin(const char* name, 
163                                      Short_t low, Short_t high) const
164 {
165   // 
166   // Make a centrality bin 
167   // 
168   // Parameters:
169   //    name  Name used for histograms
170   //    low   Low cut in percent
171   //    high  High cut in percent
172   // 
173   // Return:
174   //    A newly created centrality bin 
175   //
176   return new CentralityBin(name, low, high);
177 }
178 //________________________________________________________________________
179 void 
180 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
181 {
182   // 
183   // Set normalisation scheme 
184   // 
185   UShort_t    scheme = 0;
186   TString     twhat(what);
187   twhat.ToUpper();
188   TObjString* opt;
189   TIter       next(twhat.Tokenize(" ,|"));
190   while ((opt = static_cast<TObjString*>(next()))) { 
191     TString s(opt->GetString());
192     if      (s.IsNull()) continue;
193     Bool_t add = true;
194     switch (s[0]) { 
195     case '-': add = false; // Fall through 
196     case '+': s.Remove(0,1);  // Remove character 
197     }
198     UShort_t bit = 0;
199     if      (s.CompareTo("EVENT")     == 0) bit = kEventLevel;
200     else if (s.CompareTo("SHAPE")     == 0) bit = kShape;
201     else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
202     else if (s.CompareTo("TRIGGER")   == 0) bit = kTriggerEfficiency;
203     else if (s.CompareTo("FULL")      == 0) bit = kFull;
204     else if (s.CompareTo("NONE")      == 0) bit = kNone;
205     else 
206       Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
207     if (add) scheme |= bit;
208     else     scheme ^= bit;
209   }
210   SetNormalizationScheme(scheme);
211 }
212 //________________________________________________________________________
213 void 
214 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme) 
215 {
216   fNormalizationScheme = scheme; 
217   TString tit = "";
218   if (scheme == kFull) tit = "FULL"; 
219   else {
220     if (scheme & kEventLevel)        tit.Append("EVENT ");
221     if (scheme & kShape)             tit.Append("SHAPE ");
222     if (scheme & kBackground)        tit.Append("BACKGROUND ");
223     if (scheme & kTriggerEfficiency) tit.Append("TRIGGER");
224   }
225   tit = tit.Strip(TString::kBoth);
226   if (!fSchemeString) fSchemeString = new TNamed("scheme", "");
227   fSchemeString->SetTitle(tit);
228   fSchemeString->SetUniqueID(fNormalizationScheme);
229 }
230 //________________________________________________________________________
231 void 
232 AliBasedNdetaTask::SetTriggerMask(const char* mask)
233 {
234   // 
235   // Set the trigger maskl 
236   // 
237   // Parameters:
238   //    mask Trigger mask
239   //
240   SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
241 }
242 //________________________________________________________________________
243 void 
244 AliBasedNdetaTask::SetTriggerMask(UShort_t mask) 
245
246   fTriggerMask = mask; 
247   TString tit(AliAODForwardMult::GetTriggerString(mask));
248   tit = tit.Strip(TString::kBoth);
249   if (!fTriggerString) fTriggerString = new TNamed("trigger", "");
250   fTriggerString->SetTitle(tit);
251   fTriggerString->SetUniqueID(fTriggerMask);
252 }
253
254 //________________________________________________________________________
255 void 
256 AliBasedNdetaTask::SetShapeCorrection(const TH1* c)
257 {
258   // 
259   // Set the shape correction (a.k.a., track correction) for selected
260   // trigger(s)
261   // 
262   // Parameters:
263   //    h Correction
264   //
265   if (!c) return;
266   fShapeCorr = static_cast<TH1*>(c->Clone());
267   fShapeCorr->SetDirectory(0);
268 }
269
270 //________________________________________________________________________
271 void 
272 AliBasedNdetaTask::UserCreateOutputObjects()
273 {
274   // 
275   // Create output objects.  
276   //
277   // This is called once per slave process 
278   //
279   fSums = new TList;
280   fSums->SetName(Form("%s_sums", GetName()));
281   fSums->SetOwner();
282
283   // Automatically add 'all' centrality bin if nothing has been defined. 
284   AddCentralityBin(0, 0, 0);
285   if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) { 
286     const TArrayD* bins = fCentAxis->GetXbins();
287     Int_t          nbin = fCentAxis->GetNbins(); 
288     for (Int_t i = 0; i < nbin; i++) 
289       AddCentralityBin(i+1,  Short_t((*bins)[i]), Short_t((*bins)[i+1]));
290   }
291   fSums->Add(fCentAxis);
292
293
294   // Centrality histogram 
295   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
296   fCent->SetDirectory(0);
297   fCent->SetXTitle(0);
298   fSums->Add(fCent);
299
300   // Loop over centrality bins 
301   TIter next(fListOfCentralities);
302   CentralityBin* bin = 0;
303   while ((bin = static_cast<CentralityBin*>(next()))) 
304     bin->CreateOutputObjects(fSums);
305
306   // Check that we have an AOD input handler 
307   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
308   AliAODInputHandler* ah = 
309     dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
310   if (!ah) AliFatal("No AOD input handler set in analysis manager");
311
312   // Post data for ALL output slots >0 here, to get at least an empty histogram
313   PostData(1, fSums); 
314 }
315 //____________________________________________________________________
316 void 
317 AliBasedNdetaTask::UserExec(Option_t *) 
318 {
319   // 
320   // Process a single event 
321   // 
322   // Parameters:
323   //    option Not used
324   //
325   // Main loop
326   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
327   if (!aod) {
328     AliError("Cannot get the AOD event");
329     return;
330   }  
331   
332   TObject* obj = aod->FindListObject("Forward");
333   if (!obj) { 
334     AliWarning("No forward object found");
335     return;
336   }
337   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
338   
339   // Fill centrality histogram 
340   Float_t cent = forward->GetCentrality();
341   fCent->Fill(cent);
342
343   // Get the histogram(s) 
344   TH2D* data   = GetHistogram(aod, false);
345   TH2D* dataMC = GetHistogram(aod, true);
346
347   // Loop over centrality bins 
348   CentralityBin* allBin = 
349     static_cast<CentralityBin*>(fListOfCentralities->At(0));
350   allBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
351   
352   // Find this centrality bin 
353   if (fCentAxis && fCentAxis->GetNbins() > 0) {
354     Int_t          icent   = fCentAxis->FindBin(cent);
355     CentralityBin* thisBin = 0;
356     if (icent >= 1 && icent <= fCentAxis->GetNbins()) 
357       thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
358     if (thisBin)
359       thisBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, 
360                             data, dataMC);
361   }
362
363   // Here, we get the update 
364   if (!fSNNString) { 
365     UShort_t sNN = forward->GetSNN();
366     fSNNString = new TNamed("sNN", "");
367     fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
368     fSNNString->SetUniqueID(sNN);
369     fSums->Add(fSNNString);
370   
371     UShort_t sys = forward->GetSystem();
372     fSysString = new TNamed("sys", "");
373     fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
374     fSysString->SetUniqueID(sys);
375     fSums->Add(fSysString);
376
377     fSums->Add(fSchemeString);
378     fSums->Add(fTriggerString);
379
380     // Print();
381   }
382   
383   PostData(1, fSums);
384 }
385
386 //________________________________________________________________________
387 void 
388 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
389                                           const char* title, const char* ytitle)
390 {
391   // 
392   // Set histogram graphical options, etc. 
393   // 
394   // Parameters:
395   //    h       Histogram to modify
396   //    colour  Marker color 
397   //    marker  Marker style
398   //    title   Title of histogram
399   //    ytitle  Title on y-axis. 
400   //
401   h->SetTitle(title);
402   h->SetMarkerColor(colour);
403   h->SetMarkerStyle(marker);
404   h->SetMarkerSize(1);
405   h->SetFillStyle(0);
406   h->SetYTitle(ytitle);
407   h->GetXaxis()->SetTitleFont(132);
408   h->GetXaxis()->SetLabelFont(132);
409   h->GetXaxis()->SetNdivisions(10);
410   h->GetYaxis()->SetTitleFont(132);
411   h->GetYaxis()->SetLabelFont(132);
412   h->GetYaxis()->SetNdivisions(10);
413   h->GetYaxis()->SetDecimals();
414   h->SetStats(0);
415 }
416
417 //________________________________________________________________________
418 TH1D*
419 AliBasedNdetaTask::ProjectX(const TH2D* h, 
420                             const char* name,
421                             Int_t firstbin, 
422                             Int_t lastbin, 
423                             bool  useRoot,
424                             bool  corr,
425                             bool  error)
426 {
427   // 
428   // Project onto the X axis 
429   // 
430   // Parameters:
431   //    h         2D histogram 
432   //    name      New name 
433   //    firstbin  First bin to use 
434   //    lastbin   Last bin to use
435   //    error     Whether to calculate errors
436   // 
437   // Return:
438   //    Newly created histogram or null
439   //
440   if (!h) return 0;
441   if (useRoot) 
442     return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
443   
444   TAxis* xaxis = h->GetXaxis();
445   TAxis* yaxis = h->GetYaxis();
446   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(), 
447                           xaxis->GetXmin(), xaxis->GetXmax());
448   static_cast<const TAttLine*>(h)->Copy(*ret);
449   static_cast<const TAttFill*>(h)->Copy(*ret);
450   static_cast<const TAttMarker*>(h)->Copy(*ret);
451   ret->GetXaxis()->ImportAttributes(xaxis);
452
453   Int_t first = firstbin; 
454   Int_t last  = lastbin;
455   if (first < 0)                         first = 0;
456   else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
457   if      (last  < 0)                    last  = yaxis->GetNbins();
458   else if (last  >  yaxis->GetNbins()+1) last  = yaxis->GetNbins();
459   if (last-first < 0) { 
460     AliWarningGeneral("AliBasedNdetaTask", 
461                       Form("Nothing to project [%d,%d]", first, last));
462     return 0;
463     
464   }
465
466   // Loop over X bins 
467   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
468   Int_t ybins = (last-first+1);
469   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) { 
470     Double_t content = 0;
471     Double_t error2  = 0;
472     Int_t    nbins   = 0;
473     
474     
475     for (Int_t ybin = first; ybin <= last; ybin++) { 
476       Double_t c1 = h->GetCellContent(xbin, ybin);
477       Double_t e1 = h->GetCellError(xbin, ybin);
478
479       // Ignore empty bins 
480       if (c1 < 1e-12) continue;
481       if (e1 < 1e-12) {
482         if (error) continue; 
483         e1 = 1;
484       }
485
486       content    += c1;
487       error2     += e1*e1;
488       nbins++;
489     } // for (ybin)
490     if(content > 0 && nbins > 0) {
491       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
492 #if 0
493       AliWarningGeneral(ret->GetName(), 
494                         Form("factor @ %d is %d/%d -> %f", 
495                              xbin, ybins, nbins, factor));
496 #endif
497       if (error) { 
498         // Calculate weighted average
499         ret->SetBinContent(xbin, content * factor);
500         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
501       }
502       else 
503         ret->SetBinContent(xbin, factor * content);
504     }
505   } // for (xbin)
506   
507   return ret;
508 }
509   
510 //________________________________________________________________________
511 void 
512 AliBasedNdetaTask::Terminate(Option_t *) 
513 {
514   // 
515   // Called at end of event processing.. 
516   //
517   // This is called once in the master 
518   // 
519   // Parameters:
520   //    option Not used 
521   //
522   // Draw result to screen, or perform fitting, normalizations Called
523   // once at the end of the query
524
525   fSums = dynamic_cast<TList*> (GetOutputData(1));
526   if(!fSums) {
527     AliError("Could not retrieve TList fSums"); 
528     return; 
529   }
530   
531   fOutput = new TList;
532   fOutput->SetName(Form("%s_result", GetName()));
533   fOutput->SetOwner();
534
535   fSNNString     = static_cast<TNamed*>(fSums->FindObject("sNN"));
536   fSysString     = static_cast<TNamed*>(fSums->FindObject("sys"));
537   fCentAxis      = static_cast<TAxis*>(fSums->FindObject("centAxis"));
538   fSchemeString  = static_cast<TNamed*>(fSums->FindObject("scheme"));
539   fTriggerString = static_cast<TNamed*>(fSums->FindObject("trigger"));
540
541   if(fSysString && fSNNString && 
542      fSysString->GetUniqueID() == AliForwardUtil::kPP)
543     LoadNormalizationData(fSysString->GetUniqueID(),
544                           fSNNString->GetUniqueID());
545   // Print before we loop
546   Print();
547
548   // Loop over centrality bins 
549   TIter next(fListOfCentralities);
550   CentralityBin* bin = 0;
551   while ((bin = static_cast<CentralityBin*>(next()))) 
552     bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff,
553              fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges, 
554              fTriggerMask, GetColor(), GetMarker());
555
556   // Output collision energy string 
557   if (fSNNString) fOutput->Add(fSNNString->Clone());
558
559   // Output collision system string 
560   if (fSysString) fOutput->Add(fSysString->Clone());
561
562   // Output centrality axis 
563   if (fCentAxis) fOutput->Add(fCentAxis);
564
565   // Output trigger string 
566   if (fTriggerString) fOutput->Add(fTriggerString->Clone());
567   
568   // Normalization string 
569   if (fSchemeString) fOutput->Add(fSchemeString->Clone());
570
571   // Output vertex axis 
572   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
573   vtxAxis->SetName("vtxAxis");
574   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
575   fOutput->Add(vtxAxis);
576
577   // Output trigger efficiency and shape correction 
578   fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff)); 
579   if (fShapeCorr) fOutput->Add(fShapeCorr);
580
581   TNamed* options = new TNamed("options","");
582   TString str;
583   str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
584   str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
585   str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
586   options->SetTitle(str);
587   fOutput->Add(options);
588
589   PostData(2, fOutput);
590 }
591 //________________________________________________________________________
592 void
593 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
594 {
595   // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
596   TString type("pp");
597   TString snn("900");
598   if(energy == 7000) snn.Form("7000");
599   if(energy == 2750) snn.Form("2750"); 
600   
601   if(fShapeCorr && (fTriggerEff != 1)) {
602     AliInfo("Objects already set for normalization - no action taken"); 
603     return; 
604   }
605   
606   TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWG2/FORWARD/corrections/"
607                                 "Normalization/normalizationHists_%s_%s.root",
608                                 type.Data(),snn.Data()));
609   if(!fin) {
610     AliWarning(Form("no file for normalization of %d/%d", sys, energy));
611     return;
612   }
613
614   // Shape correction 
615   TString shapeCorName(Form("h%sNormalization", 
616                             fTriggerMask == AliAODForwardMult::kInel ? "Inel" :
617                             fTriggerMask == AliAODForwardMult::kNSD ? "NSD" :
618                             fTriggerMask == AliAODForwardMult::kInelGt0 ?
619                             "InelGt0" : "All"));
620   TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
621   if (shapeCor) SetShapeCorrection(shapeCor);
622
623   // Trigger efficiency
624   TString effName(Form("%sTriggerEff", 
625                        fTriggerMask == AliAODForwardMult::kInel ? "inel" :
626                        fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
627                        fTriggerMask == AliAODForwardMult::kInelGt0 ?
628                        "inelgt0" : "all"));
629   TParameter<float>* eff = static_cast<TParameter<float>*>(fin->Get(effName));
630   if (eff) SetTriggerEff(eff->GetVal());
631
632   // TEMPORARY FIX
633   // Rescale the shape correction by the trigger efficiency 
634   AliWarning(Form("Rescaling shape correction by trigger efficiency: "
635                   "1/E_X=1/%f", fTriggerEff));
636   fShapeCorr->Scale(1. / fTriggerEff);
637
638   // Print - out
639   if (shapeCor && eff) AliInfo("Loaded objects for normalization.");
640 }
641
642
643 //________________________________________________________________________
644 void 
645 AliBasedNdetaTask::Print(Option_t*) const
646 {
647   // 
648   // Print information 
649   // 
650   std::cout << this->ClassName() << ": " << this->GetName() << "\n"
651             << std::boolalpha 
652             << " Trigger:              " << (fTriggerString ? 
653                                              fTriggerString->GetTitle() :
654                                              "none") << "\n"
655             << " Vertex range:         [" << fVtxMin << ":" << fVtxMax << "]\n"
656             << " Rebin factor:         " << fRebin << "\n" 
657             << " Cut edges:            " << fCutEdges << "\n"
658             << " Symmertrice:          " << fSymmetrice << "\n"
659             << " Use TH2::ProjectionX: " << fUseROOTProj << "\n"
660             << " Correct for empty:    " << fCorrEmpty << "\n"
661             << " Normalization scheme: " << (fSchemeString ? 
662                                              fSchemeString->GetTitle() : 
663                                              "none") <<"\n"
664             << " Trigger efficiency:   " << fTriggerEff << "\n" 
665             << " Shape correction:     " << (fShapeCorr ? 
666                                            fShapeCorr->GetName() : 
667                                            "none") << "\n"
668             << " sqrt(s_NN):           " << (fSNNString ? 
669                                              fSNNString->GetTitle() : 
670                                            "unknown") << "\n"
671             << " Collision system:     " << (fSysString ? 
672                                              fSysString->GetTitle() : 
673                                            "unknown") << "\n"
674             << " Centrality bins:      " << (fCentAxis ? "" : "none");
675   if (fCentAxis) { 
676     Int_t           nBins = fCentAxis->GetNbins();
677     const Double_t* bins  = fCentAxis->GetXbins()->GetArray();
678     for (Int_t i = 0; i <= nBins; i++) 
679       std::cout << (i==0 ? "" : "-") << bins[i];
680   }
681   std::cout << std::noboolalpha << std::endl;
682   
683 }
684
685 //________________________________________________________________________
686 TH1D*
687 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
688 {
689   // 
690   // Make a copy of the input histogram and rebin that histogram
691   // 
692   // Parameters:
693   //    h  Histogram to rebin
694   // 
695   // Return:
696   //    New (rebinned) histogram
697   //
698   if (rebin <= 1) return 0;
699
700   Int_t nBins = h->GetNbinsX();
701   if(nBins % rebin != 0) {
702     AliWarningGeneral("AliBasedNdetaTask", 
703                       Form("Rebin factor %d is not a devisor of current number "
704                            "of bins %d in the histogram %s", 
705                            rebin, nBins, h->GetName()));
706     return 0;
707   }
708     
709   // Make a copy 
710   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d", 
711                                                h->GetName(), rebin)));
712   tmp->Rebin(rebin);
713   tmp->SetDirectory(0);
714
715   // The new number of bins 
716   Int_t nBinsNew = nBins / rebin;
717   for(Int_t i = 1;i<= nBinsNew; i++) {
718     Double_t content = 0;
719     Double_t sumw    = 0;
720     Double_t wsum    = 0;
721     Int_t    nbins   = 0;
722     for(Int_t j = 1; j<=rebin;j++) {
723       Int_t    bin = (i-1)*rebin + j;
724       Double_t c   =  h->GetBinContent(bin);
725       if (c <= 0) continue;
726       
727       if (cutEdges) {
728         if (h->GetBinContent(bin+1)<=0 || 
729             h->GetBinContent(bin-1)<=0) {
730 #if 0
731           AliWarningGeneral("AliBasedNdetaTask", 
732                             Form("removing bin %d=%f of %s (%d=%f,%d=%f)", 
733                                  bin, c, h->GetName(), 
734                                  bin+1, h->GetBinContent(bin+1), 
735                                  bin-1, h->GetBinContent(bin-1)));
736 #endif
737           continue;
738         }       
739       }
740       Double_t e =  h->GetBinError(bin);
741       Double_t w =  1 / (e*e); // 1/c/c
742       content    += c;
743       sumw       += w;
744       wsum       += w * c;
745       nbins++;
746     }
747       
748     if(content > 0 && nbins > 0) {
749       tmp->SetBinContent(i, wsum / sumw);
750       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
751     }
752   }
753   
754   return tmp;
755 }
756
757 //__________________________________________________________________
758 TH1* 
759 AliBasedNdetaTask::Symmetrice(const TH1* h)
760 {
761   // 
762   // Make an extension of @a h to make it symmetric about 0 
763   // 
764   // Parameters:
765   //    h Histogram to symmertrice 
766   // 
767   // Return:
768   //    Symmetric extension of @a h 
769   //
770   Int_t nBins = h->GetNbinsX();
771   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
772   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
773   s->Reset();
774   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
775   s->SetMarkerColor(h->GetMarkerColor());
776   s->SetMarkerSize(h->GetMarkerSize());
777   s->SetMarkerStyle(h->GetMarkerStyle()+4);
778   s->SetFillColor(h->GetFillColor());
779   s->SetFillStyle(h->GetFillStyle());
780   s->SetDirectory(0);
781
782   // Find the first and last bin with data 
783   Int_t first = nBins+1;
784   Int_t last  = 0;
785   for (Int_t i = 1; i <= nBins; i++) { 
786     if (h->GetBinContent(i) <= 0) continue; 
787     first = TMath::Min(first, i);
788     last  = TMath::Max(last,  i);
789   }
790     
791   Double_t xfirst = h->GetBinCenter(first-1);
792   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
793   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
794   for (Int_t i = f1, j=l2; i <= last; i++,j--) { 
795     s->SetBinContent(j, h->GetBinContent(i));
796     s->SetBinError(j, h->GetBinError(i));
797   }
798   // Fill in overlap bin 
799   s->SetBinContent(l2+1, h->GetBinContent(first));
800   s->SetBinError(l2+1, h->GetBinError(first));
801   return s;
802 }
803
804 //====================================================================
805 AliBasedNdetaTask::CentralityBin::CentralityBin()
806   : TNamed("", ""), 
807     fSums(0), 
808     fOutput(0),
809     fSum(0), 
810     fSumMC(0), 
811     fTriggers(0), 
812     fLow(0), 
813     fHigh(0)
814 {
815   // 
816   // Constructor 
817   //
818 }
819 //____________________________________________________________________
820 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name, 
821                                                 Short_t low, Short_t high)
822   : TNamed(name, ""), 
823     fSums(0), 
824     fOutput(0),
825     fSum(0), 
826     fSumMC(0), 
827     fTriggers(0),
828     fLow(low), 
829     fHigh(high)
830 {
831   // 
832   // Constructor 
833   // 
834   // Parameters:
835   //    name Name used for histograms (e.g., Forward)
836   //    low  Lower centrality cut in percent 
837   //    high Upper centrality cut in percent 
838   //
839   if (low <= 0 && high <= 0) { 
840     fLow  = 0; 
841     fHigh = 0;
842     SetTitle("All centralities");
843   }
844   else {
845     fLow  = low;
846     fHigh = high;
847     SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
848   }
849 }
850 //____________________________________________________________________
851 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
852   : TNamed(o), 
853     fSums(o.fSums), 
854     fOutput(o.fOutput),
855     fSum(o.fSum), 
856     fSumMC(o.fSumMC), 
857     fTriggers(o.fTriggers), 
858     fLow(o.fLow), 
859     fHigh(o.fHigh)
860 {
861   // 
862   // Copy constructor 
863   // 
864   // Parameters:
865   //    other Object to copy from 
866   //
867 }
868 //____________________________________________________________________
869 AliBasedNdetaTask::CentralityBin::~CentralityBin()
870 {
871   // 
872   // Destructor 
873   //
874   if (fSums) fSums->Delete();
875   if (fOutput) fOutput->Delete();
876 }
877
878 //____________________________________________________________________
879 AliBasedNdetaTask::CentralityBin&
880 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
881 {
882   // 
883   // Assignment operator 
884   // 
885   // Parameters:
886   //    other Object to assign from 
887   // 
888   // Return:
889   //    Reference to this 
890   //
891   SetName(o.GetName());
892   SetTitle(o.GetTitle());
893   fSums      = o.fSums;
894   fOutput    = o.fOutput;
895   fSum       = o.fSum;
896   fSumMC     = o.fSumMC;
897   fTriggers  = o.fTriggers;
898   fLow       = o.fLow;
899   fHigh      = o.fHigh;
900
901   return *this;
902 }
903 //____________________________________________________________________
904 const char* 
905 AliBasedNdetaTask::CentralityBin::GetListName() const
906 {
907   // 
908   // Get the list name 
909   // 
910   // Return:
911   //    List Name 
912   //
913   if (IsAllBin()) return "all"; 
914   return Form("cent%03d_%03d", fLow, fHigh);
915 }
916 //____________________________________________________________________
917 void
918 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
919 {
920   // 
921   // Create output objects 
922   // 
923   // Parameters:
924   //    dir   Parent list
925   //
926   fSums = new TList;
927   fSums->SetName(GetListName());
928   fSums->SetOwner();
929   dir->Add(fSums);
930
931   fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers");
932   fTriggers->SetDirectory(0);
933   fSums->Add(fTriggers);
934 }
935 //____________________________________________________________________
936 void
937 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
938 {
939   // 
940   // Create sum histogram 
941   // 
942   // Parameters:
943   //    data  Data histogram to clone 
944   //    mc    (optional) MC histogram to clone 
945   //
946   fSum = static_cast<TH2D*>(data->Clone(GetName()));
947   fSum->SetDirectory(0);
948   fSum->Reset();
949   fSums->Add(fSum);
950   
951   // If no MC data is given, then do not create MC sum histogram 
952   if (!mc) return;
953
954   fSumMC = static_cast<TH2D*>(mc->Clone(Form("%sMC", GetName())));
955   fSumMC->SetTitle(Form("%s (MC)", fSumMC->GetTitle()));
956   fSumMC->SetDirectory(0);
957   fSumMC->Reset();
958   fSums->Add(fSumMC);
959 }
960
961 //____________________________________________________________________
962 Bool_t
963 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
964                                              Int_t triggerMask,
965                                              Double_t vzMin, Double_t vzMax)
966 {
967   // 
968   // Check the trigger, vertex, and centrality
969   // 
970   // Parameters:
971   //    aod Event input 
972   // 
973   // Return:
974   //    true if the event is to be used 
975   //
976   if (!forward) return false;
977
978   // We do not check for centrality here - it's already done 
979   return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, fTriggers);
980 }
981   
982   
983 //____________________________________________________________________
984 void
985 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
986                                                Int_t triggerMask,
987                                                Double_t vzMin, Double_t vzMax,
988                                                const TH2D* data, const TH2D* mc)
989 {
990   // 
991   // Process an event
992   // 
993   // Parameters:
994   //    forward     Forward data (for trigger, vertex, & centrality)
995   //    triggerMask Trigger mask 
996   //    vzMin       Minimum IP z coordinate
997   //    vzMax       Maximum IP z coordinate
998   //    data        Data histogram 
999   //    mc          MC histogram
1000   //
1001   if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
1002   if (!data) return;
1003   if (!fSum) CreateSums(data, mc);
1004   
1005   fSum->Add(data);
1006   if (mc) fSumMC->Add(mc);
1007 }
1008
1009 //________________________________________________________________________
1010 Double_t 
1011 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1012                                                 UShort_t    scheme,
1013                                                 Double_t    trigEff,
1014                                                 Double_t&   ntotal) const
1015 {
1016   // 
1017   // Calculate normalization 
1018   // 
1019   // Parameters: 
1020   //    t       Trigger histogram
1021   //    scheme  Normaliztion scheme 
1022   //    trigEff From MC
1023   //    ntotal  On return, contains the number of events. 
1024   //
1025   ntotal               = 0;
1026   Double_t nAll        = t.GetBinContent(AliAODForwardMult::kBinAll);
1027   Double_t nB          = t.GetBinContent(AliAODForwardMult::kBinB);
1028   Double_t nA          = t.GetBinContent(AliAODForwardMult::kBinA);
1029   Double_t nC          = t.GetBinContent(AliAODForwardMult::kBinC);
1030   Double_t nE          = t.GetBinContent(AliAODForwardMult::kBinE);
1031   Double_t nOffline    = t.GetBinContent(AliAODForwardMult::kBinOffline);
1032   Double_t nTriggered  = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1033   Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1034   Double_t nAccepted   = t.GetBinContent(AliAODForwardMult::kAccepted);
1035   
1036   if (nTriggered <= 0.1) { 
1037     AliError("Number of triggered events <= 0");
1038     return -1;
1039   }
1040   if (nWithVertex <= 0.1) { 
1041     AliError("Number of events with vertex <= 0");
1042     return -1;
1043   }
1044   ntotal          = nAccepted;
1045   Double_t vtxEff = nWithVertex / nTriggered;
1046   Double_t scaler = 1;
1047   Double_t beta   = nA + nC - 2*nE;
1048
1049   if (scheme & kEventLevel) {
1050     ntotal = nAccepted / vtxEff;
1051     scaler = vtxEff;
1052     AliInfo(Form("Calculating event normalisation as\n"
1053                  " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1054                  Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex), 
1055                  ntotal, scaler));
1056             
1057     if (scheme & kBackground) {
1058       //          1            E_V             E_V
1059       //   s = --------- = ------------- = ------------ 
1060       //        1 - beta   1 - beta E_V    1 - beta N_V 
1061       //       ---  ----       --------        ---- ---
1062       //       E_V  N_V          N_V           N_V  N_T 
1063       // 
1064       //          E_V
1065       //     = ------------
1066       //        1 - beta 
1067       //            ----
1068       //             N_T 
1069       // 
1070       ntotal -= nAccepted * beta / nWithVertex;
1071       // This one is direct and correct. 
1072       // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1073       // A simpler expresion
1074       scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1075       AliInfo(Form("Calculating event normalisation as\n"
1076                    " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1077                    " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1078                    Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1079                    nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta), 
1080                    Int_t(nWithVertex), ntotal, scaler));
1081     }
1082   }
1083   if (scheme & kTriggerEfficiency) {
1084     ntotal /= trigEff;
1085     scaler *= trigEff;
1086     AliInfo(Form("Correcting for trigger efficiency:\n"
1087                  " N = 1 / E_X * N = 1 / %f * %d = %f (%f)", 
1088                  trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1089   }
1090
1091   AliInfo(Form("\n"
1092                " Total of        %9d events for %s\n"
1093                "  of these       %9d have an offline trigger\n"
1094                "  of these N_T = %9d has the selected trigger\n" 
1095                "  of these N_V = %9d has a vertex\n" 
1096                "  of these N_A = %9d were in the selected range\n"
1097                "  Triggers by hardware type:\n"
1098                "    N_b   =  %9d\n"
1099                "    N_ac  =  %9d (%d+%d)\n"
1100                "    N_e   =  %9d\n"
1101                "  Vertex efficiency:          %f\n"
1102                "  Trigger efficiency:         %f\n"
1103                "  Total number of events: N = %f\n"
1104                "  Scaler (N_A/N):             %f",
1105                Int_t(nAll), GetTitle(), Int_t(nOffline), 
1106                Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1107                Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE), 
1108                vtxEff, trigEff, ntotal, scaler));
1109   return scaler;
1110 }
1111
1112 //________________________________________________________________________
1113 void 
1114 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,  
1115                                              const char* postfix, 
1116                                              bool        rootProj, 
1117                                              bool        corrEmpty,
1118                                              const TH1*  shapeCorr,
1119                                              Double_t    scaler,
1120                                              bool        symmetrice, 
1121                                              Int_t       rebin, 
1122                                              bool        cutEdges, 
1123                                              Int_t       color, 
1124                                              Int_t       marker)
1125 {
1126   // 
1127   // Generate the dN/deta result from input 
1128   // 
1129   // Parameters: 
1130   //     sum        Sum of 2D hists 
1131   //     postfix    Post fix on names
1132   //     rootProj   Whether to use ROOT TH2::ProjectionX
1133   //     corrEmpty  Correct for empty bins 
1134   //     shapeCorr  Shape correction to use 
1135   //     scaler     Event-level normalization scaler  
1136   //     symmetrice Whether to make symmetric extensions 
1137   //     rebin      Whether to rebin
1138   //     cutEdges   Whether to cut edges when rebinning 
1139   //
1140   TH2D* copy    = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s", 
1141                                                      GetName(), postfix)));
1142   TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), 0, 0, 
1143                            rootProj, corrEmpty, false);
1144   accNorm->SetDirectory(0);
1145
1146   // ---- Scale by shape correction ----------------------------------
1147   if (shapeCorr) copy->Divide(shapeCorr);
1148   else AliInfo("No shape correction specified, or disabled");
1149   
1150   // Normalize to the acceptance -
1151   // dndeta->Divide(accNorm);
1152   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { 
1153     for (Int_t j = 1; j <= copy->GetNbinsY(); j++) { 
1154       Double_t c = copy->GetBinContent(i, j);
1155       Double_t e = copy->GetBinError(i, j);
1156       Double_t a = accNorm->GetBinContent(i);
1157       copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
1158       copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
1159     }
1160   }
1161   // --- Event-level normalization -----------------------------------
1162   copy->Scale(scaler);
1163
1164   // --- Project on X axis -------------------------------------------
1165   TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix), 
1166                           1, sum->GetNbinsY(), rootProj, corrEmpty);
1167   dndeta->SetDirectory(0);
1168   // Event-level normalization 
1169   dndeta->Scale(1., "width");
1170   copy->Scale(1., "width");
1171
1172   // --- Set some histogram attributes -------------------------------
1173   SetHistogramAttributes(dndeta,  color, marker, Form("ALICE %s", GetName()));
1174   SetHistogramAttributes(accNorm, color, marker, Form("ALICE %s normalisation", 
1175                                                       GetName())); 
1176
1177   // --- Make symmetric extensions and rebinnings --------------------
1178   if (symmetrice)   fOutput->Add(Symmetrice(dndeta));
1179   fOutput->Add(dndeta);
1180   fOutput->Add(accNorm);
1181   fOutput->Add(copy);
1182   fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1183   if (symmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1184 }  
1185
1186 //________________________________________________________________________
1187 void 
1188 AliBasedNdetaTask::CentralityBin::End(TList*      sums, 
1189                                       TList*      results,
1190                                       UShort_t    scheme,
1191                                       const TH1*  shapeCorr, 
1192                                       Double_t    trigEff,
1193                                       Bool_t      symmetrice,
1194                                       Int_t       rebin, 
1195                                       Bool_t      rootProj,
1196                                       Bool_t      corrEmpty, 
1197                                       Bool_t      cutEdges, 
1198                                       Int_t       triggerMask,
1199                                       Int_t       color, 
1200                                       Int_t       marker) 
1201 {
1202   // 
1203   // End of processing 
1204   // 
1205   // Parameters:
1206   //    sums        List of sums
1207   //    results     Output list of results
1208   //    shapeCorr   Shape correction or nil
1209   //    trigEff     Trigger efficiency 
1210   //    symmetrice  Whether to symmetrice the results
1211   //    rebin       Whether to rebin the results
1212   //    corrEmpty   Whether to correct for empty bins
1213   //    cutEdges    Whether to cut edges when rebinning
1214   //    triggerMask Trigger mask 
1215   //
1216
1217   fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
1218   if(!fSums) {
1219     AliError("Could not retrieve TList fSums"); 
1220     return; 
1221   }
1222   
1223   fOutput = new TList;
1224   fOutput->SetName(GetListName());
1225   fOutput->SetOwner();
1226   results->Add(fOutput);
1227
1228   fSum      = static_cast<TH2D*>(fSums->FindObject(GetName()));
1229   fSumMC    = static_cast<TH2D*>(fSums->FindObject(Form("%sMC", GetName())));
1230   fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
1231
1232   if (!fTriggers) { 
1233     AliError("Couldn't find histogram 'triggers' in list");
1234     return;
1235   }
1236   if (!fSum) { 
1237     AliError(Form("Couldn't find histogram '%s' in list", GetName()));
1238     return;
1239   }
1240
1241   // --- Get normalization scaler ------------------------------------
1242   Double_t ntotal = 0;
1243   Double_t epsilonT = trigEff;
1244   // TEMPORARY FIX
1245   if (triggerMask == AliAODForwardMult::kNSD) {
1246     // This is a local change 
1247     epsilonT = 0.92; 
1248     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
1249   }
1250   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal);
1251   if (scaler < 0) { 
1252     AliError("Failed to calculate normalization - bailing out");
1253     return;
1254   }
1255   fOutput->Add(fTriggers->Clone());
1256
1257   // --- Make result and store ---------------------------------------
1258   MakeResult(fSum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
1259              scaler, symmetrice, rebin, cutEdges, color, marker);
1260
1261   // --- Process result from TrackRefs -------------------------------
1262   if (fSumMC) { 
1263     MakeResult(fSumMC, "MC", rootProj, corrEmpty, 
1264                (scheme & kShape) ? shapeCorr : 0,
1265                scaler, symmetrice, rebin, cutEdges, color+2, marker);
1266   }
1267
1268   // Temporary stuff 
1269   if (!IsAllBin()) return;
1270
1271   TFile* forward = TFile::Open("forward.root", "READ");
1272   if (!forward)  { 
1273     AliWarning(Form("No forward.root file found"));
1274     return;
1275   }
1276
1277   TH1D* shapeCorrProj = 0;
1278   if (shapeCorr) {
1279     shapeCorrProj = static_cast<const TH2D*>(shapeCorr)->ProjectionX();
1280     shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY());
1281     shapeCorrProj->SetDirectory(0);
1282     fOutput->Add(shapeCorrProj);
1283   }
1284
1285   TList* official = static_cast<TList*>(forward->Get("official"));
1286   if (official) { 
1287     TH1F*  histEta  = static_cast<TH1F*>(official->FindObject("fHistEta"));
1288     if (histEta)  {
1289       TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), 
1290                             histEta->GetXaxis()->GetXmin(), 
1291                             histEta->GetXaxis()->GetXmax());
1292       for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
1293         oEta->SetBinContent(i, histEta->GetBinContent(i));
1294         oEta->SetBinError(i, histEta->GetBinError(i));
1295       }
1296       if (shapeCorrProj) oEta->Divide(shapeCorrProj);
1297       oEta->Scale(1./ntotal, "width");
1298       oEta->SetDirectory(0);
1299       oEta->SetMarkerStyle(marker+4);
1300       oEta->SetMarkerColor(color+5);
1301       fOutput->Add(oEta);
1302       fOutput->Add(Rebin(oEta, rebin, false));
1303     }
1304     else 
1305       AliWarning(Form("Couldn't find histogram fHistEta in list %s", 
1306                       official->GetName()));
1307   }
1308   else 
1309     AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName()));
1310
1311   TList* tracks = static_cast<TList*>(forward->Get("tracks"));
1312   if (tracks) { 
1313     TH1F*  histEta  = static_cast<TH1F*>(tracks->FindObject("fHistEta"));
1314     if (histEta)  {
1315       TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), 
1316                             histEta->GetXaxis()->GetXmin(), 
1317                             histEta->GetXaxis()->GetXmax());
1318       for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
1319         oEta->SetBinContent(i, histEta->GetBinContent(i));
1320         oEta->SetBinError(i, histEta->GetBinError(i));
1321       }
1322       if (shapeCorrProj) oEta->Divide(shapeCorrProj);
1323       oEta->Scale(1./ntotal, "width");
1324       oEta->SetDirectory(0);
1325       oEta->SetMarkerStyle(marker);
1326       oEta->SetMarkerColor(color+5);
1327       fOutput->Add(oEta);
1328       fOutput->Add(Rebin(oEta, rebin, false));
1329     }
1330     else 
1331       AliWarning(Form("Couldn't find histogram fHistEta in list %s", 
1332                       tracks->GetName()));
1333   }
1334   else 
1335     AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName()));
1336
1337   forward->Close();
1338 }
1339
1340 //
1341 // EOF
1342 //