AliPhysics  64f4410 (64f4410)
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  if (callback) {
261  for (Int_t i = 1; i <= fCentAxis->GetNbins(); i++) {
262  Double_t low = fCentAxis->GetBinLowEdge(i);
263  Double_t high = fCentAxis->GetBinUpEdge(i);
264  Printf("Calling call-back for %5.1f-%5.1f%%", low, high);
265  TObject* obj = (*callback)();
266  Printf("Got call back return %p", obj);
267  ModHist(obj, low, high);
268  Printf("Now adding %s (%s)to centrality list",
269  obj->GetName(), obj->ClassName());
270  fCentList->Add(obj);
271  }
272  }
273  if (fCentMeth.BeginsWith("RefMult")) {
274  // Create null-bin
275  TObject* first = (*callback)();
276  ModHist(first, 0, 0);
277  fCentList->AddFirst(first);
278 
279  // Create overflow-bin
280  TObject* last = (*callback)();
281  ModHist(last, 100, 100);
282  fCentList->AddLast(last);
283  }
284  output->Add(fCentList);
285  // output->ls();
286 
287  if (fCentAxis->GetXbins()->GetArray())
288  fCentAll = new TH1D("cent", "All Centralities",
289  fCentAxis->GetNbins(),
290  fCentAxis->GetXbins()->GetArray());
291  else
292  fCentAll = new TH1D("cent", "All Centralities",
293  fCentAxis->GetNbins(),
294  fCentAxis->GetXmin(),
295  fCentAxis->GetXmax());
296  fCentAll->SetXTitle("Centrality [%]");
297  fCentAll->SetYTitle("Events");
298  fCentAll->SetFillColor(kRed+2);
299  fCentAll->SetFillStyle(3002);
300  fCentAll->SetMarkerStyle(20);
301  fCentAll->SetMarkerColor(kRed+2);
302  fCentAll->SetLineColor(kRed+2);
303  fCentAll->SetDirectory(0);
304  fCentAll->SetMinimum(0);
305  output->Add(fCentAll);
306 
307  fCentAcc = static_cast<TH1*>(fCentAll->Clone("centAcc"));
308  fCentAcc->SetTitle("Accepted centralities");
309  fCentAcc->SetFillColor(kGreen+2);
310  fCentAcc->SetMarkerColor(kGreen+2);
311  fCentAcc->SetLineColor(kGreen+2);
312  fCentAcc->SetDirectory(0);
313  output->Add(fCentAcc);
314  // output->ls();
315 
316  Printf("End of create outputs");
317  }
324  void CreateDiagnostics(TCollection* output, TH1* centHist)
325  {
326  if (fCentMult || fCentNPart) return;
327 
328  fCentMult = 0;
329  fCentNPart = 0;
330  fCentNBin = 0;
331  Int_t maxNPart = 2*210;
332  Int_t maxNBin = 7*210;
333  if (fCentAxis->GetXbins()->GetArray()) {
334  fCentNPart = new TH2D("centNPart", "Centrality vs. N_{part}",
335  fCentAxis->GetNbins(),
336  fCentAxis->GetXbins()->GetArray(),
337  maxNPart, 0, maxNPart);
338  fCentNBin = new TH2D("centNBin", "Centrality vs. N_{bin}",
339  fCentAxis->GetNbins(),
340  fCentAxis->GetXbins()->GetArray(),
341  maxNBin, 0, maxNBin);
342  fCentB = new TH2D("centB", "Centrality vs. b",
343  fCentAxis->GetNbins(),
344  fCentAxis->GetXbins()->GetArray(),
345  200, 0, 20);
346  } else {
347  fCentNPart = new TH2D("centNPart", "Centrality vs. N_{part}",
348  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
349  fCentAxis->GetXmax(),
350  maxNPart, 0, maxNPart);
351  fCentNBin = new TH2D("centNBin", "Centrality vs. N_{bin}",
352  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
353  fCentAxis->GetXmax(),
354  maxNBin, 0, maxNBin);
355  fCentB = new TH2D("centB", "Centrality vs. b",
356  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
357  fCentAxis->GetXmax(),
358  200, 0, 20);
359  }
360 
361  fCentNPart->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
362  fCentNPart->SetYTitle("N_{part}");
363  fCentNPart->SetDirectory(0);
364 
365  fCentNBin->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
366  fCentNBin->SetYTitle("N_{bin}");
367  fCentNBin->SetDirectory(0);
368 
369  fCentB->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
370  fCentB->SetYTitle("b [fm]");
371  fCentB->SetDirectory(0);
372 
373  output->Add(fCentNPart);
374  output->Add(fCentNBin);
375  output->Add(fCentB);
376 
377  if (!centHist) return;;
378  if (fCentAxis->GetXbins()->GetArray()) {
379  fCentMult = new TH2D("centMult","Event multiplicity vs. centrality",
380  centHist->GetXaxis()->GetNbins(),
381  centHist->GetXaxis()->GetXmin(),
382  centHist->GetXaxis()->GetXmax(),
383  fCentAxis->GetNbins(),
384  fCentAxis->GetXbins()->GetArray());
385  } else {
386  fCentMult = new TH2D("centMult","Event multiplicity vs. centrality",
387  centHist->GetXaxis()->GetNbins(),
388  centHist->GetXaxis()->GetXmin(),
389  centHist->GetXaxis()->GetXmax(),
390  fCentAxis->GetNbins(),
391  fCentAxis->GetXmin(), fCentAxis->GetXmax());
392  }
393  fCentMult->SetXTitle(Form("Event multiplicity (%s)", fCentMeth.Data()));
394  fCentMult->SetYTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
395  fCentMult->SetDirectory(0);
396  output->Add(fCentMult);
397 
398  return;
399  }
412  Double_t mult,
413  Double_t b,
414  Double_t nPart,
415  Double_t nBin)
416  {
417  Int_t ret = -1;
418  fCentAll->Fill(cent);
419  if (fCentMult) fCentMult->Fill(mult, cent);
420  fCentB->Fill(cent, b);
421  if (cent < 0 || cent > 999) {
422  Warning("FillDiagnostics",
423  "Centrality is unreasonable: %f -> %f",mult, cent);
424  return ret;
425  }
426  Int_t n = fCentAxis->GetNbins();
427  ret = fCentAxis->FindBin(cent);
428  if (ret-1 == n && cent == fCentAxis->GetXmax())
429  ret = n;
430  if (ret < 1 || ret > n) {
431  Warning("ProcessHeader", "Centrality %f -> %f -> bin # %d",
432  mult, cent, ret);
433  ret = -1;
434  return ret;
435  }
436  fCentNPart->Fill(cent, nPart);
437  fCentNBin->Fill(cent, nBin);
438  if (ret == n) cent -= 0.001;
439  fCentAcc->Fill(cent);
440  return ret;
441  }
452  TObject* CentObject(Int_t bin, TClass* cls=0) const
453  {
454  if (!fCentList) {
455  Warning("CentObject", "No centrality objects defined");
456  return 0;
457  }
458  if (fMapping.GetArray()) {
459  if (bin >= fMapping.GetSize()) {
460  Warning("CentObject", "Bin # %2d out of range [%2d,%2d]",
461  bin, 1, fMapping.GetSize());
462  return 0;
463  }
464  Int_t old = bin;
465  bin = fMapping[old];
466  // Info("CentObject", "Mapped bin %d to index %d", old, bin);
467  }
468  if (bin <= 0 || bin > fCentList->GetEntries()) {
469  if (!fMapping.GetArray())
470  Warning("CentObject", "Bin # %2d out of range [%2d,%2d]",
471  bin, 1, fCentList->GetEntries());
472  return 0;
473  }
474  TObject* o = fCentList->At(bin-1);
475  if (!o) {
476  Warning("CentObject", "No centrality object defined for bin %d", bin);
477  return 0;
478  }
479  if (cls && !o->IsA()->InheritsFrom(cls)) {
480  Warning("CentObject", "Centrality object %s is not a %s, but a %s",
481  o->GetName(), cls->GetName(), o->ClassName());
482  return 0;
483  }
484  return o;
485  }
495  TH1* CentHist(Int_t bin) const
496  {
497  return static_cast<TH1*>(CentObject(bin,TH1::Class()));
498  }
509  {
510  return static_cast<TCollection*>(CentObject(bin,TCollection::Class()));
511  }
512 
522  TObject* GetOutputObject(TCollection* output, const char* name, TClass* cls)
523  {
524  TObject* o = output->FindObject(name);
525  if (!o) {
526  Warning("GetOutputObject", "Object %s not found in output", name);
527  output->ls();
528  return 0;
529  }
530  if (cls && !o->IsA()->InheritsFrom(cls)) {
531  Warning("GetOutputObject", "Object %s from output is not a %s, but a %s",
532  o->GetName(), cls->GetName(), o->ClassName());
533  return o;
534  }
535  return o;
536  }
538  const char* name,
539  Bool_t projX)
540  {
541  TH2* h = static_cast<TH2*>(GetOutputObject(output, name, TH2::Class()));
542  if (!h) return 0;
543  TProfile* p = (projX ?
544  h->ProfileX(Form("%sMean", name)) :
545  h->ProfileY(Form("%sMean", name)));
546  p->SetDirectory(0);
547  output->Add(p);
548 
549  return h;
550  }
551  Bool_t Finalize(TCollection* output, Long_t minEvents,
552  TH1* (*callback)(TObject*,Int_t))
553  {
554 
555  fCentAll = static_cast<TH1*>(GetOutputObject(output, "cent",
556  TH1::Class()));
557  fCentAcc = static_cast<TH1*>(GetOutputObject(output, "centAcc",
558  TH1::Class()));
559  fCentMult = Get2DDiag(output, "centMult", false);
560  fCentNPart = Get2DDiag(output, "centNPart",true);
561  fCentNBin = Get2DDiag(output, "centNBin", true);
562  fCentB = Get2DDiag(output, "centB", true);
563  fCentList = static_cast<TList*>(GetOutputObject(output, "byCent",
564  TList::Class()));
565  if (!fCentList || !fCentAll || !fCentAcc) {
566  Warning("Finalize", "Missing stack and histograms");
567  return false;
568  }
569  Info("Terminate", "Accepted %d/%d=%6.2f%% events",
570  int(fCentAcc->GetEntries()), int(fCentAll->GetEntries()),
571  100*fCentAcc->GetEntries()/fCentAll->GetEntries());
572 
573  if (!callback) return true;
574 
575  fMapping.Set(fCentAll->GetNbinsX()+1);
576  fMapping.Reset(-1);
577  THStack* stack = new THStack("all", "All");
578  TList* hists = fCentList; // ->GetHists();
579  TObjLink* link = hists->FirstLink();
580  Int_t bin = 1;
581  Int_t cnt = 0;
582  Long64_t sum = 0;
583  Long64_t total = 0;
584  Long64_t all = fCentAll->GetEntries();
585  while (link) {
586  TObject* o = link->GetObject();
587  if (!o) {
588  link = link->Next();
589  bin++;
590  continue;
591  }
592  Int_t m = fCentAcc->GetBinContent(bin);
593  total += m;
594  printf("%9d events in bin %s ...", m, o->GetTitle());
595  if (m < minEvents) {
596  // Too few event, remove this
597  TObjLink* tmp = link->Next();
598  hists->Remove(link);
599  link = tmp;
600  delete o;
601  Printf(" removed");
602  bin++;
603  continue;
604  }
605  sum += m;
606  TH1* h = (*callback)(o,m);
607  stack->Add(h);
608  Printf(" processed [%2d->%2d]", bin, cnt);
609  link = link->Next();
610  fMapping[bin] = cnt+1;
611  cnt++;
612  bin++;
613 #if 0
614  TH1* h = static_cast<TH1*>(o);
615  Int_t n = h->GetBinContent(0);
616  h->SetBinContent(0,0);
617  printf("%9d (%9d) events in bin %s ...", n, m, o->GetTitle());
618  if (n < minEvents) {
619  // Too few event, remove this
620  TObjLink* tmp = link->Next();
621  hists->Remove(link);
622  link = tmp;
623  delete o;
624  Printf(" removed");
625  bin++;
626  continue;
627  }
628  sum += m;
629  // Scale
630  h->Scale(1. / n, "width");
631  stack->Add(h);
632  Printf(" scaled");
633  link = link->Next();
634  bin++;
635 #endif
636  }
637  output->Add(stack);
638  Printf("ana/acc/all: %9lld/%9lld/%9lld [%6.2f%%/%6.2f%%]",
639  sum, total, all, float(100*sum)/total, float(100*total)/all);
640  for (Int_t i = 1; i < fMapping.GetSize(); i++)
641  Printf(" %3d -> %3d", i, fMapping[i]);
642  }
643  void Print(Option_t* option="") const
644  {
645  Int_t nBins = (fCentAxis?fCentAxis->GetNbins():0);
646  Printf(" Method: %s", fCentMeth.Data());
647  Printf(" Axis: %d", nBins);
648  if (nBins < 2) return;
649  printf(" ");
650  for (Int_t i = 1; i <= nBins; i++)
651  printf("%5.1f-",fCentAxis->GetBinLowEdge(i));
652  Printf("%5.1f%%", fCentAxis->GetBinUpEdge(nBins));
653  }
654  ClassDef(FastCentHelper,1);
655 };
656 #endif
657 
658 //
659 // EOF
660 //
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)())
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