AliPhysics  db95e02 (db95e02)
 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 #else
14 class TAxis;
15 class TString;
16 class TList;
17 class TH1;
18 class TH2;
19 class TCollection;
20 class TClass;
21 class TH1D;
22 #endif
23 
25 {
35  FastCentHelper(const char* method)
36  : fCentAxis(0),
37  fCentMeth(""),
38  fCentList(0),
39  fCentAll(0),
40  fCentAcc(0),
41  fCentMult(0),
42  fCentNPart(0),
43  fCentNBin(0),
44  fCentB(0)
45  {
46  TString axis("default");
47  TString meth(method);
48  TObjArray* tokens = meth.Tokenize(":");
49 
50  if (fCentMeth.IsNull()) {
51  SetCentAxis(axis);
52  return;
53  }
54  fCentMeth = tokens->At(0)->GetName();
55  if (tokens->GetEntriesFast() > 1)
56  axis = tokens->At(1)->GetName();
57 
58  if (fCentMeth.Contains("RefMult")) SetCentAxis("mult");
59  else SetCentAxis(axis);
60  }
68  void SetCentAxis(Int_t n, Double_t low, Double_t high)
69  {
70  if (fCentAxis) {
71  delete fCentAxis;
72  fCentAxis = 0;
73  }
74  fCentAxis = new TAxis(n, low, high);
75  }
82  void SetCentAxis(Int_t n, Double_t* edges)
83  {
84  if (fCentAxis) {
85  delete fCentAxis;
86  fCentAxis = 0;
87  }
88  fCentAxis = new TAxis(n, edges);
89  }
105  void SetCentAxis(const char* spec)
106  {
107  Info("SetCentAxis", "Setting centrality axis from %s", spec);
108  TString s(spec);
109  s.ToLower();
110  if (s.IsNull()) return;
111  if (s.EqualTo("pbpb") || s.EqualTo("aa") || s.EqualTo("default")) {
112  Printf("Setting centrality axis Pb-Pb");
113  Double_t aa[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
114  SetCentAxis(11, aa);
115  fCentAxis->SetUniqueID(2);
116  return;
117  }
118  if (s.EqualTo("ppb") || s.EqualTo("pbp") ||
119  s.EqualTo("pa") || s.EqualTo("ap")) {
120  Printf("Setting centrality axis p-Pb/Pb-p");
121  Double_t pa[] = { 0, 5, 10, 20, 40, 60, 80, 100 };
122  SetCentAxis(7, pa);
123  fCentAxis->SetUniqueID(3);
124  return;
125  }
126  if (s.EqualTo("pp")) {
127  Printf("Setting centrality axis pp");
128  Double_t pp[] = { 0, 0.01, 0.1, 1, 5, 10, 15, 20, 30, 40, 50, 70, 100 };
129  SetCentAxis(12, pp);
130  fCentAxis->SetUniqueID(1);
131  return;
132  }
133 
134  TObjArray* tokens = s.Tokenize(":");
135  Int_t nTokens = tokens->GetEntriesFast();
136  TArrayD edges(nTokens);
137  for (Int_t i = 0; i < nTokens; i++) {
138  TObjString* token = static_cast<TObjString*>(tokens->At(i));
139  TString& edge = token->String();
140  edges[i] = edge.Atof();
141  }
142  SetCentAxis(edges.GetSize()-1, edges.GetArray());
143  delete tokens;
144  }
153  {
154  Float_t fc = high / 100;
155  Int_t nCol = gStyle->GetNumberOfColors();
156  Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5));
157  Int_t col = gStyle->GetColorPalette(icol);
158  //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
159  return col;
160  }
169  virtual const char* HistName(Double_t low, Double_t high) const
170  {
171  return Form("h%03dd%02d_%03dd%02d",
172  Int_t(low), Int_t(100*low) %100,
173  Int_t(high), Int_t(100*high)%100);
174  }
183  virtual const char* HistTitle(Double_t low, Double_t high) const
184  {
185  if (fCentMeth.BeginsWith("RefMult")) {
186  Int_t iLow = low;
187  Int_t iHigh = high;
188  if (iLow == iHigh) {
189  if (iLow == 0) return " 0";
190  else return Form("%2d+", iLow);
191  }
192  return Form("%2d-%2d", iLow+1, iHigh);
193  }
194  return Form("%6.2f-%6.2f%%", low, high);
195  }
204  void ModHist(TH1* h, Double_t low, Double_t high)
205  {
206  Color_t col = GetCentralityColor(low, high);
207  h->SetLineColor(col);
208  h->SetLineStyle(1);
209  h->SetMarkerColor(col);
210  h->SetMarkerStyle(24);
211  h->SetFillColor(kWhite);
212  h->SetFillStyle(0);
213  h->SetName(HistName(low, high));
214  h->SetTitle(HistTitle(low, high));
215  }
222  void CreateHistos(TCollection* output, TH1D* (*callback)())
223  {
224  // Create our stack
225  // fCentStack = new THStack("stack", "Stack of all dN_{ch}/d#eta");
226  fCentList = new TList;
227  fCentList->SetName("byCent");
228  for (Int_t i = 1; i <= fCentAxis->GetNbins(); i++) {
229  Double_t low = fCentAxis->GetBinLowEdge(i);
230  Double_t high = fCentAxis->GetBinUpEdge(i);
231  TH1* hist = (*callback)();
232  ModHist(hist, low, high);
233  hist->SetMinimum(0);
234  hist->Sumw2();
235  fCentList->Add(hist);
236  }
237  if (fCentMeth.BeginsWith("RefMult")) {
238  // Create null-bin
239  TH1* first = (*callback)();
240  ModHist(first, 0, 0);
241  fCentList->AddFirst(first);
242 
243  // Create overflow-bin
244  TH1* last = (*callback)();
245  ModHist(last, 100, 100);
246  fCentList->AddLast(last);
247  }
248  output->Add(fCentList);
249 
250  if (fCentAxis->GetXbins()->GetArray())
251  fCentAll = new TH1D("cent", "All Centralities",
252  fCentAxis->GetNbins(),
253  fCentAxis->GetXbins()->GetArray());
254  else
255  fCentAll = new TH1D("cent", "All Centralities",
256  fCentAxis->GetNbins(),
257  fCentAxis->GetXmin(),
258  fCentAxis->GetXmax());
259  fCentAll->SetXTitle("Centrality [%]");
260  fCentAll->SetYTitle("Events");
261  fCentAll->SetFillColor(kRed+2);
262  fCentAll->SetFillStyle(3002);
263  fCentAll->SetMarkerStyle(20);
264  fCentAll->SetMarkerColor(kRed+2);
265  fCentAll->SetLineColor(kRed+2);
266  fCentAll->SetDirectory(0);
267  fCentAll->SetMinimum(0);
268  output->Add(fCentAll);
269 
270  fCentAcc = static_cast<TH1*>(fCentAll->Clone("centAcc"));
271  fCentAcc->SetTitle("Accepted centralities");
272  fCentAcc->SetFillColor(kGreen+2);
273  fCentAcc->SetMarkerColor(kGreen+2);
274  fCentAcc->SetLineColor(kGreen+2);
275  fCentAcc->SetDirectory(0);
276  output->Add(fCentAcc);
277  // output->ls();
278  }
285  void CreateDiagnostics(TCollection* output, TH1* centHist)
286  {
287  if (fCentMult || fCentNPart) return;
288 
289  fCentMult = 0;
290  fCentNPart = 0;
291  fCentNBin = 0;
292  Int_t maxNPart = 2*210;
293  Int_t maxNBin = 3*210;
294  if (fCentAxis->GetXbins()->GetArray()) {
295  fCentNPart = new TH2D("centNPart", "Centrality vs. N_{part}",
296  fCentAxis->GetNbins(),
297  fCentAxis->GetXbins()->GetArray(),
298  maxNPart, 0, maxNPart);
299  fCentNBin = new TH2D("centNBin", "Centrality vs. N_{bin}",
300  fCentAxis->GetNbins(),
301  fCentAxis->GetXbins()->GetArray(),
302  maxNBin, 0, maxNBin);
303  fCentB = new TH2D("centB", "Centrality vs. b",
304  fCentAxis->GetNbins(),
305  fCentAxis->GetXbins()->GetArray(),
306  200, 0, 20);
307  } else {
308  fCentNPart = new TH2D("centNPart", "Centrality vs. N_{part}",
309  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
310  fCentAxis->GetXmax(),
311  maxNPart, 0, maxNPart);
312  fCentNBin = new TH2D("centNBin", "Centrality vs. N_{bin}",
313  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
314  fCentAxis->GetXmax(),
315  maxNBin, 0, maxNBin);
316  fCentB = new TH2D("centB", "Centrality vs. b",
317  fCentAxis->GetNbins(), fCentAxis->GetXmin(),
318  fCentAxis->GetXmax(),
319  200, 0, 20);
320  }
321 
322  fCentNPart->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
323  fCentNPart->SetYTitle("N_{part}");
324  fCentNPart->SetDirectory(0);
325 
326  fCentNBin->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
327  fCentNBin->SetYTitle("N_{bin}");
328  fCentNBin->SetDirectory(0);
329 
330  fCentB->SetXTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
331  fCentB->SetYTitle("b [fm]");
332  fCentB->SetDirectory(0);
333 
334  output->Add(fCentNPart);
335  output->Add(fCentNBin);
336  output->Add(fCentB);
337 
338  if (!centHist) return;;
339  if (fCentAxis->GetXbins()->GetArray()) {
340  fCentMult = new TH2D("centMult","Event multiplicity vs. centrality",
341  centHist->GetXaxis()->GetNbins(),
342  centHist->GetXaxis()->GetXmin(),
343  centHist->GetXaxis()->GetXmax(),
344  fCentAxis->GetNbins(),
345  fCentAxis->GetXbins()->GetArray());
346  } else {
347  fCentMult = new TH2D("centMult","Event multiplicity vs. centrality",
348  centHist->GetXaxis()->GetNbins(),
349  centHist->GetXaxis()->GetXmin(),
350  centHist->GetXaxis()->GetXmax(),
351  fCentAxis->GetNbins(),
352  fCentAxis->GetXmin(), fCentAxis->GetXmax());
353  }
354  fCentMult->SetXTitle(Form("Event multiplicity (%s)", fCentMeth.Data()));
355  fCentMult->SetYTitle(Form("Centrality (%s) %%", fCentMeth.Data()));
356  fCentMult->SetDirectory(0);
357  output->Add(fCentMult);
358 
359  return;
360  }
373  Double_t mult,
374  Double_t b,
375  Double_t nPart,
376  Double_t nBin)
377  {
378  Int_t ret = -1;
379  fCentAll->Fill(cent);
380  if (fCentMult) fCentMult->Fill(mult, cent);
381  fCentB->Fill(cent, b);
382  if (cent < 0 || cent > 999) {
383  Warning("FillDiagnostics",
384  "Centrality is unreasonable: %f -> %f",mult, cent);
385  return ret;
386  }
387  Int_t n = fCentAxis->GetNbins();
388  ret = fCentAxis->FindBin(cent);
389  if (ret-1 == n && cent == fCentAxis->GetXmax())
390  ret = n;
391  if (ret < 1 || ret > n) {
392  Warning("ProcessHeader", "Centrality %f -> %f -> bin # %d",
393  mult, cent, ret);
394  ret = -1;
395  return ret;
396  }
397  fCentNPart->Fill(cent, nPart);
398  fCentNBin->Fill(cent, nBin);
399  if (ret == n) cent -= 0.001;
400  fCentAcc->Fill(cent);
401  return ret;
402  }
412  TObject* GetOutputObject(TCollection* output, const char* name, TClass* cls)
413  {
414  TObject* o = output->FindObject(name);
415  if (!o) {
416  Warning("GetOutputObject", "Object %s not found in output", name);
417  output->ls();
418  return 0;
419  }
420  if (cls && !o->IsA()->InheritsFrom(cls)) {
421  Warning("GetOutputObject", "Object %s from output is not a %s, but a %s",
422  o->GetName(), cls->GetName(), o->ClassName());
423  return o;
424  }
425  return o;
426  }
427  Bool_t Finalize(TCollection* output, Long_t minEvents)
428  {
429 
430  fCentAll = static_cast<TH1*>(GetOutputObject(output, "cent",
431  TH1::Class()));
432  fCentAcc = static_cast<TH1*>(GetOutputObject(output, "centAcc",
433  TH1::Class()));
434  fCentMult = static_cast<TH2*>(GetOutputObject(output, "centMult",
435  TH2::Class()));
436  fCentNPart = static_cast<TH2*>(GetOutputObject(output, "centNPart",
437  TH2::Class()));
438  fCentNBin = static_cast<TH2*>(GetOutputObject(output, "centNBin",
439  TH2::Class()));
440  fCentList = static_cast<TList*>(GetOutputObject(output, "byCent",
441  TList::Class()));
442  if (!fCentList || !fCentAll || !fCentAcc) {
443  Warning("Finalize", "Missing stack and histograms");
444  return false;
445  }
446  Info("Terminate", "Accepted %d/%d=%6.2f%% events",
447  int(fCentAcc->GetEntries()), int(fCentAll->GetEntries()),
448  100*fCentAcc->GetEntries()/fCentAll->GetEntries());
449 
450  THStack* stack = new THStack("all", "All");
451  TList* hists = fCentList; // ->GetHists();
452  TObjLink* link = hists->FirstLink();
453  Int_t bin = 1;
454  Long64_t sum = 0;
455  Long64_t total = 0;
456  Long64_t all = fCentAll->GetEntries();
457  while (link) {
458  TObject* o = link->GetObject();
459  if (!o) {
460  link = link->Next();
461  bin++;
462  continue;
463  }
464  TH1* h = static_cast<TH1*>(o);
465  Int_t n = h->GetBinContent(0);
466  Int_t m = fCentAcc->GetBinContent(bin);
467  total += m;
468  h->SetBinContent(0,0);
469  printf("%9d (%9d) events in bin %s ...", n, m, o->GetTitle());
470  if (n < minEvents) {
471  // Too few event, remove this
472  TObjLink* tmp = link->Next();
473  hists->Remove(link);
474  link = tmp;
475  delete o;
476  Printf(" removed");
477  bin++;
478  continue;
479  }
480  sum += m;
481  // Scale
482  h->Scale(1. / n, "width");
483  stack->Add(h);
484  Printf(" scaled");
485  link = link->Next();
486  bin++;
487  }
488  Printf("ana/acc/all: %9lld/%9lld/%9lld [%6.2f%%/%6.2f%%]",
489  sum, total, all, float(100*sum)/total, float(100*total)/all);
490  output->Add(stack);
491  }
493 };
494 #endif
495 
496 //
497 // EOF
498 //
void SetCentAxis(const char *spec)
double Double_t
Definition: External.C:58
void ModHist(TH1 *h, Double_t low, Double_t high)
TObject * GetOutputObject(TCollection *output, const char *name, TClass *cls)
long long Long64_t
Definition: External.C:43
void SetCentAxis(Int_t n, Double_t *edges)
virtual const char * HistName(Double_t low, Double_t high) const
virtual const char * HistTitle(Double_t low, Double_t high) const
void CreateHistos(TCollection *output, TH1D *(*callback)())
void SetCentAxis(Int_t n, Double_t low, Double_t high)
int Int_t
Definition: External.C:63
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
ClassDef(FastCentHelper, 1)
void CreateDiagnostics(TCollection *output, TH1 *centHist)
Bool_t Finalize(TCollection *output, Long_t minEvents)
Definition: External.C:220
Int_t CheckCentrality(Double_t cent, Double_t mult, Double_t b, Double_t nPart, Double_t nBin)
bool Bool_t
Definition: External.C:53
Definition: External.C:196