AliPhysics  8d00e07 (8d00e07)
ExtractGSEs.C
Go to the documentation of this file.
1 static Int_t PbPbBin(Double_t c1, Double_t c2)
2 {
3  Double_t c = (c1+c2) / 2;
4  if (c < 5) return 0;
5  else if (c < 10) return 1;
6  else if (c < 20) return 2;
7  else if (c < 30) return 3;
8  else if (c < 40) return 4;
9  else if (c < 50) return 5;
10  else if (c < 60) return 6;
11  else if (c < 70) return 7;
12  else if (c < 80) return 8;
13  else if (c < 90) return 9;
14  return 10;
15 }
24 static Color_t PbPbColor(Double_t c1, Double_t c2)
25 {
26  const Color_t cc[] = { kMagenta+2,
27  kBlue+2,
28  kAzure-1, // 10,
29  kCyan+2,
30  kGreen+1,
31  kSpring+5,//+10,
32  kYellow+1,
33  kOrange+5,//+10,
34  kRed+1,
35  kPink+5,//+10,
36  kBlack };
37  Int_t bin = PbPbBin(c1,c2);
38  return cc[bin];
39 }
40 //____________________________________________________________________
41 TObject*
42 GetObject(TDirectory* d, const char* name, TClass* cls=0, Bool_t verb=true)
43 {
44  if (!d) {
45  Error("GetOject", "No directory");
46  return 0;
47  }
48 
49  TObject* o = d->Get(name);
50  if (!o) {
51  if (verb)
52  Error("GetObject", "No object %s in directory %s", name, d->GetName());
53  return 0;
54  }
55 
56  if (!cls) return o;
57 
58  if (!o->IsA()->InheritsFrom(cls)) {
59  if (verb)
60  Error("GetObject", "%s from %s is not a %s, but a %s",
61  o->GetName(), d->GetName(), cls->GetName(), o->ClassName());
62  return 0;
63  }
64 
65  return o;
66 }
67 
68 //____________________________________________________________________
69 TObject*
70 GetObject(TCollection* d, const char* name, TClass* cls=0, Bool_t verb=true)
71 {
72  if (!d) {
73  Error("GetOject", "No collection");
74  return 0;
75  }
76 
77  TObject* o = d->FindObject(name);
78  if (!o) {
79  if (verb)
80  Error("GetObject", "No object %s in collection %s", name, d->GetName());
81  return 0;
82  }
83 
84  if (!cls) return o;
85 
86  if (!o->IsA()->InheritsFrom(cls)) {
87  if (verb)
88  Error("GetObject", "%s from %s is not a %s, but a %s",
89  o->GetName(), d->GetName(), cls->GetName(), o->ClassName());
90  return 0;
91  }
92 
93  return o;
94 }
95 //____________________________________________________________________
97 GetCollection(TDirectory* d, const char* name, Bool_t verb=true)
98 {
99  return static_cast<TCollection*>(GetObject(d, name, TCollection::Class(), verb));
100 }
101 //____________________________________________________________________
103 GetCollection(TCollection* d, const char* name, Bool_t verb=true)
104 {
105  return static_cast<TCollection*>(GetObject(d, name, TCollection::Class(), verb));
106 }
107 //____________________________________________________________________
108 TH2*
109 GetH2(TCollection* c, const char* name, Bool_t verb=true)
110 {
111  return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verb));
112 }
113 //____________________________________________________________________
114 TH1*
115 GetH1(TCollection* c, const char* name, Bool_t verb=true)
116 {
117  return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verb));
118 }
119 //____________________________________________________________________
120 THStack*
121 GetStack(TCollection* c, const char* name, Bool_t verb=true)
122 {
123  return static_cast<THStack*>(GetObject(c, name, THStack::Class(),verb));
124 }
125 //____________________________________________________________________
126 TAxis*
127 GetAxis(TCollection* c, const char* name, Bool_t verb=true)
128 {
129  return static_cast<TAxis*>(GetObject(c, name, TAxis::Class(),verb));
130 }
131 //____________________________________________________________________
133 GetParam(TCollection* c, const char* name, Bool_t verb=false)
134 {
135  return static_cast<TParameter<double>*>(GetObject(c, name,
137  verb));
138 }
139 
140 //====================================================================
141 void
142 AddPath(const TString& dir, Bool_t before=true)
143 {
144  const char* now = gROOT->GetMacroPath();
145  const char* fst = (before ? dir.Data() : now);
146  const char* lst = (before ? now : dir.Data());
147  gROOT->SetMacroPath(Form("%s:%s",fst,lst));
148 
149  gSystem->AddIncludePath(Form("-I%s",dir.Data()));
150 }
151 
152 //____________________________________________________________________
153 void
154 AddRelPath(const char* rel, Bool_t first=true)
155 {
156  TString alp(Form("$ALICE_PHYSICS/PWGLF/FORWARD/analysis2/%s",rel));
157  AddPath(alp, first);
158  if (gSystem->Getenv("ANA_SRC")) {
159  TString mine(Form("$ANA_SRC/%s", rel));
160  AddPath(mine, first);
161  }
162 }
163 
164 
165 
166 //====================================================================
167 TH1* Rebin(TH1* h, Int_t rebin, Bool_t cutEdges) const
168 {
169  if (rebin <= 1) return h;
170 
171  Int_t nBins = h->GetNbinsX();
172  if(nBins % rebin != 0) {
173  Warning("Rebin", "Rebin factor %d is not a devisor of current number "
174  "of bins %d in the histogram %s", rebin, nBins, h->GetName());
175  return;
176  }
177 
178  // Make a copy
179  TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
180  tmp->Rebin(rebin);
181  tmp->SetDirectory(0);
182  tmp->Reset();
183  // The new number of bins
184  Int_t nBinsNew = nBins / rebin;
185  for(Int_t i = 1;i<= nBinsNew; i++) {
186  Double_t content = 0;
187  Double_t sumw = 0;
188  Double_t wsum = 0;
189  Int_t nbins = 0;
190  for(Int_t j = 1; j<=rebin;j++) {
191  Int_t bin = (i-1)*rebin + j;
192  Double_t c = h->GetBinContent(bin);
193 
194  if (c <= 0) continue;
195 
196  if ((cutEdges)) {
197  if (h->GetBinContent(bin+1)<=0 ||
198  h->GetBinContent(bin-1)<=0) {
199  Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
200  bin, c, h->GetName(),
201  bin+1, h->GetBinContent(bin+1),
202  bin-1, h->GetBinContent(bin-1));
203  continue;
204  }
205  }
206  Double_t e = h->GetBinError(bin);
207  Double_t w = 1 / (e*e); // 1/c/c
208  content += c;
209  sumw += w;
210  wsum += w * c;
211  nbins++;
212  }
213 
214  if(content > 0 && nbins > 1 ) {
215  tmp->SetBinContent(i, wsum / sumw);
216  tmp->SetBinError(i,1./TMath::Sqrt(sumw));
217  }
218  }
219 
220  // Finally, rebin the histogram, and set new content
221  return tmp;
222 }
223 
224 
225 
226 //____________________________________________________________________
227 TObject*
229  const TString& dir,
230  const TString& name,
231  Int_t rebin,
232  Double_t eff,
233  void* oa,
234  Bool_t cutEdges=false,
235  TLegend* l=0)
236 {
237  SysErrorAdder* adder = reinterpret_cast<SysErrorAdder*>(oa);
238 
239  TCollection* sub = GetCollection(c, dir);
240  if (!sub) return 0;
241 
242  TH1* hist = GetH1(sub, name);
243  if (!hist) return 0;
244 
245  hist = Rebin(hist, rebin, cutEdges);
246  hist->SetMarkerColor(kBlack);
247  hist->SetFillColor(kBlack);
248  hist->SetLineColor(kBlack);
249 
250  if (dir.BeginsWith("cent"))
251  hist->SetName(Form("%s_%s",hist->GetName(),dir.Data()));
252 
253  Info("DoOne","Make graph from %s with eff=%f", hist->GetName(), eff);
254  GraphSysErr* gse = adder->Make(hist, l, eff, true);
255  if (!gse) return 0;
256 
258  gse->SetSumFillColor(gse->GetMarkerColor());
259  gse->SetSumFillStyle(3002);
260  gse->SetCommonSumOption(GraphSysErr::kBox);
261  gse->SetCommonSumFillColor(gse->GetMarkerColor());
262  gse->SetCommonSumFillStyle(3001);
263 
264  return gse;
265 }
266 
267 //____________________________________________________________________
268 void ModOne(TObject* o, TString& sys, UInt_t sNN, TString trg, Color_t col)
269 {
270  GraphSysErr* g = static_cast<GraphSysErr*>(o);
271  if (!sys.EqualTo("pp"))
272  g->AddQualifier("SQRT(S)/NUCLEON IN GEV", Form("%d.0", sNN), true);
273  else {
274  g->AddQualifier("SQRT(S) IN GEV", Form("%d.0", sNN), true);
275  g->AddQualifier("TRIGGER", trg);
276  }
277 
278  g->SetMarkerColor(col);
279  g->SetSumFillColor(kWhite);
280  g->SetSumFillStyle(0);
281  g->SetSumLineColor(col);
282  g->SetCommonSumFillColor(col);
283  g->SetCommonSumFillStyle(1001);
284  g->SetCommonSumLineColor(col);
285  g->Print("sys qual key");
286 }
287 
288 //____________________________________________________________________
289 TH1*
290 DrawOne(TObject* o, Option_t* opt, Double_t& min, Double_t& max)
291 {
292  GraphSysErr* g = static_cast<GraphSysErr*>(o);
293  g->Draw(opt);
294  Double_t mn, mx;
295  g->GetMinMax(opt, mn, mx);
296  min = TMath::Min(min, mn);
297  max = TMath::Max(max, mx);
298  TMultiGraph* mg = g->GetMulti();
299  // Info("DrawOne", "Got Multigraph %p", mg);
300  TH1* hist = (mg ? mg->GetHistogram() : 0);
301  // Info("DrawOne", "Got frame histogram %p", hist);
302  return hist;
303 }
304 
305 
306 //____________________________________________________________________
307 TList*
308 ExtractGSEs(const char* filename="forward_dndeta.root",
309  Int_t rebin=5,
310  Double_t eff=1,
311  Bool_t raw=false,
312  Bool_t cutEdges=false,
313  const char* name="Forward")
314 {
315  if (!gROOT->GetClass("GraphSysErr")) {
316  AddPath("$HOME/GraphSysErr");
317  AddRelPath("gse",false);
318 
319  gROOT->LoadMacro("GraphSysErr.C+g");
320  }
321  if (!gROOT->GetClass("SysErrorAdder")) {
322  AddRelPath("dndeta", true);
323  gROOT->LoadMacro("SysErrorAdder.C+g");
324  }
325 
326  TFile* file = TFile::Open(filename,"READ");
327  if (!file) {
328  Warning("ExtractGSEs", "Failed to open %s", filename);
329  return 0;
330  }
331 
332  TCollection* results = GetCollection(file, Form("%sdNdetaResults",name));
333  if (!results) return 0;
334 
335 
336  // TH1* empH = GetH1 (results, "empirical");
337  TObject* empO = GetObject(results, "empirical");
338  Bool_t emp = (empO != 0 && !raw);
339  TObject* osNN = GetObject(results, "sNN");
340  TObject* osys = GetObject(results, "sys");
341  TObject* otrg = GetObject(results, "trigger");
342  TObject* ocentM = GetObject(results, "centEstimator");
343  TAxis* centA = GetAxis (results, "centAxis");
344  TString sys = osys->GetTitle();
345  UShort_t sNN = osNN->GetUniqueID();
346  TString mth = (ocentM && centA ? ocentM->GetTitle() : "");
347  if (mth.EqualTo("none",TString::kIgnoreCase) ||
348  (centA && centA->GetNbins() < 1)) mth = "";
349  TString trg = (!mth.IsNull() && centA ? "CENT" : otrg->GetTitle());
350  TString hname = Form("dndeta%s%s",name,emp?"Emp":"");
351 
352  Printf("Read settings\n"
353  " Empirical: %s (%p)\n"
354  " System: %s (%p)\n"
355  " Energy: %d (%s - %p)\n"
356  " Trigger: %s (%s - %p)\n"
357  " Efficiency: %6.4f\n"
358  " Centrality: %s (%p)\n"
359  " Axis: %p",
360  (emp ? "yes" : "no"), empO,
361  sys.Data(), osys,
362  sNN, (osNN ? osNN->GetTitle() : ""), osNN,
363  trg.Data(), eff, otrg->GetTitle(), otrg, mth.Data(),
364  ocentM, centA);
365  if (centA) Printf(" %d bins between %f and %f",
366  centA->GetNbins(), centA->GetXmin(), centA->GetXmax());
367 
368  // Possibly change trigger
369  if (trg.EqualTo("MBOR") && eff > 0 && eff < 1) trg="INEL";
370  if (trg.EqualTo("V0AND") && eff > 0 && eff < 1) trg="NSD";
371  if (trg.EqualTo("VISX") && eff > 0 && eff < 1) trg="NSD";
372  TString strg(trg);
373  if (strg.Contains("INEL>0")) trg = "INELGt0";
374 
375 
376  SysErrorAdder* adder = SysErrorAdder::Create(trg,sys,sNN,mth);
377  if (!adder) {
378  Warning("ExtractGSEs", "Failed to make adder");
379  return 0;
380  }
381 
382  TCanvas* c = new TCanvas("c", "C");
383  c->SetTopMargin(0.01);
384  c->SetRightMargin(0.20);
385  TLegend* l = new TLegend(0.8,0.1,0.99,0.99,
386  Form("%s @ %s (%s)",
387  sys.Data(), osNN->GetTitle(), trg.Data()));
388  l->SetFillColor(0);
389  l->SetFillStyle(0);
390  l->SetBorderSize(0);
391 
392 
393  TList* ret = new TList;
394  Bool_t first = true;
395  TH1* frame = 0;
396  Double_t min = 100000, max = -1000000;
397  if (!centA || centA->GetNbins() < 1 || mth.IsNull()) {
398  Info("ExtractGSEs", "Doing pp-like extraction"
399  " all bin, %s, rebin=%d, eff=%f, %p, %s",
400  hname.Data(), rebin, eff, adder,
401  (cutEdges ? "cut edges" : "edges stay"));
402  TObject* all = DoOne(results,"all",hname,rebin,eff,adder,cutEdges,l);
403  if (!all) {
404  Warning("ExtractGSEs", "Nothing returned from DoOne(\"all\"...)");
405  return 0;
406  }
407  GraphSysErr* gg = static_cast<GraphSysErr*>(all);
408  ModOne(gg, sys, sNN, trg, gg->GetMarkerColor());
409  gg->SetName(strg);
410  ret->Add(all);
411  frame = DrawOne(all, "SUM QUAD AXIS", min, max);
412  }
413  else {
414  for (Int_t i = 1; i <= centA->GetNbins(); i++) {
415  Double_t low = centA->GetBinLowEdge(i);
416  Double_t high = centA->GetBinUpEdge(i);
417  TString dir;
418  dir.Form("cent%03dd%02d_%03dd%02d",
419  Int_t(low), Int_t(low *100)%100,
420  Int_t(high), Int_t(high*100)%100);
421  TObject* g = DoOne(results, dir, hname, rebin, eff, adder, cutEdges,
422  (first ? l : 0));
423  if (!g) continue;
424  ret->Add(g);
425 
426  GraphSysErr* gse = static_cast<GraphSysErr*>(g);
427  Color_t col = PbPbColor(low, high);
428  ModOne(gse, sys, sNN, trg, col);
429  if (first) frame = DrawOne(g, "SUM QUAD AXIS", min, max);
430  else DrawOne(g, "SUM QUAD", min, max);
431  TLegendEntry* e = l->AddEntry("", Form("%3d-%3d%%",
432  int(low), int(high)),"F");
433  e->SetFillColor(col);
434  e->SetFillStyle(1001);
435  e->SetLineColor(col);
436  e->SetMarkerColor(col);
437  first = false;
438  }
439  }
440  if (!frame) {
441  Warning("ExtractGSEs", "No frame given");
442  }
443  else {
444  frame->SetMinimum(0.9*min);
445  frame->SetMaximum(1.1*max);
446  }
447  l->Draw();
448 
449  TString outName;
450  if (mth.EqualTo("V0M", TString::kIgnoreCase) &&
451  sys.EqualTo("PbPb",TString::kIgnoreCase))
452  mth = "";
453  if (raw) mth = "RAW";
454  outName.Form("%s_%05d_%s%s.root",sys.Data(),sNN,strg.Data(),mth.Data());
455  TFile* out = TFile::Open(outName,"RECREATE");
456  ret->Write("container",TObject::kSingleKey);
457  out->Write();
458 
459  Info("ExtractGSEs", "Written to %s", outName.Data());
460  return ret;
461 }
462 
TH1 * DrawOne(TObject *o, Option_t *opt, Double_t &min, Double_t &max)
Definition: ExtractGSEs.C:290
const char * filename
Definition: TestFCM.C:1
const Color_t cc[]
Definition: DrawKs.C:1
static Int_t PbPbBin(Double_t c1, Double_t c2)
Definition: ExtractGSEs.C:1
static SysErrorAdder * Create(const TString &t, const TString &s, UShort_t e, const TString &c)
double Double_t
Definition: External.C:58
void Draw(Option_t *option="")
Definition: GraphSysErr.C:791
virtual GraphSysErr * Make(TH1 *h, TLegend *l, Double_t eff=1, Bool_t verb=false)
TH1 * Rebin(TH1 *h, Int_t rebin, Bool_t cutEdges) const
Definition: ExtractGSEs.C:167
TSystem * gSystem
TObject * DoOne(TCollection *c, const TString &dir, const TString &name, Int_t rebin, Double_t eff, void *oa, Bool_t cutEdges=false, TLegend *l=0)
Definition: ExtractGSEs.C:228
void SetCommonSumLineColor(Color_t color)
Definition: GraphSysErr.C:4503
TCollection * GetCollection(TDirectory *d, const char *name, Bool_t verb=true)
Definition: ExtractGSEs.C:97
TH2 * GetH2(TCollection *c, const char *name, Bool_t verb=true)
Definition: ExtractGSEs.C:109
TCanvas * c
Definition: TestFitELoss.C:172
void AddQualifier(const TString &key, const TString &value, Bool_t replace=false)
Definition: GraphSysErr.C:4759
THStack * GetStack(TCollection *c, const char *name, Bool_t verb=true)
Definition: ExtractGSEs.C:121
void AddRelPath(const char *rel, Bool_t first=true)
Definition: ExtractGSEs.C:154
TH1 * GetH1(TCollection *c, const char *name, Bool_t verb=true)
Definition: ExtractGSEs.C:115
static Color_t PbPbColor(Double_t c1, Double_t c2)
Definition: ExtractGSEs.C:24
int Int_t
Definition: External.C:63
TAxis * GetAxis(TCollection *c, const char *name, Bool_t verb=true)
Definition: ExtractGSEs.C:127
unsigned int UInt_t
Definition: External.C:33
void GetMinMax(Option_t *option, Double_t &ymin, Double_t &ymax) const
Definition: GraphSysErr.C:3985
TParameter< double > * GetParam(TCollection *c, const char *name, Bool_t verb=false)
Definition: ExtractGSEs.C:133
virtual void Print(Option_t *option="R") const
Definition: GraphSysErr.C:502
void SetSumFillStyle(Style_t style)
Definition: GraphSysErr.C:4485
TMultiGraph * GetMulti(Option_t *option="")
Definition: GraphSysErr.C:924
void SetSumOption(EDrawOption_t opt)
Definition: GraphSysErr.C:4449
TList * ExtractGSEs(const char *filename="forward_dndeta.root", Int_t rebin=5, Double_t eff=1, Bool_t raw=false, Bool_t cutEdges=false, const char *name="Forward")
Definition: ExtractGSEs.C:308
Definition: External.C:220
TFile * file
TList with histograms for a given trigger.
Int_t rebin
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
void SetSumFillColor(Color_t color)
Definition: GraphSysErr.C:4479
TObject * GetObject(TDirectory *d, const char *name, TClass *cls=0, Bool_t verb=true)
Definition: ExtractGSEs.C:42
const Int_t nbins
void SetCommonSumFillStyle(Style_t style)
Definition: GraphSysErr.C:4527
bool Bool_t
Definition: External.C:53
void ModOne(TObject *o, TString &sys, UInt_t sNN, TString trg, Color_t col)
Definition: ExtractGSEs.C:268
void AddPath(const TString &dir, Bool_t before=true)
Definition: ExtractGSEs.C:142
void SetSumLineColor(Color_t color)
Definition: GraphSysErr.C:4461
void SetCommonSumFillColor(Color_t color)
Definition: GraphSysErr.C:4521
Definition: External.C:196
TDirectoryFile * dir