AliPhysics  ebc57ae (ebc57ae)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FastCentHelper.C
Go to the documentation of this file.
1 #ifndef FASTCENTHELPER_H
2 #define FASTCENTHELPER_H
3 #ifndef __CINT__
4 # include <TAxis.h>
5 # include <TString.h>
6 # include <TList.h>
7 # include <TH1.h>
8 # include <TH2.h>
9 # include <TError.h>
10 # include <THStack.h>
11 # include <TClass.h>
12 # include <TStyle.h>
13 # include <TProfile.h>
14 #else
15 class TAxis;
16 class TString;
17 class TList;
18 class TH1;
19 class TH2;
20 class TCollection;
21 class TClass;
22 class TH1D;
23 class TArrayI;
24 #endif
25 
27 {
38  FastCentHelper(const char* method)
39  : fCentAxis(0),
40  fCentMeth(method),
41  fCentList(0),
42  fCentAll(0),
43  fCentAcc(0),
44  fCentMult(0),
45  fCentNPart(0),
46  fCentNBin(0),
47  fCentB(0),
48  fMapping(0)
49  {
50  TString axis("default");
51  TString meth(method);
52  TObjArray* tokens = meth.Tokenize(":");
53 
54  if (fCentMeth.IsNull()) {
55  SetCentAxis(axis);
56  return;
57  }
58  fCentMeth = tokens->At(0)->GetName();
59  if (tokens->GetEntriesFast() > 1)
60  axis = tokens->At(1)->GetName();
61 
62  if (fCentMeth.Contains("RefMult")) SetCentAxis("mult");
63  else SetCentAxis(axis);
64  }
72  void SetCentAxis(Int_t n, Double_t low, Double_t high)
73  {
74  if (fCentAxis) {
75  delete fCentAxis;
76  fCentAxis = 0;
77  }
78  fCentAxis = new TAxis(n, low, high);
79  }
86  void SetCentAxis(Int_t n, Double_t* edges)
87  {
88  if (fCentAxis) {
89  delete fCentAxis;
90  fCentAxis = 0;
91  }
92  fCentAxis = new TAxis(n, edges);
93  }
109  void SetCentAxis(const char* spec)
110  {
111  Info("SetCentAxis", "Setting centrality axis from %s", spec);
112  TString s(spec);
113  s.ToLower();
114  if (s.IsNull()) return;
115  if (s.EqualTo("pbpb") || s.EqualTo("aa") || s.EqualTo("default")) {
116  Printf("Setting centrality axis Pb-Pb");
117  Double_t aa[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
118  SetCentAxis(11, aa);
119  fCentAxis->SetUniqueID(2);
120  return;
121  }
122  if (s.EqualTo("ppb") || s.EqualTo("pbp") ||
123  s.EqualTo("pa") || s.EqualTo("ap")) {
124  Printf("Setting centrality axis p-Pb/Pb-p");
125  Double_t pa[] = { 0, 5, 10, 20, 40, 60, 80, 100 };
126  SetCentAxis(7, pa);
127  fCentAxis->SetUniqueID(3);
128  return;
129  }
130  if (s.EqualTo("pp")) {
131  Printf("Setting centrality axis pp");
132  Double_t pp[] = { 0, 0.01, 0.1, 1, 5, 10, 15, 20, 30, 40, 50, 70, 100 };
133  SetCentAxis(12, pp);
134  fCentAxis->SetUniqueID(1);
135  return;
136  }
137 
138  TObjArray* tokens = s.Tokenize(":");
139  Int_t nTokens = tokens->GetEntriesFast();
140  TArrayD edges(nTokens);
141  for (Int_t i = 0; i < nTokens; i++) {
142  TObjString* token = static_cast<TObjString*>(tokens->At(i));
143  TString& edge = token->String();
144  edges[i] = edge.Atof();
145  }
146  SetCentAxis(edges.GetSize()-1, edges.GetArray());
147  delete tokens;
148  }
157  {
158  Float_t fc = high / 100;
159  Int_t nCol = gStyle->GetNumberOfColors();
160  Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
161  Int_t col = gStyle->GetColorPalette(icol);
162  //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
163  return col;
164  }
173  virtual const char* HistName(Double_t low, Double_t high) const
174  {
175  return Form("h%03dd%02d_%03dd%02d",
176  Int_t(low), Int_t(100*low) %100,
177  Int_t(high), Int_t(100*high)%100);
178  }
179  virtual const char* HistName(Int_t i) const
180  {
181  if (!fCentAxis) return "?";
182  return HistName(fCentAxis->GetBinLowEdge(i), fCentAxis->GetBinUpEdge(i));
183  }
192  virtual const char* HistTitle(Double_t low, Double_t high) const
193  {
194  if (fCentMeth.BeginsWith("RefMult")) {
195  Int_t iLow = low;
196  Int_t iHigh = high;
197  if (iLow == iHigh) {
198  if (iLow == 0) return " 0";
199  else return Form("%2d+", iLow);
200  }
201  return Form("%2d-%2d", iLow+1, iHigh);
202  }
203  return Form("%6.2f-%6.2f%%", low, high);
204  }
213  void ModHist(TObject* o, Double_t low, Double_t high, Int_t lvl=0)
214  {
215  if (!o) return;
216  if (lvl == 0) {
217  Printf("Base level, setting names on %p", o);
218  if (o->IsA()->InheritsFrom(TNamed::Class())) {
219  TNamed* n = static_cast<TNamed*>(o);
220  n->SetName(HistName(low, high));
221  n->SetTitle(HistTitle(low, high));
222  }
223  }
224  if (o->IsA()->InheritsFrom(TCollection::Class())) {
225  Printf("Got a collection, will loop on that");
226  TCollection* c = static_cast<TCollection*>(o);
227  c->SetName(HistName(low, high));
228 
229  TObject* oo = 0;
230  TIter next(c);
231  while ((oo = next())) ModHist(oo, low, high, lvl+1);
232  return;
233  }
234  Color_t col = GetCentralityColor(low, high);
235  if (o->IsA()->InheritsFrom(TH1::Class())) {
236  Printf("Got a histogram, will modify that");
237  TH1* h = static_cast<TH1*>(o);
238  h->SetLineColor(col);
239  h->SetLineStyle(1);
240  h->SetMarkerColor(col);
241  h->SetMarkerStyle(24);
242  h->SetFillColor(kWhite);
243  h->SetFillStyle(0);
244  h->SetMinimum(0);
245  h->Sumw2();
246  }
247  }
254  void CreateHistos(TCollection* output, TObject* (*callback)())
255  {
256  // Create our stack
257  // fCentStack = new THStack("stack", "Stack of all dN_{ch}/d#eta");
258  fCentList = new TList;
259  fCentList->SetName("byCent");
260  for (Int_t i = 1; i <= fCentAxis->GetNbins(); i++) {
261  Double_t low = fCentAxis->GetBinLowEdge(i);
262  Double_t high = fCentAxis->GetBinUpEdge(i);
263  Printf("Calling call-back for %5.1f-%5.1f%%", low, high);
264  TObject* obj = (*callback)();
265  Printf("Got call back return %p", obj);
266  ModHist(obj, low, high);
267  Printf("Now adding %s (%s)to centrality list",
268  obj->GetName(), obj->ClassName());
269  fCentList->Add(obj);
270  }
271  if (fCentMeth.BeginsWith("RefMult")) {
272  // Create null-bin
273  TObject* first = (*callback)();
274  ModHist(first, 0, 0);
275  fCentList->AddFirst(first);
276 
277  // Create overflow-bin
278  TObject* last = (*callback)();
279  ModHist(last, 100, 100);
280  fCentList->AddLast(last);
281  }
282  output->Add(fCentList);
283  // output->ls();
284 
285  if (fCentAxis->GetXbins()->GetArray())
286  fCentAll = new TH1D("cent", "All Centralities",
287  fCentAxis->GetNbins(),
288  fCentAxis->GetXbins()->GetArray());
289  else
290  fCentAll = new TH1D("cent", "All Centralities",
291  fCentAxis->GetNbins(),
292  fCentAxis->GetXmin(),
293  fCentAxis->GetXmax());
294  fCentAll->SetXTitle("Centrality [%]");
295  fCentAll->SetYTitle("Events");
296  fCentAll->SetFillColor(kRed+2);
297  fCentAll->SetFillStyle(3002);
298  fCentAll->SetMarkerStyle(20);
299  fCentAll->SetMarkerColor(kRed+2);
300  fCentAll->SetLineColor(kRed+2);
301  fCentAll->SetDirectory(0);
302  fCentAll->SetMinimum(0);
303  output->Add(fCentAll);
304 
305  fCentAcc = static_cast<TH1*>(fCentAll->Clone("centAcc"));
306  fCentAcc->SetTitle("Accepted centralities");
307  fCentAcc->SetFillColor(kGreen+2);
308  fCentAcc->SetMarkerColor(kGreen+2);
309  fCentAcc->SetLineColor(kGreen+2);
310  fCentAcc->SetDirectory(0);
311  output->Add(fCentAcc);
312  // output->ls();
313 
314  Printf("End of create outputs");
315  }
322  void CreateDiagnostics(TCollection* output, TH1* centHist)
323  {
324  if (fCentMult || fCentNPart) return;
325 
326  fCentMult = 0;
327  fCentNPart = 0;
328  fCentNBin = 0;
329  Int_t maxNPart = 2*210;
330  Int_t maxNBin = 7*210;
331  if (fCentAxis->GetXbins()->GetArray()) {
332  fCentNPart = new TH2D("centNPart", "Centrality vs. N_{part}",
333  fCentAxis->GetNbins(),
334  fCentAxis->GetXbins()->GetArray(),
335  maxNPart, 0, maxNPart);
336  fCentNBin = new TH2D("centNBin", "Centrality vs. N_{bin}",
337  fCentAxis->GetNbins(),
338  fCentAxis->GetXbins()->GetArray(),
339  maxNBin, 0, maxNBin);
340  fCentB = new TH2D("centB", "Centrality vs. b",
341  fCentAxis->GetNbins(),
342  fCentAxis->GetXbins()->GetArray(),
343  200, 0, 20);
344  } else {
345  fCentNPart = new TH2D("centNPart", "Centrality vs. N_{part}",
346  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
347  fCentAxis->GetXmax(),
348  maxNPart, 0, maxNPart);
349  fCentNBin = new TH2D("centNBin", "Centrality vs. N_{bin}",
350  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
351  fCentAxis->GetXmax(),
352  maxNBin, 0, maxNBin);
353  fCentB = new TH2D("centB", "Centrality vs. b",
354  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
355  fCentAxis->GetXmax(),
356  200, 0, 20);
357  }
358 
359  fCentNPart->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
360  fCentNPart->SetYTitle("N_{part}");
361  fCentNPart->SetDirectory(0);
362 
363  fCentNBin->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
364  fCentNBin->SetYTitle("N_{bin}");
365  fCentNBin->SetDirectory(0);
366 
367  fCentB->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
368  fCentB->SetYTitle("b [fm]");
369  fCentB->SetDirectory(0);
370 
371  output->Add(fCentNPart);
372  output->Add(fCentNBin);
373  output->Add(fCentB);
374 
375  if (!centHist) return;;
376  if (fCentAxis->GetXbins()->GetArray()) {
377  fCentMult = new TH2D("centMult","Event multiplicity vs. centrality",
378  centHist->GetXaxis()->GetNbins(),
379  centHist->GetXaxis()->GetXmin(),
380  centHist->GetXaxis()->GetXmax(),
381  fCentAxis->GetNbins(),
382  fCentAxis->GetXbins()->GetArray());
383  } else {
384  fCentMult = new TH2D("centMult","Event multiplicity vs. centrality",
385  centHist->GetXaxis()->GetNbins(),
386  centHist->GetXaxis()->GetXmin(),
387  centHist->GetXaxis()->GetXmax(),
388  fCentAxis->GetNbins(),
389  fCentAxis->GetXmin(), fCentAxis->GetXmax());
390  }
391  fCentMult->SetXTitle(Form("Event multiplicity (%s)", fCentMeth.Data()));
392  fCentMult->SetYTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
393  fCentMult->SetDirectory(0);
394  output->Add(fCentMult);
395 
396  return;
397  }
410  Double_t mult,
411  Double_t b,
412  Double_t nPart,
413  Double_t nBin)
414  {
415  Int_t ret = -1;
416  fCentAll->Fill(cent);
417  if (fCentMult) fCentMult->Fill(mult, cent);
418  fCentB->Fill(cent, b);
419  if (cent < 0 || cent > 999) {
420  Warning("FillDiagnostics",
421  "Centrality is unreasonable: %f -> %f",mult, cent);
422  return ret;
423  }
424  Int_t n = fCentAxis->GetNbins();
425  ret = fCentAxis->FindBin(cent);
426  if (ret-1 == n && cent == fCentAxis->GetXmax())
427  ret = n;
428  if (ret < 1 || ret > n) {
429  Warning("ProcessHeader", "Centrality %f -> %f -> bin # %d",
430  mult, cent, ret);
431  ret = -1;
432  return ret;
433  }
434  fCentNPart->Fill(cent, nPart);
435  fCentNBin->Fill(cent, nBin);
436  if (ret == n) cent -= 0.001;
437  fCentAcc->Fill(cent);
438  return ret;
439  }
450  TObject* CentObject(Int_t bin, TClass* cls=0) const
451  {
452  if (!fCentList) {
453  Warning("CentObject", "No centrality objects defined");
454  return 0;
455  }
456  if (fMapping.GetArray()) {
457  if (bin >= fMapping.GetSize()) {
458  Warning("CentObject", "Bin # %2d out of range [%2d,%2d]",
459  bin, 1, fMapping.GetSize());
460  return 0;
461  }
462  Int_t old = bin;
463  bin = fMapping[old];
464  // Info("CentObject", "Mapped bin %d to index %d", old, bin);
465  }
466  if (bin <= 0 || bin > fCentList->GetEntries()) {
467  if (!fMapping.GetArray())
468  Warning("CentObject", "Bin # %2d out of range [%2d,%2d]",
469  bin, 1, fCentList->GetEntries());
470  return 0;
471  }
472  TObject* o = fCentList->At(bin-1);
473  if (!o) {
474  Warning("CentObject", "No centrality object defined for bin %d", bin);
475  return 0;
476  }
477  if (cls && !o->IsA()->InheritsFrom(cls)) {
478  Warning("CentObject", "Centrality object %s is not a %s, but a %s",
479  o->GetName(), cls->GetName(), o->ClassName());
480  return 0;
481  }
482  return o;
483  }
493  TH1* CentHist(Int_t bin) const
494  {
495  return static_cast<TH1*>(CentObject(bin,TH1::Class()));
496  }
507  {
508  return static_cast<TCollection*>(CentObject(bin,TCollection::Class()));
509  }
510 
520  TObject* GetOutputObject(TCollection* output, const char* name, TClass* cls)
521  {
522  TObject* o = output->FindObject(name);
523  if (!o) {
524  Warning("GetOutputObject", "Object %s not found in output", name);
525  output->ls();
526  return 0;
527  }
528  if (cls && !o->IsA()->InheritsFrom(cls)) {
529  Warning("GetOutputObject", "Object %s from output is not a %s, but a %s",
530  o->GetName(), cls->GetName(), o->ClassName());
531  return o;
532  }
533  return o;
534  }
536  const char* name,
537  Bool_t projX)
538  {
539  TH2* h = static_cast<TH2*>(GetOutputObject(output, name, TH2::Class()));
540  if (!h) return 0;
541  TProfile* p = (projX ?
542  h->ProfileX(Form("%sMean", name)) :
543  h->ProfileY(Form("%sMean", name)));
544  p->SetDirectory(0);
545  output->Add(p);
546 
547  return h;
548  }
549  Bool_t Finalize(TCollection* output, Long_t minEvents,
550  TH1* (*callback)(TObject*,Int_t))
551  {
552 
553  fCentAll = static_cast<TH1*>(GetOutputObject(output, "cent",
554  TH1::Class()));
555  fCentAcc = static_cast<TH1*>(GetOutputObject(output, "centAcc",
556  TH1::Class()));
557  fCentMult = Get2DDiag(output, "centMult", false);
558  fCentNPart = Get2DDiag(output, "centNPart",true);
559  fCentNBin = Get2DDiag(output, "centNBin", true);
560  fCentB = Get2DDiag(output, "centB", true);
561  fCentList = static_cast<TList*>(GetOutputObject(output, "byCent",
562  TList::Class()));
563  if (!fCentList || !fCentAll || !fCentAcc) {
564  Warning("Finalize", "Missing stack and histograms");
565  return false;
566  }
567  Info("Terminate", "Accepted %d/%d=%6.2f%% events",
568  int(fCentAcc->GetEntries()), int(fCentAll->GetEntries()),
569  100*fCentAcc->GetEntries()/fCentAll->GetEntries());
570 
571  fMapping.Set(fCentAll->GetNbinsX()+1);
572  fMapping.Reset(-1);
573  THStack* stack = new THStack("all", "All");
574  TList* hists = fCentList; // ->GetHists();
575  TObjLink* link = hists->FirstLink();
576  Int_t bin = 1;
577  Int_t cnt = 0;
578  Long64_t sum = 0;
579  Long64_t total = 0;
580  Long64_t all = fCentAll->GetEntries();
581  while (link) {
582  TObject* o = link->GetObject();
583  if (!o) {
584  link = link->Next();
585  bin++;
586  continue;
587  }
588  Int_t m = fCentAcc->GetBinContent(bin);
589  total += m;
590  printf("%9d events in bin %s ...", m, o->GetTitle());
591  if (m < minEvents) {
592  // Too few event, remove this
593  TObjLink* tmp = link->Next();
594  hists->Remove(link);
595  link = tmp;
596  delete o;
597  Printf(" removed");
598  bin++;
599  continue;
600  }
601  sum += m;
602  TH1* h = (*callback)(o,m);
603  stack->Add(h);
604  Printf(" processed [%2d->%2d]", bin, cnt);
605  link = link->Next();
606  fMapping[bin] = cnt+1;
607  cnt++;
608  bin++;
609 #if 0
610  TH1* h = static_cast<TH1*>(o);
611  Int_t n = h->GetBinContent(0);
612  h->SetBinContent(0,0);
613  printf("%9d (%9d) events in bin %s ...", n, m, o->GetTitle());
614  if (n < minEvents) {
615  // Too few event, remove this
616  TObjLink* tmp = link->Next();
617  hists->Remove(link);
618  link = tmp;
619  delete o;
620  Printf(" removed");
621  bin++;
622  continue;
623  }
624  sum += m;
625  // Scale
626  h->Scale(1. / n, "width");
627  stack->Add(h);
628  Printf(" scaled");
629  link = link->Next();
630  bin++;
631 #endif
632  }
633  Printf("ana/acc/all: %9lld/%9lld/%9lld [%6.2f%%/%6.2f%%]",
634  sum, total, all, float(100*sum)/total, float(100*total)/all);
635  output->Add(stack);
636  for (Int_t i = 1; i < fMapping.GetSize(); i++)
637  Printf(" %3d -> %3d", i, fMapping[i]);
638  }
639  void Print(Option_t* option="") const
640  {
641  Int_t nBins = (fCentAxis?fCentAxis->GetNbins():0);
642  Printf(" Method: %s", fCentMeth.Data());
643  Printf(" Axis: %d", nBins);
644  if (nBins < 2) return;
645  printf(" ");
646  for (Int_t i = 1; i <= nBins; i++)
647  printf("%5.1f-",fCentAxis->GetBinLowEdge(i));
648  Printf("%5.1f%%", fCentAxis->GetBinUpEdge(nBins));
649  }
651 };
652 #endif
653 
654 //
655 // EOF
656 //
void SetCentAxis(const char *spec)
double Double_t
Definition: External.C:58
TObject * GetOutputObject(TCollection *output, const char *name, TClass *cls)
long long Long64_t
Definition: External.C:43
Bool_t Finalize(TCollection *output, Long_t minEvents, TH1 *(*callback)(TObject *, Int_t))
void SetCentAxis(Int_t n, Double_t *edges)
TH1 * CentHist(Int_t bin) const
TH2 * Get2DDiag(TCollection *output, const char *name, Bool_t projX)
virtual const char * HistName(Double_t low, Double_t high) const
virtual const char * HistTitle(Double_t low, Double_t high) const
TCanvas * c
Definition: TestFitELoss.C:172
AliStack * stack
virtual const char * HistName(Int_t i) const
void SetCentAxis(Int_t n, Double_t low, Double_t high)
int Int_t
Definition: External.C:63
TObject * CentObject(Int_t bin, TClass *cls=0) const
FastCentHelper(const char *method)
float Float_t
Definition: External.C:68
Int_t GetCentralityColor(Double_t, Double_t high) const
Int_t method
Definition: External.C:228
Definition: External.C:212
void CreateHistos(TCollection *output, TObject *(*callback)())
ClassDef(FastCentHelper, 1)
void ModHist(TObject *o, Double_t low, Double_t high, Int_t lvl=0)
TCollection * CentCollection(Int_t bin) const
void Print(Option_t *option="") const
void CreateDiagnostics(TCollection *output, TH1 *centHist)
Definition: External.C:220
Int_t CheckCentrality(Double_t cent, Double_t mult, Double_t b, Double_t nPart, Double_t nBin)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Definition: External.C:196