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