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