AliPhysics  32b88a8 (32b88a8)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliBasedNdetaTask.cxx
Go to the documentation of this file.
1 //====================================================================
2 #include "AliBasedNdetaTask.h"
3 #include <TMath.h>
4 #include <TH2D.h>
5 #include <TH1D.h>
6 #include <TF1.h>
7 #include <THStack.h>
8 #include <TList.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 "AliAODMultEventClass.h"
16 #include <TFile.h>
17 #include <TStyle.h>
18 #include <TROOT.h>
19 #include <TParameter.h>
20 #include <TColor.h>
21 
22 //____________________________________________________________________
24  : AliBaseAODTask(),
25  fCorrEmpty(true),
26  fUseROOTProj(false),
27  fTriggerEff(1),
28  fTriggerEff0(1),
29  fListOfCentralities(0),
30  fNormalizationScheme(kFull),
31  fFinalMCCorrFile(""),
32  fSatelliteVertices(0),
33  fEmpiricalCorrection(0),
34  fMeanVsC(0),
35  fSeenCent(0),
36  fTakenCent(0),
37  fCentMethod("default"),
38  fAnaUtil(),
39  fUseUtilPileup(false),
40  fIpzReweight(0),
41  fCacheCent(-10),
42  fCacheQual(0xFFFF)
43 {
44  //
45  // Constructor
46  //
47  DGUARD(fDebug,3,"Default CTOR of AliBasedNdetaTask");
48 }
49 
50 //____________________________________________________________________
52  : AliBaseAODTask(Form("%sdNdeta", name),"AliBasedNdetaTask"),
53  fCorrEmpty(true),
54  fUseROOTProj(false),
55  fTriggerEff(1),
56  fTriggerEff0(1),
57  fListOfCentralities(0),
58  fNormalizationScheme(kFull),
59  fFinalMCCorrFile(""),
60  fSatelliteVertices(0),
61  fEmpiricalCorrection(0),
62  fMeanVsC(0),
63  fSeenCent(0),
64  fTakenCent(0),
65  fCentMethod("default"),
66  fAnaUtil(),
67  fUseUtilPileup(false),
68  fIpzReweight(0),
69  fCacheCent(-10),
70  fCacheQual(0xFFFF)
71 {
72  //
73  // Constructor
74  //
75  DGUARD(fDebug, 3,"Named CTOR of AliBasedNdetaTask: %s", name);
76 
78  fMinIpZ = -10;
79  fMaxIpZ = +10;
81  fListOfCentralities->SetName("centralityBins");
82  // fListOfCentralities->SetTitle("Centrality bins");
83 
84  // Set the normalisation scheme
86 }
87 
88 
89 //____________________________________________________________________
91 {
92  //
93  // Destructor
94  //
95  DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask");
96 }
97 
98 //________________________________________________________________________
99 void
101 {
102  AliAnalysisTaskSE::SetDebugLevel(lvl);
103  for (Int_t i = 0; i < fListOfCentralities->GetEntries(); i++) {
104  CentralityBin* bin =
105  static_cast<CentralityBin*>(fListOfCentralities->At(i));
106  bin->SetDebugLevel(lvl);
107  }
108 }
109 
110 //________________________________________________________________________
111 void
113 {
114  //
115  // Add a centrality bin
116  //
117  // Parameters:
118  // low Low cut
119  // high High cut
120  //
121  DGUARD(fDebug,3,"Add a centrality bin [%6.2f,%6.2f] @ %d", low, high, at);
122  CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
123  if (!bin) {
124  Error("AddCentralityBin",
125  "Failed to create centrality bin for %s [%6.2f,%6.2f] @ %d",
126  GetName(), low, high, at);
127  return;
128  }
130  bin->SetDebugLevel(fDebug);
131  fListOfCentralities->AddAtAndExpand(bin, at);
132 }
133 
134 //________________________________________________________________________
137  Float_t low, Float_t high) const
138 {
139  //
140  // Make a centrality bin
141  //
142  // Parameters:
143  // name Name used for histograms
144  // low Low cut in percent
145  // high High cut in percent
146  //
147  // Return:
148  // A newly created centrality bin
149  //
150  DGUARD(fDebug,3,"Make a centrality bin %s [%6.2f,%6.2f]", name, low, high);
151  return new CentralityBin(name, low, high);
152 }
153 
154 #define TESTAPPEND(SCHEME,BIT,STRING) \
155  do { if (!(SCHEME & BIT)) break; \
156  if (!s.IsNull()) s.Append(","); s.Append(STRING); } while(false)
157 
158 //________________________________________________________________________
159 const Char_t*
161 {
162  // Create a string from normalization scheme bits
163  static TString s;
164  s = "";
165 
166  if (scheme == kNone)
167  return s.Data();
168  if (scheme == kFull) {
169  s = "FULL";
170  return s.Data();
171  }
172  TESTAPPEND(scheme, kEventLevel, "EVENT");
173  TESTAPPEND(scheme, kBackground, "BACKGROUND");
174  TESTAPPEND(scheme, kTriggerEfficiency, "TRIGGER");
175  TESTAPPEND(scheme, kZeroBin, "ZEROBIN");
176 
177  return s.Data();
178 }
179 //________________________________________________________________________
180 UShort_t
182 {
183  UShort_t scheme = 0;
184  TString twhat(what);
185  twhat.ToUpper();
186  if (twhat.EqualTo("DEFAULT"))
188  TObjString* opt;
189  TObjArray* token = twhat.Tokenize(" ,|");
190  TIter next(token);
191  while ((opt = static_cast<TObjString*>(next()))) {
192  TString s(opt->GetString());
193  if (s.IsNull()) continue;
194  Bool_t add = true;
195  switch (s[0]) {
196  case '-': add = false; // Fall through
197  case '+': s.Remove(0,1); // Remove character
198  }
199  UShort_t bit = 0;
200  if (s.EqualTo("SHAPE")) {
201  AliWarningGeneral("AliBasedNdetaTask",
202  Form("SHAPE correction no longer supported (%s)",
203  what));
204  continue;
205  }
206  if (s.CompareTo("EVENT") == 0) bit = kEventLevel;
207  else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
208  else if (s.CompareTo("TRIGGER") == 0) bit = kTriggerEfficiency;
209  else if (s.CompareTo("FULL") == 0) bit = kFull;
210  else if (s.CompareTo("NONE") == 0) bit = kNone;
211  else if (s.CompareTo("ZEROBIN") == 0) bit = kZeroBin;
212  else
213  ::Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
214  if (add) scheme |= bit;
215  else scheme ^= bit;
216  }
217  delete token;
218  return scheme;
219 }
220 //________________________________________________________________________
221 void
223 {
224  //
225  // Set normalisation scheme
226  //
227  DGUARD(fDebug,3,"Set the normalization scheme: %s", what);
229 }
230 //________________________________________________________________________
231 void
233 {
234  DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme);
235  fNormalizationScheme = scheme;
236 }
237 //____________________________________________________________________
238 Bool_t
240 {
241  Int_t id = GetCentMethodID(method);
242  if (id < -1) {
243  AliErrorF("Unknown centrality estimator: %s", method.Data());
244  fCentMethod = "";
245  return false;
246  }
247  if (id < 0) {
248  // Do not use any estimator
249  AliInfoF("No centrality estimator: \"%s\"", method.Data());
250  fCentMethod = "";
251  return false;
252  }
253 
254  TString meth = GetCentMethod(id);
255  if (fName.Contains("Forward", TString::kIgnoreCase) && meth.Contains("FMD"))
256  AliWarningF("Centrality estimator %s used by %s - beware of auto-corr",
257  meth.Data(), fName.Data());
258  else if (fName.Contains("Central", TString::kIgnoreCase) &&
259  (meth.Contains("CL0") || meth.Contains("TKL")))
260  AliWarningF("Centrality estimator %s used by %s - beware of auto-corr",
261  meth.Data(), fName.Data());
262 
263  fCentMethod = meth;
264  AliInfoF("Centrality estimator set to %s", fCentMethod.Data());
265  return true;
266 }
267 
268 //________________________________________________________________________
269 Int_t
271 {
272  Int_t ret = -2;
273  TString m(meth);
274  m.ToUpper();
275  if (m.EqualTo("NONE") || m.EqualTo("NO") || m.EqualTo("FALSE"))
276  ret = kCentNone;
277  else if (m.IsNull()) ret = kCentDefault;
278  else if (m.BeginsWith("DEFAULT")) ret = kCentDefault;
279  else if (m.BeginsWith("ZEMVSZDC")) ret = kCentZEMvsZDC;
280  else if (m.BeginsWith("TKLVSV0M")) ret = kCentTklvsV0M;
281  else if (m.BeginsWith("V0MVSFMD")) ret = kCentV0MvsFMD;
282  else if (m.BeginsWith("NPA")) ret = kCentNPA;
283  else if (m.BeginsWith("ZPC")) ret = kCentZPC;
284  else if (m.BeginsWith("ZPA")) ret = kCentZPA;
285  else if (m.BeginsWith("ZNC")) ret = kCentZNC;
286  else if (m.BeginsWith("ZNA")) ret = kCentZNA;
287  else if (m.BeginsWith("CND")) ret = kCentCND;
288  else if (m.BeginsWith("CL1")) ret = kCentCL1;
289  else if (m.BeginsWith("CL0")) ret = kCentCL0;
290  else if (m.BeginsWith("TKL")) ret = kCentTkl;
291  else if (m.BeginsWith("TRK")) ret = kCentTrk;
292  else if (m.BeginsWith("FMD")) ret = kCentFMD;
293  else if (m.BeginsWith("V0C")) ret = kCentV0C;
294  else if (m.BeginsWith("V0A123")) ret = kCentV0A123;
295  else if (m.BeginsWith("V0A")) ret = kCentV0A;
296  else if (m.BeginsWith("V0M")) ret = kCentV0M;
297  else if (m.BeginsWith("MULTV0A")) ret = kMultV0A;
298  else if (m.BeginsWith("MULTV0M")) ret = kMultV0M;
299  else if (m.BeginsWith("MULTV0C")) ret = kMultV0C;
300  else if (m.BeginsWith("MULT")) ret = kMult;
301  if (m.Contains("TRUE")) ret |= kCentTrue;
302  if (m.Contains("EQ")) ret |= kCentEq;
303 
304  return ret;
305 }
306 //________________________________________________________________________
307 const char*
309 {
310  static TString ret("");
311  UShort_t base = (id & 0xFF);
312  switch (base) {
313  case kCentNone: ret = "none"; break;
314  case kCentDefault: ret = ""; break;
315  case kCentV0M: ret = "V0M"; break;
316  case kCentV0A: ret = "V0A"; break;
317  case kCentV0A123: ret = "V0A123"; break;
318  case kCentV0C: ret = "V0C"; break;
319  case kCentFMD: ret = "FMD"; break;
320  case kCentTrk: ret = "TRK"; break;
321  case kCentTkl: ret = "TKL"; break;
322  case kCentCL0: ret = "CL0"; break;
323  case kCentCL1: ret = "CL1"; break;
324  case kCentCND: ret = "CND"; break;
325  case kCentZNA: ret = "ZNA"; break;
326  case kCentZNC: ret = "ZNC"; break;
327  case kCentZPA: ret = "ZPA"; break;
328  case kCentZPC: ret = "ZPC"; break;
329  case kCentNPA: ret = "NPA"; break;
330  case kCentV0MvsFMD: ret = "V0MvsFMD"; break;
331  case kCentTklvsV0M: ret = "TKLvsV0M"; break;
332  case kCentZEMvsZDC: ret = "ZEMvsZDC"; break;
333  case kMult: ret = "MULT"; break;
334  case kMultV0A: ret = "MULTV0A"; break;
335  case kMultV0M: ret = "MULTV0M"; break;
336  case kMultV0C: ret = "MULTV0C"; break;
337  default: ret = ""; break;
338  }
339  Bool_t tru = id & kCentTrue;
340  Bool_t eq = id & kCentEq;
341  if (eq) {
342  if (!tru) ret.Append("eq");
343  else ret.Append("Eq");
344  }
345  if (tru) ret.Append("true");
346 
347  return ret.Data();
348 }
349 
350 
351 //________________________________________________________________________
352 void
354 {
355  if (fListOfCentralities->GetEntries() > 0) return;
356 
357  // Automatically add 'all' centrality bin if nothing has been defined.
358  AddCentralityBin(0, 0, 0);
359 
360  // Check if the centrality method was set to none, and in that case
361  // remove the centrality axis.
362  if (fCentMethod.EqualTo("none", TString::kIgnoreCase)) {
363  fCentAxis.Set(0,0,0);
364  TH1* h = static_cast<TH1*>(fSums->FindObject(fCentAxis.GetName()));
365  if (h) {
366  Info("InitializeCentBins",
367  "No centrality, forcing centrality axis to null");
368  h->GetXaxis()->Set(0,0,0);
369  h->Rebuild();
370  }
371  }
372  if (HasCentrality()) {
373  const TArrayD* bins = fCentAxis.GetXbins();
374  Int_t nbin = fCentAxis.GetNbins();
375  for (Int_t i = 0; i < nbin; i++)
376  AddCentralityBin(i+1, (*bins)[i], (*bins)[i+1]);
377  AddCentralityBin(nbin+1, 0, 101);
378  }
379 }
380 
381 //________________________________________________________________________
382 Bool_t
384 {
385  //
386  // Create output objects.
387  //
388  // This is called once per slave process
389  //
390  DGUARD(fDebug,1,"Create user ouput object");
391 
392  // Modify axis of centrality histograms
393  fCent->SetXTitle(Form("Centrality (%s) [%%]", fCentMethod.Data()));
394  fAccCent->SetXTitle(Form("Centrality (%s) [%%]", fCentMethod.Data()));
395 
396  fSums->Add(AliForwardUtil::MakeParameter("empirical",
397  fEmpiricalCorrection != 0));
399  fSums->Add(AliForwardUtil::MakeParameter("centEstimator",
401  if (fIpzReweight) fSums->Add(fIpzReweight->Clone("ipZw"));
402  // fSums->Add(new TNamed("centEstimator", fCentMethod.Data()));
403 
404  // Make our centrality bins
406 
407  // Loop over centrality bins
408  TIter next(fListOfCentralities);
409  CentralityBin* bin = 0;
410  while ((bin = static_cast<CentralityBin*>(next())))
412 
413 
414  if (HasCentrality()) {
415  if (fCentAxis.GetXbins()->GetArray()) {
416  Int_t nbin = fCentAxis.GetNbins();
417  TArrayD a(nbin+1, fCentAxis.GetXbins()->GetArray());
418  a.Set(nbin+2);
419  a[nbin+1] = fCentAxis.GetXmax()+fCentAxis.GetBinWidth(nbin+1);
420 
421  fSeenCent = new TH1D("centSeen", "Centralities seen",
422  a.GetSize()-1, a.GetArray());
423  fMeanVsC=new TH2D("sumVsC",
424  "Integral vs centrality",
425  a.GetSize()-1, a.GetArray(),
426  400, 0, 20000);
427  }
428  else {
429  fSeenCent = new TH1D("centSeen", "Centralities seen",
430  fCentAxis.GetNbins()+1,
431  fCentAxis.GetXmin(),
432  fCentAxis.GetXmax()+fCentAxis.GetBinWidth(1));
433  fMeanVsC=new TH2D("sumVsC",
434  "Integral vs centrality",
435  fCentAxis.GetNbins()+1,
436  fCentAxis.GetXmin(),
437  fCentAxis.GetXmax()+fCentAxis.GetBinWidth(1),
438  400, 0, 20000);
439  }
440  }
441  else {
442  fSeenCent = new TH1D("centSeen", "Null",1, 0, 100);
443  fMeanVsC = new TH2D("sumVsC", "Null",1,0,20000,1,0,100);
444  }
445  fSeenCent->SetDirectory(0);
446  fSeenCent->SetXTitle(Form("Centrality (%s) [%%]",fCentMethod.Data()));
447  fSeenCent->SetYTitle("Events");
448  fSeenCent->SetFillStyle(3004);
449  fSeenCent->SetFillColor(kYellow-2);
450  fMeanVsC->SetDirectory(0);
451  fMeanVsC->SetXTitle(fSeenCent->GetXaxis()->GetTitle());
452  fMeanVsC->SetYTitle("#int signal");
453 
454  fTakenCent = static_cast<TH1D*>(fSeenCent->Clone("centTaken"));
455  fTakenCent->SetTitle("Centralities taken by bins");
456  fTakenCent->SetFillColor(kCyan-2);
457  fTakenCent->SetFillStyle(3005);
458  fTakenCent->SetDirectory(0);
459 
460  fSums->Add(fSeenCent);
461  fSums->Add(fTakenCent);
462  fSums->Add(fMeanVsC);
463 
464  // fSums->ls();
465  return true;
466 }
467 
468 //____________________________________________________________________
469 Bool_t
471 {
472  if (!fCentMethod.BeginsWith("MULT"))
474  else
477 
478  // Here, we always return true, as the centrality bins will do
479  // their own checks on the events - this is needed for event
480  // normalization etc.
481  return true;
482 }
483 
484 //____________________________________________________________________
485 Double_t
487  AliAODForwardMult* forward,
488  Int_t& qual)
489 {
490  DGUARD(fDebug,2,"Getting centrality from event of object: %s",
491  fCentMethod.Data());
492  if (fCacheCent > -1) {
493  // In case we already got the centrality, don't do it again
494  DMSG(fDebug,1,"Returning cached value: %5.1f%%", fCent);
495  qual = fCacheQual;
496  return fCacheCent;
497  }
498  qual = 0;
499  Float_t cent = AliBaseAODTask::GetCentrality(event,forward,qual);
500  DMSG(fDebug,1,"Centrality stored in AOD forward: %5.1f%%", cent);
501  if (!fCentMethod.IsNull()) {
502  // Clear bad cent if already set
503  cent = AliForwardUtil::GetCentrality(event,fCentMethod,qual,(fDebug > 1));
504  if (cent < 0) cent = -.5; // Bad centrality
505  else if (TMath::Abs(cent-100) < 1.1) cent = 100; // Special centralities
506  DMSG(fDebug,1,"Centrality from mult: %5.1f%% (%d)", cent, qual);
507  }
508  fCacheQual = qual;
509  return fCacheCent = cent;
510 }
511 //____________________________________________________________________
512 Double_t
514  AliAODForwardMult* forward)
515 {
516  Int_t qual = 0;
517  Double_t cent = GetCentrality(event, forward, qual);
518  if (qual > 0) forward->SetTriggerBits(AliAODForwardMult::kCentNoCalib);
519  return cent;
520 }
521 
522 //____________________________________________________________________
523 Bool_t
525 {
526  //
527  // Process a single event
528  //
529  // Parameters:
530  // option Not used
531  //
532  // Main loop
533  DGUARD(fDebug,1,"Analyse the AOD event");
534  if (fUseUtilPileup && fAnaUtil.IsPileUpEvent(&aod)) return false;
535 
536  AliAODForwardMult* forward = GetForward(aod);
537  if (!forward) return false;
538 
539  // Fill centrality histogram
540 
541  Double_t vtx = forward->GetIpZ();
542  TH2D* data = GetHistogram(aod, false);
543  TH2D* dataMC = GetHistogram(aod, true);
544  if (!data) return false;
545 
546  CheckEventData(vtx, data, dataMC);
547 
548  if (!ApplyEmpiricalCorrection(forward,data))
549  return false;
550 
551  Double_t ipzW = 1;
552  if (fIpzReweight) {
553  ipzW = fIpzReweight->Eval(vtx);
554  DMSG(fDebug,5,"IPz=%f -> Weight %f", vtx, ipzW);
555  }
556 
557  Bool_t isZero = ((fNormalizationScheme & kZeroBin) &&
559  Bool_t taken = false;
560  // Loop over centrality bins
561  CentralityBin* allBin =
562  static_cast<CentralityBin*>(fListOfCentralities->At(0));
563  if (allBin->ProcessEvent(forward, fTriggerMask, isZero, fMinIpZ, fMaxIpZ,
564  data, dataMC, fFilterMask, ipzW)) taken = true;
565 
566  // Find this centrality bin
567  if (HasCentrality()) {
568  taken = false;
569  // After this call, the if the event isn't covered by the
570  // centrality calibrations, a flag will be set in the forward
571  // object.
572  Double_t cent = GetCentrality(aod, forward);
573  fMeanVsC->Fill(cent, data->Integral());
574 
575  // fSeenCent->Fill(cent);
576  DMSG(fDebug,1,"Got event centrality %f (%s)", cent,
579  "out-of-calib" : "in-calib");
580 
581  Int_t icent = fCentAxis.FindBin(cent);
582  if (icent == (fCentAxis.GetNbins()+1) &&
583  TMath::Abs(cent-fCentAxis.GetXmax()) < 1e-6)
584  // Make sure that upper edge is analysed
585  icent = fCentAxis.GetNbins();
586  CentralityBin* thisBin = 0;
587  if (icent >= 1 && icent <= fCentAxis.GetNbins())
588  thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
589  if (thisBin && thisBin->ProcessEvent(forward, fTriggerMask,
590  isZero, fMinIpZ, fMaxIpZ,
591  data, dataMC, fFilterMask, ipzW))
592  taken = true;
593  if (taken) fTakenCent->Fill(cent);
594 
595  Int_t nbins = fCentAxis.GetNbins();
596  CentralityBin* fullBin =
597  static_cast<CentralityBin*>(fListOfCentralities->At(nbins+1));
598  if (fullBin && fullBin->ProcessEvent(forward, fTriggerMask, isZero,
599  fMinIpZ, fMaxIpZ, data, dataMC,
600  fFilterMask, ipzW))
601  fSeenCent->Fill(cent);
602  }
603 
604  return taken;
605 }
606 
607 //________________________________________________________________________
609  TH2*,
610  TH2*)
611 {
612 }
613 
614 //________________________________________________________________________
615 void
617  const char* title, const char* ytitle)
618 {
619  //
620  // Set histogram graphical options, etc.
621  //
622  // Parameters:
623  // h Histogram to modify
624  // colour Marker color
625  // marker Marker style
626  // title Title of histogram
627  // ytitle Title on y-axis.
628  //
629  h->SetTitle(title);
630  h->SetMarkerColor(colour);
631  h->SetMarkerStyle(marker);
632  h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1);
633  h->SetFillStyle(0);
634  TString ytit;
635  if (ytitle && ytitle[0] != '\0') ytit = ytitle;
636  ytit = "#frac{1}{#it{N}}#frac{d#it{N}_{ch}}{d#it{#eta}}";
637  h->SetYTitle(ytit);
638  h->GetXaxis()->SetTitleFont(132);
639  h->GetXaxis()->SetLabelFont(132);
640  h->GetXaxis()->SetNdivisions(10);
641  h->GetYaxis()->SetTitleFont(132);
642  h->GetYaxis()->SetLabelFont(132);
643  h->GetYaxis()->SetNdivisions(10);
644  h->GetYaxis()->SetDecimals();
645  h->SetStats(0);
646 }
647 
648 //________________________________________________________________________
649 void
651 {
652  // Normalize to the acceptance -
653  // dndeta->Divide(accNorm);
654  for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
655  Double_t a = norm->GetBinContent(i);
656  for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
657  if (a <= 0) {
658  copy->SetBinContent(i,j,0);
659  copy->SetBinError(i,j,0);
660  continue;
661  }
662  Double_t c = copy->GetBinContent(i, j);
663  Double_t e = copy->GetBinError(i, j);
664  copy->SetBinContent(i, j, c / a);
665  copy->SetBinError(i, j, e / a);
666  }
667  }
668 }
669 //________________________________________________________________________
670 void
672 {
673  // Normalize to the acceptance -
674  // dndeta->Divide(accNorm);
675  for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
676  Double_t a = norm->GetBinContent(i);
677  if (a <= 0) {
678  copy->SetBinContent(i,0);
679  copy->SetBinError(i,0);
680  continue;
681  }
682  Double_t c = copy->GetBinContent(i);
683  Double_t e = copy->GetBinError(i);
684  copy->SetBinContent(i, c / a);
685  copy->SetBinError(i, e / a);
686  }
687 }
688 
689 //________________________________________________________________________
690 TH1D*
692  const char* name,
693  Int_t firstbin,
694  Int_t lastbin,
695  bool useRoot,
696  bool corr,
697  bool error)
698 {
699  //
700  // Project onto the X axis
701  //
702  // Parameters:
703  // h 2D histogram
704  // name New name
705  // firstbin First bin to use
706  // lastbin Last bin to use
707  // error Whether to calculate errors
708  //
709  // Return:
710  // Newly created histogram or null
711  //
712  if (!h) return 0;
713  if (useRoot)
714  return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
715 
716  const TAxis* xaxis = h->GetXaxis();
717  const TAxis* yaxis = h->GetYaxis();
718  TH1D* ret = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
719  xaxis->GetXmin(), xaxis->GetXmax());
720  static_cast<const TAttLine*>(h)->Copy(*ret);
721  static_cast<const TAttFill*>(h)->Copy(*ret);
722  static_cast<const TAttMarker*>(h)->Copy(*ret);
723  ret->GetXaxis()->ImportAttributes(xaxis);
724 
725  Int_t first = firstbin;
726  Int_t last = lastbin;
727  if (first < 0) first = 1;
728  else if (first >= yaxis->GetNbins()+2) first = yaxis->GetNbins()+1;
729  if (last < 0) last = yaxis->GetNbins();
730  else if (last >= yaxis->GetNbins()+2) last = yaxis->GetNbins()+1;
731  if (last-first < 0) {
732  AliWarningGeneral("AliBasedNdetaTask",
733  Form("Nothing to project [%d,%d]", first, last));
734  return 0;
735 
736  }
737 
738  // Loop over X bins
739  //DMSG(fDebug,3,"Projecting bins [%d,%d] of %s", first, last, h->GetName()));
740  Int_t ybins = (last-first+1);
741  for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
742  Double_t content = 0;
743  Double_t error2 = 0;
744  Int_t nbins = 0;
745 
746 
747  for (Int_t ybin = first; ybin <= last; ybin++) {
748  Double_t c1 = h->GetBinContent(h->GetBin(xbin, ybin));
749  Double_t e1 = h->GetBinError(h->GetBin(xbin, ybin));
750 
751  // Ignore empty bins
752  if (c1 < 1e-12) continue;
753  if (e1 < 1e-12) {
754  if (error) continue;
755  e1 = 1;
756  }
757 
758  content += c1;
759  error2 += e1*e1;
760  nbins++;
761  } // for (ybin)
762  if(content > 0 && nbins > 0) {
763  Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
764 #if 0
765  AliWarningGeneral(ret->GetName(),
766  Form("factor @ %d is %d/%d -> %f",
767  xbin, ybins, nbins, factor));
768 #endif
769  if (error) {
770  // Calculate weighted average
771  ret->SetBinContent(xbin, content * factor);
772  ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
773  }
774  else
775  ret->SetBinContent(xbin, factor * content);
776  }
777  } // for (xbin)
778 
779  return ret;
780 }
781 
782 //________________________________________________________________________
783 Bool_t
785 {
786  //
787  // Called at end of event processing..
788  //
789  // This is called once in the master
790  //
791  // Parameters:
792  // option Not used
793  //
794  // Draw result to screen, or perform fitting, normalizations Called
795  // once at the end of the query
796  DGUARD(fDebug,1,"Process final merged results");
797 
798  UShort_t sNN;
799  UShort_t sys;
800  ULong_t trig;
801  ULong_t filter;
802  UShort_t scheme;
803  Int_t centID;
804  AliForwardUtil::GetParameter(fSums->FindObject("sNN"), sNN);
805  AliForwardUtil::GetParameter(fSums->FindObject("sys"), sys);
806  AliForwardUtil::GetParameter(fSums->FindObject("trigger"), trig);
807  AliForwardUtil::GetParameter(fSums->FindObject("filter"), filter);
808  AliForwardUtil::GetParameter(fSums->FindObject("scheme"), scheme);
809  AliForwardUtil::GetParameter(fSums->FindObject("centEstimator"), centID);
810 
811  TH1* cH = static_cast<TH1*>(fSums->FindObject("centAxis"));
812  Info("", "centAxis: %p (%s)", cH, (cH ? cH->ClassName() : "null"));
813  // TAxis* cA = static_cast<TAxis*>(fSums->FindObject("centAxis"));
814  TAxis* cA = (cH ? cH->GetXaxis() : 0);
815  if (cA) cA->Copy(fCentAxis);
816  fCentAxis.SetName("centAxis");
817  fCentAxis.SetTitle("Centrality [%]");
818 
819  // (Re)Create our centrality bins
821 
822  // Print before we loop
823  Print();
824 
825  // Loop over centrality bins
826  TIter next(fListOfCentralities);
827  CentralityBin* bin = 0;
828  gStyle->SetPalette(1);
829  THStack* dndetaStack = new THStack("dndeta", "dN/d#eta");
830  THStack* dndetaMCStack = new THStack("dndetaMC", "dN_{ch}/d#eta");
831  THStack* dndetaEmpStack = new THStack("dndetaEmp", "dN_{ch}/d#eta");
832  THStack* leftRightStack = new THStack("leftRight", "#eta>0/#eta<0");
833 
834  TList* mclist = 0;
835  TList* truthlist = 0;
836 
837  if (fFinalMCCorrFile.Contains(".root")) {
838  TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
839  if(ftest) {
840  mclist = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
841  truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
842  }
843  else
844  AliWarning("MC analysis file invalid - no final MC correction possible");
845  }
846  Int_t style = GetMarker();
847  Int_t color = GetColor();
848 
849  DMSG(fDebug,3,"Marker style=%d, color=%d", style, color);
850  while ((bin = static_cast<CentralityBin*>(next()))) {
853  style, color, mclist, truthlist);
854  if (HasCentrality() && bin->IsAllBin())
855  // If we have centrality bins, do not add the min-bias
856  // distribution to the output stack.
857  continue;
858  TH1* dndeta = bin->GetResult("");
859  TH1* dndetaMC = bin->GetResult("MC", false);
860  TH1* dndetaEmp = bin->GetResult("Emp", false);
861  TH1* leftRight = Asymmetry(dndetaEmp ? dndetaEmp : dndeta);
862  // if (leftRight) bin->fOutput->Add(leftRight);
863  DMSG(fDebug,2,"Results: bare=%p mcbare=%p", dndeta, dndetaMC);
864  if (dndeta) dndetaStack->Add(dndeta);
865  if (dndetaMC) dndetaMCStack->Add(dndetaMC);
866  if (dndetaEmp) dndetaEmpStack->Add(dndetaEmp);
867  if (leftRight) leftRightStack->Add(leftRight);
868  }
869  // Output the stack
870  fResults->Add(dndetaStack);
871 
872  // If available, output track-ref stack
873  if (!dndetaMCStack->GetHists() ||
874  dndetaMCStack->GetHists()->GetEntries() <= 0) {
875  // AliWarning("No MC histograms found");
876  delete dndetaMCStack;
877  dndetaMCStack = 0;
878  }
879  if (dndetaMCStack) fResults->Add(dndetaMCStack);
880 
881  // If available, output track-ref stack
882  if (!dndetaEmpStack->GetHists() ||
883  dndetaEmpStack->GetHists()->GetEntries() <= 0) {
884  // AliWarning("No MC histograms found");
885  delete dndetaEmpStack;
886  dndetaEmpStack = 0;
887  }
888  if (dndetaEmpStack) fResults->Add(dndetaEmpStack);
889 
890 
891  // If available, output track-ref stack
892  if (!leftRightStack->GetHists() ||
893  leftRightStack->GetHists()->GetEntries() <= 0) {
894  // AliWarning("No MC histograms found");
895  delete leftRightStack;
896  leftRightStack = 0;
897  }
898  if (leftRightStack) fResults->Add(leftRightStack);
899 
900  // Output collision energy string
901  if (sNN > 0) {
902  TNamed* sNNObj = new TNamed("sNN",
904  sNNObj->SetUniqueID(sNN);
905  fResults->Add(sNNObj);
906  }
907 
908  // Output collision system string
909  if (sys > 0) {
910  TNamed* sysObj = new TNamed("sys",
912  sysObj->SetUniqueID(sys);
913  fResults->Add(sysObj);
914  }
915 
916  // Output centrality axis
917  TNamed* centEstimator = new TNamed("centEstimator", fCentMethod.Data());
918  centEstimator->SetUniqueID(centID);
919  fResults->Add(&fCentAxis);
920  fResults->Add(centEstimator);
921 
922  // Output trigger string
923  if (trig) {
924  TNamed* maskObj = new TNamed("trigger",
926  maskObj->SetUniqueID(trig);
927  fResults->Add(maskObj);
928  }
929  // Output filter string
930  if (filter) {
931  TNamed* maskObj = new TNamed("filter",
932  AliAODForwardMult::GetTriggerString(filter, " | "));
933  maskObj->SetUniqueID(filter);
934  fResults->Add(maskObj);
935  }
936 
937  // Normalization string
938  if (scheme > 0) {
939  TNamed* schemeObj = new TNamed("scheme",
940  NormalizationSchemeString(scheme));
941  schemeObj->SetUniqueID(scheme);
942  fResults->Add(schemeObj);
943  }
944 
945  // Output vertex axis
946  TAxis* vtxAxis = new TAxis(1,fMinIpZ,fMaxIpZ);
947  vtxAxis->SetName("vtxAxis");
948  vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fMinIpZ,fMaxIpZ));
949  fResults->Add(vtxAxis);
950 
951  // Output trigger efficiency
954 
955  TNamed* options = new TNamed("options","");
956  TString str;
957  str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
958  str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
959  options->SetTitle(str);
960  fResults->Add(options);
961 
962  return true;
963 }
964 
965 #define PF(N,V,...) \
966  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
967 #define PFB(N,FLAG) \
968  do { \
969  AliForwardUtil::PrintName(N); \
970  std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
971  } while(false)
972 #define PFV(N,VALUE) \
973  do { \
974  AliForwardUtil::PrintName(N); \
975  std::cout << (VALUE) << std::endl; } while(false)
976 #if 0
977 namespace {
978  void appendBit(TString& str, const char* bit)
979  {
980  if (!str.IsNull()) str.Append("|");
981  str.Append(bit);
982  }
983 }
984 #endif
985 
986 //________________________________________________________________________
987 void
989 {
990  //
991  // Print information
992  //
993  AliBaseAODTask::Print(option);
995 
996  gROOT->IncreaseDirLevel();
997  PFB("Use AnaUtil pile-up", fUseUtilPileup);
998  PFB("Use TH2::ProjectionX", fUseROOTProj);
999  PFB("Correct for empty", fCorrEmpty);
1000  PFV("Normalization scheme", schemeString );
1001  PFV("Trigger efficiency", fTriggerEff);
1002  PFV("Bin-0 Trigger efficiency", fTriggerEff0);
1003  PFV("Centrality estimator", (fCentMethod.IsNull() ?
1004  "-default-" : fCentMethod.Data()));
1005  PFB("Re-weight from IPz", fIpzReweight!=0);
1006  if (fIpzReweight) fIpzReweight->Print();
1007 
1008 
1009  TString opt(option);
1010  opt.ToUpper();
1011  if (opt.Contains("R") &&
1013  fListOfCentralities->GetEntries() > 0) {
1014  TIter next(fListOfCentralities);
1015  TObject* bin = 0;
1016  while ((bin = next())) bin->Print(option);
1017  }
1018  gROOT->DecreaseDirLevel();
1019 }
1020 
1021 //__________________________________________________________________
1022 Int_t
1024 {
1025  Int_t base = bits & (0xFE);
1026  Bool_t hollow = bits & kHollow;
1027  switch (base) {
1028  case kCircle: return (hollow ? 24 : 20);
1029  case kSquare: return (hollow ? 25 : 21);
1030  case kUpTriangle: return (hollow ? 26 : 22);
1031  case kDownTriangle: return (hollow ? 32 : 23);
1032  case kDiamond: return (hollow ? 27 : 33);
1033  case kCross: return (hollow ? 28 : 34);
1034  case kStar: return (hollow ? 30 : 29);
1035  }
1036  return 1;
1037 }
1038 //__________________________________________________________________
1039 UShort_t
1041 {
1042  UShort_t bits = 0;
1043  switch (style) {
1044  case 24: case 25: case 26: case 27: case 28: case 30: case 32:
1045  bits |= kHollow; break;
1046  }
1047  switch (style) {
1048  case 20: case 24: bits |= kCircle; break;
1049  case 21: case 25: bits |= kSquare; break;
1050  case 22: case 26: bits |= kUpTriangle; break;
1051  case 23: case 32: bits |= kDownTriangle; break;
1052  case 27: case 33: bits |= kDiamond; break;
1053  case 28: case 34: bits |= kCross; break;
1054  case 29: case 30: bits |= kStar; break;
1055  }
1056  return bits;
1057 }
1058 //__________________________________________________________________
1059 Int_t
1061 {
1062  UShort_t bits = GetMarkerBits(style);
1063  Int_t ret = GetMarkerStyle(bits ^ kHollow);
1064  return ret;
1065 }
1066 
1067 //====================================================================
1068 void
1070 {
1071  DGUARD(fDebug,1,"Initializing sums with %s", data->GetName());
1072  TString n(GetHistName(0));
1073  TString n0(GetHistName(1));
1074  const char* postfix = GetTitle();
1075 
1076  fSum = static_cast<TH2D*>(data->Clone(n));
1077  if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix));
1078  fSum->SetDirectory(0);
1079  fSum->SetMarkerColor(col);
1080  fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid));
1081  fSum->Reset();
1082  list->Add(fSum);
1083 
1084  fSum0 = static_cast<TH2D*>(data->Clone(n0));
1085  if (postfix)
1086  fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix));
1087  else
1088  fSum0->SetTitle(Form("%s 0-bin", data->GetTitle()));
1089  fSum0->SetDirectory(0);
1090  fSum0->SetMarkerColor(col);
1091  fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow));
1092  fSum0->Reset();
1093  list->Add(fSum0);
1094 
1095  fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5);
1096  fEvents->SetDirectory(0);
1097  fEvents->SetFillColor(kRed+1);
1098  fEvents->SetFillStyle(3002);
1099  fEvents->GetXaxis()->SetBinLabel(1, "Non-zero");
1100  fEvents->GetXaxis()->SetBinLabel(2, "Zero");
1101  list->Add(fEvents);
1102 }
1103 
1104 //____________________________________________________________________
1105 TString
1107  Int_t what, const char* post)
1108 {
1109  TString n;
1110  switch (what) {
1111  case 0: n = "sum"; break;
1112  case 1: n = "sum0"; break;
1113  case 2: n = "events"; break;
1114  }
1115  if (post && post[0] != '\0') n.Append(post);
1116  return n;
1117 }
1118 
1119 //____________________________________________________________________
1120 TString
1122 {
1123  return GetHistName(GetName(), what, GetTitle());
1124 }
1125 
1126 //____________________________________________________________________
1127 void
1128 AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero, Double_t weight)
1129 {
1130  DGUARD(fDebug,2,"Adding %s to sums", data->GetName());
1131  if (isZero) fSum0->Add(data, weight);
1132  else fSum->Add(data, weight);
1133  fEvents->Fill(isZero ? 1 : 0);
1134 }
1135 
1136 //____________________________________________________________________
1137 TH2D*
1139  Double_t& ntotal,
1140  Double_t epsilon0,
1141  Double_t epsilon,
1142  Int_t marker,
1143  Bool_t rootProj,
1144  Bool_t corrEmpty) const
1145 {
1146  DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName());
1147 
1148  // The return value `ret' is not scaled in anyway
1149  TH2D* ret = static_cast<TH2D*>(fSum->Clone(fSum->GetName()));
1150  ret->SetDirectory(0);
1151  Int_t n = Int_t(fEvents->GetBinContent(1));
1152  Int_t n0 = Int_t(fEvents->GetBinContent(2));
1153  ntotal = n;
1154  if (n0 > 0) {
1155  ret->Reset();
1156  DMSG(fDebug,1,
1157  "Adding histograms %s(%d) and %s(%d) w/weights %f and %f resp.",
1158  fSum0->GetName(), n, fSum->GetName(), n0, 1./epsilon,1./epsilon0);
1159  ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon);
1160  ntotal = n / epsilon + n0 / epsilon0;
1161  }
1162 
1163  TList* out = new TList;
1164  out->SetOwner();
1165  const char* postfix = GetTitle();
1166  if (!postfix) postfix = "";
1167  out->SetName(Form("partial%s", postfix));
1168  output->Add(out);
1169 
1170  // Now make copies, normalize them, and store in output list
1171  // Note, these are the only ones normalized here
1172  // These are mainly for diagnostics
1173  TH2D* sumCopy = static_cast<TH2D*>(fSum->Clone("sum"));
1174  TH2D* sum0Copy = static_cast<TH2D*>(fSum0->Clone("sum0"));
1175  TH2D* retCopy = static_cast<TH2D*>(ret->Clone("sumAll"));
1176  sumCopy->SetMarkerStyle(FlipHollowStyle(marker));
1177  sumCopy->SetDirectory(0);
1178  sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4));
1179  sum0Copy->SetDirectory(0);
1180  retCopy->SetMarkerStyle(marker);
1181  retCopy->SetDirectory(0);
1182 
1183  Int_t nY = fSum->GetNbinsY();
1184  Int_t o = 0; // nY+1;
1185  TH1D* norm = ProjectX(fSum, "norm", o, o, rootProj, corrEmpty, false);
1186  TH1D* norm0 = ProjectX(fSum0, "norm0", o, o, rootProj, corrEmpty, false);
1187  TH1D* normAll = ProjectX(ret, "normAll", o, o, rootProj, corrEmpty, false);
1188  norm->SetTitle("#eta coverage - >0-bin");
1189  norm0->SetTitle("#eta coverage - 0-bin");
1190  normAll->SetTitle("#eta coverage");
1191  norm->SetDirectory(0);
1192  norm0->SetDirectory(0);
1193  normAll->SetDirectory(0);
1194 
1195  TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1,nY,rootProj,corrEmpty);
1196  TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1,nY,rootProj,corrEmpty);
1197  TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1,nY,rootProj,corrEmpty);
1198  sumCopyPx-> SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sumCopy->GetTitle()));
1199  sum0CopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sum0Copy->GetTitle()));
1200  retCopyPx-> SetTitle(Form("#sum_{i}^{N_{#phi}}%s", retCopy->GetTitle()));
1201  sumCopyPx-> SetDirectory(0);
1202  sum0CopyPx->SetDirectory(0);
1203  retCopyPx-> SetDirectory(0);
1204 
1205  TH1D* phi = ProjectX(fSum, "phi", nY+1,nY+1,rootProj,corrEmpty,false);
1206  TH1D* phi0 = ProjectX(fSum0, "phi0", nY+1,nY+1,rootProj,corrEmpty,false);
1207  TH1D* phiAll = ProjectX(ret, "phiAll", nY+1,nY+1,rootProj,corrEmpty,false);
1208  phi ->SetTitle("#phi acceptance from dead strips - >0-bin");
1209  phi0 ->SetTitle("#phi acceptance from dead strips - 0-bin");
1210  phiAll->SetTitle("#phi acceptance from dead strips");
1211  phi ->SetDirectory(0);
1212  phi0 ->SetDirectory(0);
1213  phiAll->SetDirectory(0);
1214 
1215  const TH1D* cov = (corrEmpty ? norm : phi);
1216  const TH1D* cov0 = (corrEmpty ? norm0 : phi0);
1217  const TH1D* covAll = (corrEmpty ? normAll : phiAll);
1218 
1219  // Here, we scale to the coverage (or phi acceptance)
1220  ScaleToCoverage(sumCopy, cov);
1221  ScaleToCoverage(sum0Copy, cov0);
1222  ScaleToCoverage(retCopy, covAll);
1223 
1224  // Scale our 1D histograms
1225  sumCopyPx ->Scale(1., "width");
1226  sum0CopyPx->Scale(1., "width");
1227  retCopyPx ->Scale(1., "width");
1228 
1229  DMSG(fDebug,2,"Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(),
1230  sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum());
1231 
1232  // Scale the normalization - they should be 1 at the maximum
1233  norm ->Scale(n > 0 ? 1. / n : 1);
1234  norm0 ->Scale(n0 > 0 ? 1. / n0 : 1);
1235  normAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1236 
1237  // Scale the normalization - they should be 1 at the maximum
1238  phi ->Scale(n > 0 ? 1. / n : 1);
1239  phi0 ->Scale(n0 > 0 ? 1. / n0 : 1);
1240  phiAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1241 
1242  out->Add(sumCopy);
1243  out->Add(sum0Copy);
1244  out->Add(retCopy);
1245  out->Add(sumCopyPx);
1246  out->Add(sum0CopyPx);
1247  out->Add(retCopyPx);
1248  out->Add(norm);
1249  out->Add(norm0);
1250  out->Add(normAll);
1251  out->Add(phi);
1252  out->Add(phi0);
1253  out->Add(phiAll);
1254 
1255  if (fDebug >= 1) {
1256  if (n0 > 0)
1257  DMSG(fDebug,1,"Returning (1/%f * %s + 1/%f * %s), "
1258  "1/%f * %d + 1/%f * %d = %d",
1259  epsilon0, fSum0->GetName(), epsilon, fSum->GetName(),
1260  epsilon0, n0, epsilon, n, int(ntotal));
1261  else
1262  DMSG(fDebug,1, "Returning %s, %d", fSum->GetName(), int(ntotal));
1263  }
1264 
1265 #if 0
1266  for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
1267  Double_t nc = sum->GetBinContent(i, 0);
1268  Double_t nc0 = sum0->GetBinContent(i, 0);
1269  ret->SetBinContent(i, 0, nc + nc0); // Just count events
1270  }
1271 #endif
1272 
1273  return ret;
1274 }
1275 //____________________________________________________________________
1276 void
1278 {
1279  PFV("dN/deta sum bin", GetName());
1280  gROOT->IncreaseDirLevel();
1281  PF("Normal sum", "%s (%p)", GetHistName(0).Data(), fSum);
1282  PF("0-bin sum", "%s (%p)", GetHistName(1).Data(), fSum0);
1283  PF("Event count","%s (%p)", GetHistName(2).Data(), fEvents);
1284  gROOT->DecreaseDirLevel();
1285 }
1286 
1287 //====================================================================
1289  : TNamed("", ""),
1290  fSums(0),
1291  fOutput(0),
1292  fSum(0),
1293  fSumMC(0),
1294  fTriggers(0),
1295  fStatus(0),
1296  fLow(0),
1297  fHigh(0),
1298  fDoFinalMCCorrection(false),
1299  fSatelliteVertices(false),
1300  fDebug(0)
1301 {
1302  //
1303  // Constructor
1304  //
1305  DGUARD(fDebug,3,"Default CTOR of AliBasedNdeta::CentralityBin");
1306 }
1307 #define TRUNC(X) (Int_t(X) + Float_t(Int_t(X*100)%100)/100)
1308 
1309 //____________________________________________________________________
1311  Float_t low, Float_t high)
1312  : TNamed(name, ""),
1313  fSums(0),
1314  fOutput(0),
1315  fSum(0),
1316  fSumMC(0),
1317  fTriggers(0),
1318  fStatus(0),
1319  fLow(TRUNC(low)),
1320  fHigh(TRUNC(high)),
1321  fDoFinalMCCorrection(false),
1322  fSatelliteVertices(false),
1323  fDebug(0)
1324 {
1325  //
1326  // Constructor
1327  //
1328  // Parameters:
1329  // name Name used for histograms (e.g., Forward)
1330  // low Lower centrality cut in percent
1331  // high Upper centrality cut in percent
1332  //
1333  DGUARD(fDebug,3,"Named CTOR of AliBasedNdeta::CentralityBin: "
1334  "%s [%6.2f,%6.2f]",name, fLow, fHigh);
1335  if (low <= 0 && high <= 0) {
1336  fLow = 0;
1337  fHigh = 0;
1338  SetTitle("All centralities");
1339  }
1340  else {
1341  fLow = TRUNC(low);
1342  fHigh = TRUNC(high);
1343  SetTitle(Form("Centrality bin from %6.2f%% to %6.2f%%", low, high));
1344  }
1345 }
1346 //____________________________________________________________________
1348  : TNamed(o),
1349  fSums(o.fSums),
1350  fOutput(o.fOutput),
1351  fSum(o.fSum),
1352  fSumMC(o.fSumMC),
1353  fTriggers(o.fTriggers),
1354  fStatus(o.fStatus),
1355  fLow(o.fLow),
1356  fHigh(o.fHigh),
1357  fDoFinalMCCorrection(o.fDoFinalMCCorrection),
1359  fDebug(o.fDebug)
1360 {
1361  //
1362  // Copy constructor
1363  //
1364  // Parameters:
1365  // other Object to copy from
1366  //
1367  DGUARD(fDebug,3,"Copy CTOR of AliBasedNdeta::CentralityBin");
1368 }
1369 //____________________________________________________________________
1371 {
1372  //
1373  // Destructor
1374  //
1375  DGUARD(fDebug,3,"DTOR of AliBasedNdeta::CentralityBin");
1376 
1377  // if (fSums) fSums->Delete();
1378  // if (fOutput) fOutput->Delete();
1379 }
1380 
1381 //____________________________________________________________________
1384 {
1385  //
1386  // Assignment operator
1387  //
1388  // Parameters:
1389  // other Object to assign from
1390  //
1391  // Return:
1392  // Reference to this
1393  //
1394  DGUARD(fDebug,3,"Centrality bin assignment");
1395  if (&o == this) return *this;
1396  SetName(o.GetName());
1397  SetTitle(o.GetTitle());
1398  fSums = o.fSums;
1399  fOutput = o.fOutput;
1400  fSum = o.fSum;
1401  fSumMC = o.fSumMC;
1402  fTriggers = o.fTriggers;
1403  fStatus = o.fStatus;
1404  fLow = o.fLow;
1405  fHigh = o.fHigh;
1406  fDoFinalMCCorrection = o.fDoFinalMCCorrection;
1408 
1409  return *this;
1410 }
1411 #if 0
1412 namespace {
1413  Color_t Brighten(Color_t origNum, Int_t nTimes=2)
1414  {
1415  TColor* col = gROOT->GetColor(origNum);
1416  if (!col) return origNum;
1417  Int_t origR = Int_t(0xFF * col->GetRed());
1418  Int_t origG = Int_t(0xFF * col->GetGreen());
1419  Int_t origB = Int_t(0xFF * col->GetBlue());
1420  Int_t off = nTimes*0x33;
1421  Int_t newR = TMath::Min((origR+off),0xff);
1422  Int_t newG = TMath::Min((origG+off),0xff);
1423  Int_t newB = TMath::Min((origB+off),0xff);
1424  Int_t newNum = TColor::GetColor(newR, newG, newB);
1425  return newNum;
1426  }
1427 }
1428 #endif
1429 //____________________________________________________________________
1430 Int_t
1432 {
1433  if (IsAllBin()) return fallback;
1434  Float_t fc = /*(fLow+double(fHigh-fLow)/2)*/ fHigh / 100;
1435  Int_t nCol = gStyle->GetNumberOfColors();
1436  Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
1437  Int_t col = gStyle->GetColorPalette(icol);
1438 #if 0
1439  Color_t orig = col;
1440  col = Brighten(orig);
1441 #endif
1442  return col;
1443 }
1444 
1445 //____________________________________________________________________
1446 const char*
1448 {
1449  //
1450  // Get the list name
1451  //
1452  // Return:
1453  // List Name
1454  //
1455  if (IsAllBin()) return "all";
1456  return Form("cent%03dd%02d_%03dd%02d",
1457  Int_t(fLow), Int_t(fLow*100) % 100,
1458  Int_t(fHigh), Int_t(fHigh*100) % 100);
1459 
1460 }
1461 //____________________________________________________________________
1462 void
1464 {
1465  //
1466  // Create output objects
1467  //
1468  // Parameters:
1469  // dir Parent list
1470  //
1471  DGUARD(fDebug,1,"Create centrality bin output objects");
1472  fSums = new TList;
1473  fSums->SetName(GetListName());
1474  fSums->SetOwner();
1475  dir->Add(fSums);
1476 
1478  fTriggers->SetDirectory(0);
1479 
1480  fStatus = AliAODForwardMult::MakeStatusHistogram("status");
1481  fStatus->SetDirectory(0);
1482 
1483  fSums->Add(fTriggers);
1484  fSums->Add(fStatus);
1485 
1486  fSums->Add(new TParameter<float>("low", fLow, 'f'));
1487  fSums->Add(new TParameter<float>("high", fHigh, 'f'));
1488 }
1489 //____________________________________________________________________
1490 void
1492 {
1493  fDebug = lvl;
1494  if (fSum) fSum->fDebug = lvl;
1495  if (fSumMC) fSumMC->fDebug = lvl;
1496 }
1497 
1498 //____________________________________________________________________
1499 Bool_t
1501 {
1502  const char* post = (mc ? "MC" : "");
1503  TString sn = Sum::GetHistName(GetName(),0,post);
1504  TString sn0 = Sum::GetHistName(GetName(),1,post);
1505  TString ev = Sum::GetHistName(GetName(),2,post);
1506  TH2D* sum = static_cast<TH2D*>(list->FindObject(sn));
1507  TH2D* sum0 = static_cast<TH2D*>(list->FindObject(sn0));
1508  TH1I* events = static_cast<TH1I*>(list->FindObject(ev));
1509  if (!sum || !sum0 || !events) {
1510  if (!mc)
1511  AliWarningF("Failed to find one or more histograms: "
1512  "%s (%p) %s (%p) %s (%p)",
1513  sn.Data(), sum,
1514  sn0.Data(), sum0,
1515  ev.Data(), events);
1516  return false;
1517  }
1518  Sum* ret = new Sum(GetName(), post);
1519  ret->fSum = sum;
1520  ret->fSum0 = sum0;
1521  ret->fEvents = events;
1522  ret->fDebug = fDebug;
1523  if (mc) fSumMC = ret;
1524  else fSum = ret;
1525 
1526  return true;
1527 }
1528 
1529 //____________________________________________________________________
1530 void
1532 {
1533  //
1534  // Create sum histogram
1535  //
1536  // Parameters:
1537  // data Data histogram to clone
1538  // mc (optional) MC histogram to clone
1539  //
1540  DGUARD(fDebug,1,"Create centrality bin sums from %s",
1541  data ? data->GetName() : "(null)");
1542  if (data) {
1543  fSum = new Sum(GetName(),"");
1544  fSum->Init(fSums, data, GetColor());
1545  fSum->fDebug = fDebug;
1546  }
1547 
1548  // If no MC data is given, then do not create MC sum histogram
1549  if (!mc) return;
1550 
1551  fSumMC = new Sum(GetName(), "MC");
1552  fSumMC->Init(fSums, mc, GetColor());
1553  fSumMC->fDebug = fDebug;
1554 }
1555 
1556 //____________________________________________________________________
1557 Bool_t
1559  Int_t triggerMask,
1560  Double_t vzMin,
1561  Double_t vzMax,
1562  Int_t filter)
1563 {
1564  //
1565  // Check the trigger, vertex, and centrality
1566  //
1567  // Parameters:
1568  // aod Event input
1569  //
1570  // Return:
1571  // true if the event is to be used
1572  //
1573  if (!forward) return false;
1574 
1575  DGUARD(fDebug,2,"Check the event");
1576  // We do not check for centrality here - it's already done
1577  Bool_t ret = forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0,
1578  fTriggers, fStatus, filter);
1579  DMSG(fDebug, 2, "%s", (ret ? "Accepted" : "Rejected"));
1580  return ret;
1581 }
1582 
1583 
1584 //____________________________________________________________________
1585 Bool_t
1587  UInt_t triggerMask,
1588  Bool_t isZero,
1589  Double_t vzMin,
1590  Double_t vzMax,
1591  const TH2D* data,
1592  const TH2D* mc,
1593  UInt_t filter,
1594  Double_t weight)
1595 {
1596  //
1597  // Process an event
1598  //
1599  // Parameters:
1600  // forward Forward data (for trigger, vertex, & centrality)
1601  // triggerMask Trigger mask
1602  // vzMin Minimum IP z coordinate
1603  // vzMax Maximum IP z coordinate
1604  // data Data histogram
1605  // mc MC histogram
1606  // weight Event weight
1607  //
1608  DGUARD(fDebug,1,"Process one event for %s a given centrality bin "
1609  "[%5.1f%%,%5.1f%%)", data ? data->GetName() : "(null)", fLow, fHigh);
1610  if (!CheckEvent(forward, triggerMask, vzMin, vzMax, filter))
1611  return false;
1612  if (!data) return false;
1613  if (!fSum) CreateSums(data, mc);
1614 
1615  fSum->Add(data, isZero, weight);
1616  if (mc) fSumMC->Add(mc, isZero, weight);
1617 
1618  return true;
1619 }
1620 
1621 //________________________________________________________________________
1622 Double_t
1624  UShort_t scheme,
1625  Double_t trigEff,
1626  Double_t& ntotal,
1627  TString* text) const
1628 {
1629  //
1630  // Calculate normalization
1631  //
1632  // Parameters:
1633  // t Trigger histogram
1634  // scheme Normaliztion scheme
1635  // trigEff From MC
1636  // ntotal On return, contains the number of events.
1637  //
1638  DGUARD(fDebug,1,"Normalize centrality bin %s [%6.2f-%6.2f%%] with %s",
1639  GetName(), fLow, fHigh, t.GetName());
1640  Double_t nAll = t.GetBinContent(AliAODForwardMult::kBinAll);
1641  Double_t nB = t.GetBinContent(AliAODForwardMult::kBinB);
1642  Double_t nA = t.GetBinContent(AliAODForwardMult::kBinA);
1643  Double_t nC = t.GetBinContent(AliAODForwardMult::kBinC);
1644  Double_t nE = t.GetBinContent(AliAODForwardMult::kBinE);
1645  Double_t nOffline = t.GetBinContent(AliAODForwardMult::kBinOffline);
1646  Double_t nTriggered = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1647  Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1648  Double_t nAccepted = ntotal;
1649  ntotal = 0;
1650 
1651  if (nTriggered <= 0.1) {
1652  AliError("Number of triggered events <= 0");
1653  return -1;
1654  }
1655  if (nWithVertex <= 0.1) {
1656  AliError("Number of events with vertex <= 0");
1657  return -1;
1658  }
1659  ntotal = nAccepted;
1660  Double_t vtxEff = nWithVertex / nTriggered;
1661  Double_t scaler = 1;
1662  Double_t beta = nA + nC - 2*nE;
1663 
1664 
1665  TString rhs("N = N_acc");
1666  if (!(scheme & kZeroBin)) {
1667  if (scheme & kEventLevel) {
1668  ntotal = nAccepted / vtxEff;
1669  scaler = vtxEff;
1670  DMSG(fDebug,0,"Calculating event normalisation as\n"
1671  " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1672  Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex),
1673  ntotal, scaler);
1674  if (scheme & kBackground) {
1675  // 1 E_V E_V
1676  // s = --------- = ------------- = ------------
1677  // 1 - beta 1 - beta E_V 1 - beta N_V
1678  // --- ---- -------- ---- ---
1679  // E_V N_V N_V N_V N_T
1680  //
1681  // E_V
1682  // = ------------
1683  // 1 - beta
1684  // ----
1685  // N_T
1686  //
1687  ntotal -= nAccepted * beta / nWithVertex;
1688  // This one is direct and correct.
1689  // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1690  // A simpler expresion
1691  scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1692  DMSG(fDebug,0,"Calculating event normalisation as\n"
1693  " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1694  " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1695  Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1696  nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta),
1697  Int_t(nWithVertex), ntotal, scaler);
1698  rhs.Append("(1/eps_V - beta/N_vtx)");
1699  } // Background
1700  else
1701  rhs.Append("/eps_V");
1702  } // Event-level
1703  if (scheme & kTriggerEfficiency) {
1704  Int_t old = Int_t(ntotal);
1705  ntotal /= trigEff;
1706  scaler *= trigEff;
1707  DMSG(fDebug,0,"Correcting for trigger efficiency:\n"
1708  " N = 1 / E_X * N = 1 / %f * %d = %f (%f)",
1709  trigEff, old, ntotal, scaler);
1710  rhs.Append("/eps_T");
1711  } // Trigger efficiency
1712  }
1713  else {
1714  // Calculate as
1715  //
1716  // N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta)
1717  // = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1718  // = N_A (1 + 1/E_X (1/E_V - 1 - beta / N_V))
1719  //
1720  // s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1721  // = N_V / (N_V + 1/E_X (N_T - N_V - beta))
1722  //
1723  if (!(scheme & kBackground)) beta = 0;
1724  ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1
1725  - beta / nWithVertex));
1726  scaler = nWithVertex / (nWithVertex +
1727  1/trigEff * (nTriggered-nWithVertex-beta));
1728  DMSG(fDebug,0,"Calculating event normalisation as\n"
1729  " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1730  " N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = "
1731  "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)",
1732  Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1733  Int_t(nAccepted), trigEff, Int_t(nTriggered),
1734  Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex),
1735  ntotal, scaler);
1736  rhs.Append("(1+1/eps_T(1/eps_V-1-beta/N_vtx))");
1737  }
1738 
1739  if (text) {
1740  text->Append(Form("%-40s = %d\n", "N_all", UInt_t(nAll)));
1741  text->Append(Form("%-40s = %d\n", "N_acc", UInt_t(nAccepted)));
1742  text->Append(Form("%-40s = %d\n", "N_trg", UInt_t(nTriggered)));
1743  text->Append(Form("%-40s = %d\n", "N_vtx", UInt_t(nWithVertex)));
1744  text->Append(Form("%-40s = %d\n", "N_B", UInt_t(nB)));
1745  text->Append(Form("%-40s = %d\n", "N_A", UInt_t(nA)));
1746  text->Append(Form("%-40s = %d\n", "N_C", UInt_t(nC)));
1747  text->Append(Form("%-40s = %d\n", "N_E", UInt_t(nE)));
1748  text->Append(Form("%-40s = %d\n", "beta = N_A + N_C - 2N_E",UInt_t(beta)));
1749  text->Append(Form("%-40s = %f\n", "eps_V = N_vtx/N_trg",vtxEff));
1750  text->Append(Form("%-40s = %f\n", "eps_T", trigEff));
1751  text->Append(Form("%-40s = %f\n", rhs.Data(), ntotal));
1752  }
1753  TString tN = t.GetXaxis()->GetBinLabel(AliAODForwardMult::kWithTrigger);
1754  tN.ReplaceAll("w/Selected trigger (","");
1755  tN.ReplaceAll(")", "");
1756  DMSG(fDebug,0,"\n"
1757  " Total of %9d events for %s\n"
1758  " of these %9d have an offline trigger\n"
1759  " of these N_T = %9d has the selected trigger (%s)\n"
1760  " of these N_V = %9d has a vertex\n"
1761  " of these N_A = %9d were in the selected range\n"
1762  " Triggers by hardware type:\n"
1763  " N_b = %9d\n"
1764  " N_ac = %9d (%d+%d)\n"
1765  " N_e = %9d\n"
1766  " Vertex efficiency: %f\n"
1767  " Trigger efficiency: %f\n"
1768  " Total number of events: N = %f\n"
1769  " Scaler (N_A/N): %f\n"
1770  " %25s = %f",
1771  Int_t(nAll), GetTitle(), Int_t(nOffline),
1772  Int_t(nTriggered), tN.Data(),
1773  Int_t(nWithVertex), Int_t(nAccepted),
1774  Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE),
1775  vtxEff, trigEff, ntotal, scaler, rhs.Data(), ntotal);
1776  return scaler;
1777 }
1778 
1779 //________________________________________________________________________
1780 const char*
1782 {
1783  static TString n;
1784  n = GetName();
1785  n.ReplaceAll("dNdeta", "");
1786  n.Prepend("dndeta");
1787  n.Append(postfix);
1788  return n.Data();
1789 }
1790 //________________________________________________________________________
1791 TH1*
1793  Bool_t verbose) const
1794 {
1795  if (!fOutput) {
1796  AliWarningF("No output list defined in %s [%6.2f,%6.2f]", GetName(),
1797  fLow, fHigh);
1798  return 0;
1799  }
1800  TString n = GetResultName(postfix);
1801  TObject* o = fOutput->FindObject(n.Data());
1802  if (!o) {
1803  if (verbose)
1804  AliWarningF("Object %s not found in output list of %s",
1805  n.Data(), GetName());
1806  return 0;
1807  }
1808  return static_cast<TH1*>(o);
1809 }
1810 
1811 //________________________________________________________________________
1812 void
1814  const char* postfix,
1815  bool rootProj,
1816  bool corrEmpty,
1817  Double_t scaler,
1818  Int_t marker,
1819  Int_t color,
1820  TList* mclist,
1821  TList* truthlist)
1822 {
1823  //
1824  // Generate the dN/deta result from input
1825  //
1826  // Parameters:
1827  // sum Sum of 2D hists
1828  // postfix Post fix on names
1829  // rootProj Whether to use ROOT TH2::ProjectionX
1830  // corrEmpty Correct for empty bins
1831  // scaler Event-level normalization scaler
1832  //
1833  DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
1834  TString base(GetName());
1835  base.ReplaceAll("dNdeta", "");
1836  base.Append(postfix);
1837  TH2D* copy = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s",
1838  base.Data())));
1839 
1840  TH1D* accNorm = 0;
1841  Int_t nY = sum->GetNbinsY();
1842  // Hack HHD Hans test code to manually remove FMD2 dead channel (but
1843  // it is on outer?)
1844  //
1845  // cholm comment: The original hack has been moved to
1846  // AliForwarddNdetaTask::CheckEventData - this simplifies things a
1847  // great deal, and we could in principle use the new phi acceptance.
1848  //
1849  // However, since we may have filtered out the dead sectors in the
1850  // AOD production already, we can't be sure we can recalculate the
1851  // phi acceptance correctly, so for now, we rely on fCorrEmpty being set.
1852  Int_t o = (corrEmpty ? 0 : nY+1);
1853  accNorm = ProjectX(sum, Form("norm%s",base.Data()), o, o,
1854  rootProj, corrEmpty, false);
1855  accNorm->SetDirectory(0);
1856 
1857  // --- Normalize to the coverage -----------------------------------
1858  if (corrEmpty) {
1859  ScaleToCoverage(copy, accNorm);
1860  // --- Event-level normalization ---------------------------------
1861  copy->Scale(scaler);
1862  }
1863 
1864  // --- Project on X axis -------------------------------------------
1865  TH1D* dndeta = ProjectX(copy, Form("dndeta%s",base.Data()),
1866  1, nY, rootProj, corrEmpty);
1867  dndeta->SetDirectory(0);
1868  // Event-level normalization
1869  if (!corrEmpty) {
1870  ScaleToCoverage(dndeta, accNorm);
1871  dndeta->Scale(scaler);
1872  }
1873  dndeta->Scale(1., "width");
1874  copy->Scale(1., "width");
1875 
1876  TH1D* dndetaMCCorrection = 0;
1877  TH1D* dndetaMCtruth = 0;
1878  TList* centlist = 0;
1879  TList* truthcentlist = 0;
1880 
1881  // --- Possible final correction to <MC analysis> / <MC truth> -----
1882  // we get the rebinned distribution for satellite to make the correction
1883  TString rebinSuf(fSatelliteVertices ? "_rebin05" : "");
1884  if(mclist) {
1885  centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
1886  if(centlist)
1887  dndetaMCCorrection =
1888  static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s",
1889  base.Data(),
1890  rebinSuf.Data())));
1891  }
1892  if (truthlist) {
1893  truthcentlist = static_cast<TList*>(truthlist->FindObject(GetListName()));
1894  if (truthcentlist)
1895  // TODO here new is "dndetaTruth"
1896  dndetaMCtruth =
1897  static_cast<TH1D*>(truthcentlist->FindObject(Form("dndetaMCTruth%s",
1898  rebinSuf.Data())));
1899  }
1900 
1901  if (dndetaMCCorrection && dndetaMCtruth) {
1902  AliInfo("Correcting with final MC correction");
1903  TString testString(dndetaMCCorrection->GetName());
1904 
1905  // We take the measured MC dN/deta and divide with truth
1906  dndetaMCCorrection->Divide(dndetaMCtruth);
1907  dndetaMCCorrection->SetTitle("Final MC correction");
1908  dndetaMCCorrection->SetName("finalMCCorr");
1909  for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) {
1910  if(dndetaMCCorrection->GetBinContent(m) < 0.5 ||
1911  dndetaMCCorrection->GetBinContent(m) > 1.75) {
1912  dndetaMCCorrection->SetBinContent(m,1.);
1913  dndetaMCCorrection->SetBinError(m,0.1);
1914  }
1915  }
1916  // Applying the correction
1917  if (!fSatelliteVertices)
1918  // For non-satellites we took the same binning, so we do a straight
1919  // division
1920  dndeta->Divide(dndetaMCCorrection);
1921  else {
1922  // For satelitte events, we took coarser binned histograms, so
1923  // we need to do a bit more
1924  for(Int_t m = 1; m <= dndeta->GetNbinsX(); m++) {
1925  if(dndeta->GetBinContent(m) <= 0.01 ) continue;
1926 
1927  Double_t eta = dndeta->GetXaxis()->GetBinCenter(m);
1928  Int_t bin = dndetaMCCorrection->GetXaxis()->FindBin(eta);
1929  Double_t mccorr = dndetaMCCorrection->GetBinContent(bin);
1930  Double_t mcerror = dndetaMCCorrection->GetBinError(bin);
1931  if (mccorr < 1e-6) {
1932  dndeta->SetBinContent(m, 0);
1933  dndeta->SetBinError(m, 0);
1934  }
1935  Double_t value = dndeta->GetBinContent(m);
1936  Double_t error = dndeta->GetBinError(m);
1937  Double_t sumw2 = (error * error * mccorr * mccorr +
1938  mcerror * mcerror * value * value);
1939  dndeta->SetBinContent(m,value/mccorr) ;
1940  dndeta->SetBinError(m,TMath::Sqrt(sumw2)/mccorr/mccorr);
1941  }
1942  }
1943  }
1944  else
1945  DMSG(fDebug,1,"No final MC correction applied");
1946 
1947  // --- Set some histogram attributes -------------------------------
1948  TString post;
1949  Int_t rColor = GetColor(color);
1950  if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
1951  SetHistogramAttributes(dndeta, rColor, marker,
1952  Form("ALICE %s%s", GetName(), post.Data()));
1953  SetHistogramAttributes(accNorm, rColor, marker,
1954  Form("ALICE %s normalisation%s",
1955  GetName(), post.Data()));
1956 
1957  // --- Add the histograms to output list ---------------------------
1958  fOutput->Add(dndeta);
1959  fOutput->Add(accNorm);
1960  fOutput->Add(copy);
1961  if (dndetaMCCorrection) fOutput->Add(dndetaMCCorrection);
1962 
1963  // HHD Test of dN/deta in phi bins add flag later?
1964  //
1965  // cholm comment: We disable this for now
1966 #if 0
1967  for (Int_t nn=1; nn <= sum->GetNbinsY(); nn++) {
1968  TH1D* dndeta_phi = ProjectX(copy, Form("dndeta%s_phibin%d",
1969  base.Data(), nn),
1970  nn, nn, rootProj, corrEmpty);
1971  dndeta_phi->SetDirectory(0);
1972  // Event-level normalization
1973  dndeta_phi->Scale(TMath::Pi()/10., "width");
1974 
1975  if(centlist)
1976  dndetaMCCorrection =
1977  static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s_phibin%d",
1978  base.Data(),nn)));
1979  if(truthcentlist)
1980  dndetaMCtruth
1981  = static_cast<TH1D*>(truthcentlist->FindObject("dndetaMCTruth"));
1982 
1983  if (dndetaMCCorrection && dndetaMCtruth) {
1984  AliInfo("Correcting with final MC correction");
1985  TString testString(dndetaMCCorrection->GetName());
1986  dndetaMCCorrection->Divide(dndetaMCtruth);
1987  dndetaMCCorrection->SetTitle(Form("Final_MC_correction_phibin%d",nn));
1988  dndetaMCCorrection->SetName(Form("Final_MC_correction_phibin%d",nn));
1989  for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) {
1990  if(dndetaMCCorrection->GetBinContent(m) < 0.25 ||
1991  dndetaMCCorrection->GetBinContent(m) > 1.75) {
1992  dndetaMCCorrection->SetBinContent(m,1.);
1993  dndetaMCCorrection->SetBinError(m,0.1);
1994  }
1995  }
1996  //Applying the correction
1997  dndeta_phi->Divide(dndetaMCCorrection);
1998  }
1999  fOutput->Add(dndeta_phi);
2000  if(dndetaMCCorrection) fOutput->Add(dndetaMCCorrection);
2001  } // End of phi
2002 #endif
2003 }
2004 
2005 //________________________________________________________________________
2006 void
2008  TList* results,
2009  UShort_t scheme,
2010  Double_t trigEff,
2011  Double_t trigEff0,
2012  Bool_t rootProj,
2013  Bool_t corrEmpty,
2014  Int_t triggerMask,
2015  Int_t marker,
2016  Int_t color,
2017  TList* mclist,
2018  TList* truthlist)
2019 {
2020  //
2021  // End of processing
2022  //
2023  // Parameters:
2024  // sums List of sums
2025  // results Output list of results
2026  // trigEff Trigger efficiency
2027  // corrEmpty Whether to correct for empty bins
2028  // triggerMask Trigger mask
2029  //
2030  DGUARD(fDebug,1,"End centrality bin procesing");
2031 
2032  fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
2033  if(!fSums) {
2034  AliError("Could not retrieve TList fSums");
2035  return;
2036  }
2037 
2038  fOutput = new TList;
2039  fOutput->SetName(GetListName());
2040  fOutput->SetOwner();
2041  results->Add(fOutput);
2042 
2043  if (!fSum) {
2044  if (!ReadSum(fSums, false)) {
2045  AliInfo("This task did not produce any output");
2046  return;
2047  }
2048  }
2049  if (!fSumMC) ReadSum(fSums, true);
2050 
2051  fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
2052 
2053  if (!fTriggers) {
2054  AliError("Couldn't find histogram 'triggers' in list");
2055  return;
2056  }
2057 
2058  // --- Get normalization scaler ------------------------------------
2059  Double_t epsilonT = trigEff;
2060  Double_t epsilonT0 = trigEff0;
2061  DMSG(fDebug,2,"Using epsilonT=%f, epsilonT0=%f for 0x%x",
2062  epsilonT, epsilonT0, triggerMask);
2063 
2064  // Get our histograms
2065  Double_t nSum = 0;
2066  TH2D* sum = fSum->CalcSum(fOutput, nSum, epsilonT0, epsilonT,
2067  marker, rootProj, corrEmpty);
2068  Double_t nSumMC = 0;
2069  TH2D* sumMC = 0;
2070  if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC,
2071  epsilonT0, epsilonT, marker,
2072  rootProj, corrEmpty);
2073  if (!sum) {
2074  AliError("Failed to get sum from summer - bailing out");
2075  return;
2076  }
2077 
2078  TString text;
2079  Double_t ntotal = nSum;
2080  Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
2081  if (scaler < 0) {
2082  AliError("Failed to calculate normalization - bailing out");
2083  return;
2084  }
2085  fOutput->Add(fTriggers->Clone());
2086  fOutput->Add(new TNamed("normCalc", text.Data()));
2087  fOutput->Add(new TParameter<float>("low", fLow, 'f'));
2088  fOutput->Add(new TParameter<float>("high",fHigh, 'f'));
2089 
2090  // --- Make result and store ---------------------------------------
2091  MakeResult(sum, "", rootProj, corrEmpty, scaler, marker, color,
2092  mclist, truthlist);
2093 
2094  // --- Process result from TrackRefs -------------------------------
2095  if (sumMC)
2096  MakeResult(sumMC, "MC", rootProj, corrEmpty, scaler,
2097  GetMarkerStyle(GetMarkerBits(marker)+4), color,
2098  mclist, truthlist);
2099 
2100  // Temporary stuff
2101  // if (!IsAllBin()) return;
2102 
2103 }
2104 
2105 //____________________________________________________________________
2106 void
2108 {
2109  PFV("Centrality bin", GetTitle());
2110  gROOT->IncreaseDirLevel();
2111  PF("Range", "%6.2f - %6.2f%%", fLow, fHigh);
2112  PFB("All bin", IsAllBin());
2113  PFB("Final MC", fDoFinalMCCorrection);
2114  PFB("Satellite collisions", fSatelliteVertices);
2115 
2116  TString opt(option);
2117  opt.ToUpper();
2118  if (opt.Contains("R")) {
2119  if (fSum) fSum->Print(option);
2120  if (fSumMC) fSumMC->Print(option);
2121  }
2122 
2123  gROOT->DecreaseDirLevel();
2124 }
2125 
2126 //====================================================================
2127 Bool_t
2129  TH2D* data)
2130 {
2131  if (!fEmpiricalCorrection || !data)
2132  return true;
2133  Float_t zvertex=aod->GetIpZ();
2134  Int_t binzvertex=fEmpiricalCorrection->GetXaxis()->FindBin(zvertex);
2135  if(binzvertex<1||binzvertex>fEmpiricalCorrection->GetNbinsX())
2136  return false;
2137  for (int i=1;i<=data->GetNbinsX();i++) {
2138  Int_t bincorrection=fEmpiricalCorrection->GetYaxis()
2139  ->FindBin(data->GetXaxis()->GetBinCenter(i));
2140  if(bincorrection<1||bincorrection>fEmpiricalCorrection->GetNbinsY())
2141  return false;
2142  Float_t correction=fEmpiricalCorrection
2143  ->GetBinContent(binzvertex,bincorrection);
2144  if(correction<0.001) {
2145  data->SetBinContent(i,0,0);
2146  data->SetBinContent(i,data->GetNbinsY()+1,0);
2147  }
2148  for(int j=1;j<=data->GetNbinsY();j++) {
2149  if (data->GetBinContent(i,j)>0.0) {
2150  data->SetBinContent(i,j,data->GetBinContent(i,j)*correction);
2151  data->SetBinError(i,j,data->GetBinError(i,j)*correction);
2152  }
2153  }
2154  }
2155  return true;
2156 }
2157 
2158 //____________________________________________________________________
2159 TH1*
2161 {
2162  if (!h) return 0;
2163 
2164  TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
2165  // Int_t oBins = h->GetNbinsX();
2166  // Double_t high = h->GetXaxis()->GetXmax();
2167  // Double_t low = h->GetXaxis()->GetXmin();
2168  // Double_t dBin = (high - low) / oBins;
2169  // Int_t tBins = Int_t(2*high/dBin+.5);
2170  // ret->SetBins(tBins, -high, high);
2171  ret->SetDirectory(0);
2172  ret->Reset();
2173  ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
2174  ret->SetYTitle("Right/Left");
2175  Int_t nBins = h->GetNbinsX();
2176  for (Int_t i = 1; i <= nBins; i++) {
2177  Double_t x = h->GetBinCenter(i);
2178  if (x > 0) break;
2179 
2180  Double_t c1 = h->GetBinContent(i);
2181  Double_t e1 = h->GetBinError(i);
2182  if (c1 <= 0) continue;
2183 
2184  Int_t j = h->FindBin(-x);
2185  if (j <= 0 || j > nBins) continue;
2186 
2187  Double_t c2 = h->GetBinContent(j);
2188  Double_t e2 = h->GetBinError(j);
2189 
2190  Double_t c12 = c1*c1;
2191  Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
2192 
2193  Int_t k = ret->FindBin(x);
2194  ret->SetBinContent(k, c2/c1);
2195  ret->SetBinError(k, e);
2196  }
2197 
2198  return ret;
2199 }
2200 
2201 //
2202 // EOF
2203 //
Int_t color[]
virtual Bool_t CheckEvent(const AliAODForwardMult &forward)
#define PFB(N, FLAG)
const char * GetResultName(const char *postfix="") const
return jsonbuilder str().c_str()
void Print(Option_t *option="R") const
double Double_t
Definition: External.C:58
static void ScaleToCoverage(TH2D *copy, const TH1D *norm)
void SetTriggerBits(UInt_t bits)
static const char * CenterOfMassEnergyString(UShort_t cms)
const char * title
Definition: MakeQAPdf.C:26
static UShort_t ParseNormalizationScheme(const Char_t *what)
Double_t GetCentrality(AliAODEvent &event, AliAODForwardMult *forward, Int_t &qual)
#define TRUNC(X)
virtual Bool_t Event(AliAODEvent &aod)
char Char_t
Definition: External.C:18
TList * list
static TH1I * MakeTriggerHistogram(const char *name="triggers", UInt_t mask=0)
virtual void InitializeCentBins()
#define DMSG(L, N, F,...)
static TH1 * Asymmetry(TH1 *h)
TH2D * CalcSum(TList *o, Double_t &ntotal, Double_t zeroEff, Double_t otherEff=1, Int_t marker=20, Bool_t rootXproj=false, Bool_t corrEmpty=true) const
TCanvas * c
Definition: TestFitELoss.C:172
virtual void MakeResult(const TH2D *sum, const char *postfix, bool rootProj, bool corrEmpty, Double_t scaler, Int_t marker, Int_t color, TList *mclist, TList *truthlist)
virtual void CheckEventData(Double_t vtx, TH2 *data, TH2 *mcData)
virtual void SetDebugLevel(Int_t level)
Float_t GetIpZ() const
virtual void Print(Option_t *option="") const
static Int_t GetMarkerStyle(UShort_t bits)
Int_t cH
Definition: Combine.C:26
static Int_t FlipHollowStyle(Int_t style)
virtual Bool_t Finalize()
AliAnalysisUtils fAnaUtil
Per-event per bin.
int Int_t
Definition: External.C:63
void SetNormalizationScheme(UShort_t scheme)
void Init(TList *list, const TH2D *data, Int_t col)
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
static Bool_t IsTriggerBits(UInt_t bits, UInt_t trg)
float Float_t
Definition: External.C:68
virtual Bool_t CheckEvent(const AliAODForwardMult *forward, Int_t triggerMask, Double_t vzMin, Double_t vzMax, Int_t filter)
#define PF(N, V,...)
Various utilities used in PWGLF/FORWARD.
Definition: External.C:228
Definition: External.C:212
static void GetParameter(TObject *o, UShort_t &value)
void Print(Option_t *option="") const
virtual Int_t GetColor() const
virtual void CreateSums(const TH2D *data, const TH2D *mc)
#define PFV(N, VALUE)
unsigned long ULong_t
Definition: External.C:38
virtual Bool_t ProcessEvent(const AliAODForwardMult *forward, UInt_t triggerMask, Bool_t isZero, Double_t vzMin, Double_t vzMax, const TH2D *data, const TH2D *mc, UInt_t filter, Double_t weight)
TObjArray * fListOfCentralities
static const Char_t * GetTriggerString(UInt_t mask, const char *sep="&")
#define DGUARD(L, N, F,...)
virtual Bool_t Book()
CentralityBin & operator=(const CentralityBin &other)
static Int_t GetCentMethodID(const TString &meth)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
static const char * GetCentMethod(UShort_t id)
virtual CentralityBin * MakeCentralityBin(const char *name, Float_t low, Float_t high) const
static TObject * MakeParameter(const char *name, UShort_t value)
void SetSatelliteVertices(Bool_t satVtx)
const char * fwd
virtual TH2D * GetHistogram(const AliAODEvent &aod, Bool_t mc=false)=0
static TH1I * MakeStatusHistogram(const char *name="status")
virtual Bool_t ReadSum(TList *list, bool mc=false)
void AddCentralityBin(UShort_t at, Float_t low, Float_t high)
Bool_t ApplyEmpiricalCorrection(const AliAODForwardMult *aod, TH2D *data)
TH1 * GetResult(const char *postfix="", Bool_t verbose=true) const
static Float_t GetCentrality(const AliVEvent &event, const TString &method, Int_t &qual, Bool_t verbose=false)
Definition: External.C:220
static TH1D * ProjectX(const TH2D *h, const char *name, Int_t firstbin, Int_t lastbin, bool useROOT=false, bool corr=true, bool error=true)
void Add(const TH2D *data, Bool_t isZero, Double_t weight)
AliAODForwardMult * GetForward(const AliAODEvent &aod, Bool_t mc=false, Bool_t verb=true)
virtual void End(TList *sums, TList *results, UShort_t scheme, Double_t trigEff, Double_t trigEff0, Bool_t rootProj, Bool_t corrEmpty, Int_t triggerMask, Int_t marker, Int_t color, TList *mclist, TList *truthlist)
unsigned short UShort_t
Definition: External.C:28
virtual Bool_t CheckEvent(const AliAODForwardMult &forward)
const char Option_t
Definition: External.C:48
static TString GetHistName(const char *name, Int_t what=0, const char *post=0)
#define TESTAPPEND(SCHEME, BIT, STRING)
static const char * CollisionSystemString(UShort_t sys)
virtual void CreateOutputObjects(TList *dir, Int_t mask)
virtual Int_t GetMarker() const
const Int_t nbins
bool Bool_t
Definition: External.C:53
virtual Double_t GetCentrality(AliAODEvent &event, AliAODForwardMult *forward, Int_t &qual)
virtual void Print(Option_t *option="") const
Bool_t CheckEvent(UInt_t triggerMask=kInel, Double_t vzMin=-10, Double_t vzMax=10, Double_t cMin=0, Double_t cMax=100, TH1 *hist=0, TH1 *status=0, UInt_t filterMask=kDefaultFilter) const
Int_t GetColor(Int_t fallback=kRed+2) const
Definition: External.C:196
static UShort_t GetMarkerBits(Int_t style)
Bool_t HasCentrality() const
Bool_t SetCentralityMethod(const TString &method)
static const Char_t * NormalizationSchemeString(UShort_t scheme)
virtual Double_t Normalization(const TH1I &t, UShort_t scheme, Double_t trgEff, Double_t &ntotal, TString *text) const
static void SetHistogramAttributes(TH1D *h, Int_t colour, Int_t marker, const char *title, const char *ytitle=0)