AliPhysics  cc1c0ba (cc1c0ba)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CentSysErr.C
Go to the documentation of this file.
1 #include <TFile.h>
2 #include <TCollection.h>
3 #include <TList.h>
4 #include <TClass.h>
5 #include <THStack.h>
6 #include <TAxis.h>
7 #include <TH1.h>
8 #include <TH2.h>
9 #include <TString.h>
10 #include <TError.h>
11 #include <TCanvas.h>
12 #include <TLegend.h>
13 #include <TLegendEntry.h>
14 #include <TStyle.h>
15 #include <TParameter.h>
16 #include <TLatex.h>
17 #include <TGraphAsymmErrors.h>
18 #include <TF1.h>
19 #include <TMath.h>
20 
21 //____________________________________________________________________
22 TObject*
23 GetObject(TDirectory* d, const char* name, TClass* cls=0, Bool_t verb=true)
24 {
25  if (!d) {
26  Error("GetOject", "No directory");
27  return 0;
28  }
29 
30  TObject* o = d->Get(name);
31  if (!o) {
32  if (verb)
33  Error("GetObject", "No object %s in directory %s", name, d->GetName());
34  return 0;
35  }
36 
37  if (!cls) return o;
38 
39  if (!o->IsA()->InheritsFrom(cls)) {
40  if (verb)
41  Error("GetObject", "%s from %s is not a %s, but a %s",
42  o->GetName(), d->GetName(), cls->GetName(), o->ClassName());
43  return 0;
44  }
45 
46  return o;
47 }
48 
49 //____________________________________________________________________
50 TObject*
51 GetObject(TCollection* d, const char* name, TClass* cls=0, Bool_t verb=true)
52 {
53  if (!d) {
54  Error("GetOject", "No collection");
55  return 0;
56  }
57 
58  TObject* o = d->FindObject(name);
59  if (!o) {
60  if (verb)
61  Error("GetObject", "No object %s in collection %s", name, d->GetName());
62  return 0;
63  }
64 
65  if (!cls) return o;
66 
67  if (!o->IsA()->InheritsFrom(cls)) {
68  if (verb)
69  Error("GetObject", "%s from %s is not a %s, but a %s",
70  o->GetName(), d->GetName(), cls->GetName(), o->ClassName());
71  return 0;
72  }
73 
74  return o;
75 }
76 //____________________________________________________________________
77 TCollection*
78 GetCollection(TDirectory* d, const char* name)
79 {
80  return static_cast<TCollection*>(GetObject(d, name, TCollection::Class()));
81 }
82 //____________________________________________________________________
83 TCollection*
84 GetCollection(TCollection* d, const char* name)
85 {
86  return static_cast<TCollection*>(GetObject(d, name, TCollection::Class()));
87 }
88 //____________________________________________________________________
89 TH2*
90 GetH2(TCollection* c, const char* name)
91 {
92  return static_cast<TH2*>(GetObject(c, name, TH2::Class()));
93 }
94 //____________________________________________________________________
95 TH1*
96 GetH1(TCollection* c, const char* name)
97 {
98  return static_cast<TH1*>(GetObject(c, name, TH1::Class()));
99 }
100 //____________________________________________________________________
101 THStack*
102 GetStack(TCollection* c, const char* name, Bool_t verb=true)
103 {
104  return static_cast<THStack*>(GetObject(c, name, THStack::Class(),verb));
105 }
106 //____________________________________________________________________
107 TParameter<double>*
108 GetParam(TCollection* c, const char* name, Bool_t verb=false)
109 {
110  return static_cast<TParameter<double>*>(GetObject(c, name,
111  TParameter<double>::Class(),
112  verb));
113 }
114 
115 
116 //____________________________________________________________________
117 void
118 ProcessOne(const char* meth,
119  Int_t style,
120  TList* stacks,
121  TList* ratios,
122  TList* mins,
123  TList* maxs,
124  TList* avgs,
125  TAxis*& centAxis)
126 {
127  TFile* file = TFile::Open(Form("PbPb_2760_dndeta_%s_nosec_20140512_1628/"
128  "forward_dndeta.root", meth), "READ");
129  TCollection* res = GetCollection(file, "ForwarddNdetaResults");
130  Bool_t first = false;
131  if (!centAxis) {
132  centAxis = static_cast<TAxis*>(GetObject(res, "centAxis", TAxis::Class()));
133  Info("ProcessOne", "Got centrality axis %p", centAxis);
134  first = true;
135  }
136 
137  Info("ProcessOne", "Doing stuff for %s (%sfirst)",
138  meth, (first ? "" : "not "));
139  for (Int_t i = 1; i <= centAxis->GetNbins(); i++) {
140  Double_t c1 = centAxis->GetBinLowEdge(i);
141  Double_t c2 = centAxis->GetBinUpEdge(i);
142  if (meth[0] == 'Z' && c1 > 30) break;
143 
144  TString binName(Form("cent%03dd%02d_%03dd%02d",
145  Int_t(c1), Int_t(c1 *100) % 100,
146  Int_t(c2), Int_t(c2 *100) % 100));
147  TCollection* bin = GetCollection(res, binName);
148  if (!bin) continue;
149 
150  TH1* here = GetH1(bin, "dndetaForward");
151  if (!here) continue;
152 
153  THStack* stack = GetStack(stacks, binName, false);
154  THStack* ratio = GetStack(ratios, binName, false);
155  TParameter<double>* pmin = GetParam(mins, binName, false);
156  TParameter<double>* pmax = GetParam(maxs, binName, false);
157  TParameter<double>* pavg = GetParam(avgs, binName, false);
158  if (!stack) {
159  Info("ProcessOne", "Creating stacks %s", binName.Data());
160  stack = new THStack(binName, Form("%3d-%3d%%", Int_t(c1), Int_t(c2)));
161  stack->SetUniqueID(here->GetMarkerColor());
162  stacks->AddAt(stack, i-1);
163  ratio = static_cast<THStack*>(stack->Clone());
164  ratios->AddAt(ratio, i-1);
165  pmin = new TParameter<double>(binName, +1000.);
166  pmax = new TParameter<double>(binName, -1000.);
167  pavg = new TParameter<double>(binName, -1000.);
168  mins->AddAt(pmin, i-1);
169  maxs->AddAt(pmax, i-1);
170  avgs->AddAt(pavg, i-1);
171  }
172  here->SetTitle(meth);
173  here->SetMarkerStyle(style);
174  stack->Add(here);
175 
176 
177  if (first) continue;
178 
179  TH1* denom = static_cast<TH1*>(stack->GetHists()->At(0));
180  TH1* num = static_cast<TH1*>(here->Clone(Form("%s/%s",
181  meth, denom->GetName())));
182  num->SetDirectory(0);
183  num->Divide(denom);
184 
185  Double_t min = pmin->GetVal();
186  Double_t max = pmax->GetVal();
187  Double_t sum = 0;
188  Int_t cnt = 0;
189  for (Int_t j = 1; j <= num->GetNbinsX(); j++) {
190  Double_t c = num->GetBinContent(j);
191  if (c < 1e-6) continue;
192  min = TMath::Min(c, min);
193  max = TMath::Max(c, max);
194  sum += TMath::Abs(1-c);
195  cnt++;
196  num->SetBinContent(j, c+(centAxis->GetNbins()-i));
197  }
198  ratio->Add(num);
199  pmin->SetVal(min);
200  pmax->SetVal(max);
201  sum /= cnt;
202  pavg->SetVal(TMath::Max(pavg->GetVal(), sum));
203 
204  Printf("Method %10s %16s: min=%7.5f max=%7.5f avg=%7.5f",
205  meth,stack->GetTitle(), min, max, sum);
206 
207  }
208 }
209 
210 
211 void
212 CentSysErr(Bool_t div=false, Bool_t zem=false)
213 {
214  const char* dndeta = "1/#it{N} d#it{N}_{ch}/d#it{#eta}";
215  TList* stacks = new TList;
216  TList* ratios = new TList;
217  TList* mins = new TList;
218  TList* maxs = new TList;
219  TList* avgs = new TList;
220  TAxis* centAxis = 0;
221  TString all = "";
222 
223  const char* meths[] = { "V0M", "CL1", "TRK", (zem ? "ZEMVSZDC" : 0), 0 };
224  Int_t styles[] = { 20, 21, 22, 23, 0 };
225 
226  const char** pmeth = meths;
227  Int_t* pstyle = styles;
228 
229  TCanvas* c = new TCanvas("c", "C", 1200, 800);
230  c->SetTopMargin(0.01);
231  c->SetRightMargin(0.01);
232  c->Divide(1,2,0,0);
233  c->GetPad(1)->SetRightMargin(0.01);
234  c->GetPad(2)->SetRightMargin(0.01);
235  gStyle->SetOptTitle(0);
236 
237  TLegend* mleg = new TLegend(0.35, (meths[3] ? 0.05 : 0.3), 0.65, 0.95);
238  mleg->SetFillStyle(0);
239  mleg->SetBorderSize(0);
240  mleg->SetTextFont(42);
241  TLegend* leg = new TLegend(0.35, 0.15, 0.55, 0.95);
242  leg->SetFillStyle(0);
243  leg->SetBorderSize(0);
244  leg->SetTextFont(42);
245  leg->SetNColumns(2);
246 
247 
248  while (*pmeth) {
249  ProcessOne(*pmeth, *pstyle, stacks, ratios, mins, maxs, avgs, centAxis);
250  TLegendEntry* e = mleg->AddEntry("", *pmeth, "p");
251  e->SetMarkerStyle(*pstyle);
252  all.Append(Form("_%s", *pmeth));
253  pmeth++;
254  pstyle++;
255 
256  }
257 
258  TGraphAsymmErrors* trend = new TGraphAsymmErrors(centAxis->GetNbins());
259  trend->SetMarkerStyle(20);
260  trend->SetFillColor(kBlue-9);
261  trend->SetFillStyle(3001);
262  trend->SetMarkerColor(kBlue-7);
263  trend->SetLineColor(kBlue-9);
264  TIter nextS(stacks);
265  TIter nextR(ratios);
266  TIter nextN(mins);
267  TIter nextM(maxs);
268  TIter nextA(avgs);
269  Bool_t first = true;
270  THStack* stack = 0;
271  THStack* ratio = 0;
272  TParameter<double>* pmin = 0;
273  TParameter<double>* pmax = 0;
274  TParameter<double>* pavg = 0;
275  Int_t cnt = 0;
276  Double_t least = +1000;
277  Double_t largest = -1000;
278  while ((stack = static_cast<THStack*>(nextS())) &&
279  (ratio = static_cast<THStack*>(nextR())) &&
280  (pmin = static_cast<TParameter<double>*>(nextN())) &&
281  (pmax = static_cast<TParameter<double>*>(nextM())) &&
282  (pavg = static_cast<TParameter<double>*>(nextA()))
283  ) {
284  c->cd(1);
285  stack->Draw(first ? "nostack" : "nostack same");
286 
287  Int_t y = centAxis->GetNbins()-cnt++;
288  UInt_t col = stack->GetUniqueID();
289  TLegendEntry* e = leg->AddEntry("",
290  Form("%10s (+%d)",
291  stack->GetTitle(), y-1),
292  "f");
293  e->SetFillColor(col);
294  e->SetFillStyle(1001);
295  e->SetLineColor(col);
296  c->cd(2);
297  ratio->Draw(first ? "nostack" : "nostack same");
298 
299  Double_t min = 100*(1 - pmin->GetVal());
300  Double_t max = 100*(pmax->GetVal()-1);
301  Double_t avg = 100*pavg->GetVal();
302  Double_t c1 = centAxis->GetBinLowEdge(cnt);
303  Double_t c2 = centAxis->GetBinUpEdge(cnt);
304  Double_t cc = (c2+c1)/2;
305  trend->SetPoint(cnt-1, cc, avg);
306  trend->SetPointError(cnt-1, cc-c1, c2-cc,
307  (min > 0 ? TMath::Abs(avg-min) : 0), max-avg);
308  least = TMath::Min(avg, least);
309  largest = TMath::Max(avg, largest);
310 
311  TString errs = Form("#pm%4.1f%%", avg);
312  if (div) {
313  if (min > 0) errs = Form("{}^{+%4.1f}_{-%4.1f}%%", max, min);
314  else errs = Form("+%4.1f%%", TMath::Max(min,max));
315  }
316  if (leg->GetNColumns() > 1) {
317  e = leg->AddEntry("", errs, "");
318  e->SetFillColor(kWhite);
319  e->SetFillStyle(0);
320  e->SetLineColor(kWhite);
321  e->SetMarkerColor(kWhite);
322  }
323  else {
324  TLatex* ltx = new TLatex(5.9, y, errs);
325  ltx->SetTextAlign(32);
326  ltx->SetTextFont(42);
327  ltx->SetTextSize(0.04);
328  ltx->SetTextColor(col+1);
329  ltx->Draw();
330  }
331 
332  if (!first) continue;
333 
334  ratio->SetMaximum(centAxis->GetNbins()+1.3);
335  ratio->SetMinimum(.1);
336 
337  Double_t ts = 0.06;
338  TH1* hr = ratio->GetHistogram();
339  hr->SetXTitle("#it{#eta}");
340  hr->SetYTitle("Ratio to V0M");
341  hr->GetXaxis()->SetTitleSize(ts);
342  hr->GetYaxis()->SetTitleSize(ts);
343  hr->GetXaxis()->SetLabelSize(ts);
344  hr->GetYaxis()->SetLabelSize(ts);
345  hr->GetXaxis()->SetTitleOffset(1-4*ts);
346  hr->GetYaxis()->SetTitleOffset(1-5*ts);
347  hr->GetXaxis()->SetNdivisions(10);
348  hr->GetYaxis()->SetNdivisions(10);
349 
350  TH1* hs = stack->GetHistogram();
351  hs->SetXTitle("#it{#eta}");
352  hs->SetYTitle(dndeta);
353  hs->GetXaxis()->SetTitleSize(ts);
354  hs->GetYaxis()->SetTitleSize(ts);
355  hs->GetXaxis()->SetLabelSize(ts);
356  hs->GetYaxis()->SetLabelSize(ts);
357  hs->GetXaxis()->SetTitleOffset(1-4*ts);
358  hs->GetYaxis()->SetTitleOffset(1-5*ts);
359  hs->GetXaxis()->SetNdivisions(10);
360  hs->GetYaxis()->SetNdivisions(10);
361 
362 
363  first = false;
364  }
365 
366  c->cd(1);
367  mleg->Draw();
368  c->cd(2);
369  leg->Draw();
370 
371  c->Print(Form("cent_syserr_%s%s.pdf", (div ? "div" : "avg"), all.Data()));
372  c->SaveAs(Form("cent_syserr_%s%s.root", (div ? "div" : "avg"), all.Data()));
373 
374  TCanvas* aux = new TCanvas("aux", "aix");
375  aux->SetTopMargin(0.01);
376  aux->SetRightMargin(0.01);
377 
378  trend->SetTitle();
379  trend->Draw("a3p");
380  trend->GetHistogram()->SetXTitle("Centrality [%]");
381  trend->GetHistogram()->SetYTitle(Form("#delta#left[%s#right] [%%]",dndeta));
382 
383  TF1* f = new TF1("f", "[0]*x*x+[1]", 0, 100);
384  f->SetLineColor(kMagenta+2);
385  Double_t cl = centAxis->GetBinCenter(1);
386  Double_t ch = centAxis->GetBinCenter(centAxis->GetNbins());
387  f->SetParameters(largest/ch/ch, least);
388  f->Draw("same");
389 
390  TLatex* ll = new TLatex(0.98, 0.98,
391  Form("Sys. Uncertainty on %s from Centrality",
392  dndeta));
393  ll->SetTextAlign(33);
394  ll->SetTextSize(0.05);
395  ll->SetTextFont(42);
396  ll->SetNDC();
397  ll->Draw();
398 
399  TLatex* fl = new TLatex(0.2, 0.8,
400  Form("#delta#approx%5.1fc^{2}+%5.1f",
401  largest/ch/ch*100*100,
402  least));
403  fl->SetTextAlign(13);
404  fl->SetTextSize(0.04);
405  fl->SetTextFont(42);
406  fl->SetNDC();
407  fl->Draw();
408 
409  aux->Print(Form("cent_syserr_%s%s_trend.pdf",
410  (div ? "div" : "avg"), all.Data()));
411 }
412 
413 
414 
415 
cc
const Color_t cc[]
Definition: DrawKs.C:1
Double_t
double Double_t
Definition: External.C:58
GetObject
TObject * GetObject(TDirectory *d, const char *name, TClass *cls=0, Bool_t verb=true)
Definition: CentSysErr.C:23
CentSysErr
void CentSysErr(Bool_t div=false, Bool_t zem=false)
Definition: CentSysErr.C:212
GetParam
TParameter< double > * GetParam(TCollection *c, const char *name, Bool_t verb=false)
Definition: CentSysErr.C:108
c
TCanvas * c
Definition: TestFitELoss.C:172
stack
AliStack * stack
Definition: AliMCTreeTools.cxx:72
Int_t
int Int_t
Definition: External.C:63
UInt_t
unsigned int UInt_t
Definition: External.C:33
GetH1
TH1 * GetH1(TCollection *c, const char *name)
Definition: CentSysErr.C:96
ProcessOne
void ProcessOne(const char *meth, Int_t style, TList *stacks, TList *ratios, TList *mins, TList *maxs, TList *avgs, TAxis *&centAxis)
Definition: CentSysErr.C:118
GetH2
TH2 * GetH2(TCollection *c, const char *name)
Definition: CentSysErr.C:90
TH2
Definition: External.C:220
file
TFile * file
TList with histograms for a given trigger.
Definition: DrawAnaCaloTrackQA.C:89
GetStack
THStack * GetStack(TCollection *c, const char *name, Bool_t verb=true)
Definition: CentSysErr.C:102
Bool_t
bool Bool_t
Definition: External.C:53
GetCollection
TCollection * GetCollection(TDirectory *d, const char *name)
Definition: CentSysErr.C:78
TH1
Definition: External.C:196