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