AliPhysics  323987d (323987d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 TH1*
269 DrawOne(TObject* o, Option_t* opt, Double_t& min, Double_t& max)
270 {
271  GraphSysErr* g = static_cast<GraphSysErr*>(o);
272  g->Draw(opt);
273  Double_t mn, mx;
274  g->GetMinMax(opt, mn, mx);
275  min = TMath::Min(min, mn);
276  max = TMath::Max(max, mx);
277  TMultiGraph* mg = g->GetMulti();
278  // Info("DrawOne", "Got Multigraph %p", mg);
279  TH1* hist = (mg ? mg->GetHistogram() : 0);
280  // Info("DrawOne", "Got frame histogram %p", hist);
281  return hist;
282 }
283 
284 
285 //____________________________________________________________________
286 TList*
287 ExtractGSEs(const char* filename="forward_dndeta.root",
288  Int_t rebin=5,
289  Double_t eff=1,
290  Bool_t raw=false,
291  Bool_t cutEdges=false,
292  const char* name="Forward")
293 {
294  if (!gROOT->GetClass("GraphSysErr")) {
295  AddPath("$HOME/GraphSysErr");
296  AddRelPath("gse",false);
297 
298  gROOT->LoadMacro("GraphSysErr.C+g");
299  }
300  if (!gROOT->GetClass("SysErrorAdder")) {
301  AddRelPath("dndeta", true);
302  gROOT->LoadMacro("SysErrorAdder.C+g");
303  }
304 
305  TFile* file = TFile::Open(filename,"READ");
306  if (!file) {
307  Warning("ExtractGSEs", "Failed to open %s", filename);
308  return 0;
309  }
310 
311  TCollection* results = GetCollection(file, Form("%sdNdetaResults",name));
312  if (!results) return 0;
313 
314 
315  TH1* empH = GetH1 (results, "empirical");
316  Bool_t emp = (empH != 0 && !raw);
317  TObject* osNN = GetObject(results, "sNN");
318  TObject* osys = GetObject(results, "sys");
319  TObject* otrg = GetObject(results, "trigger");
320  TObject* ocentM = GetObject(results, "centEstimator");
321  TAxis* centA = GetAxis (results, "centAxis");
322  TString sys = osys->GetTitle();
323  UShort_t sNN = osNN->GetUniqueID();
324  TString mth = (ocentM && centA ? ocentM->GetTitle() : "");
325  if (mth.EqualTo("none",TString::kIgnoreCase) ||
326  (centA && centA->GetNbins() < 1)) mth = "";
327  TString trg = (!mth.IsNull() && centA ? "CENT" : otrg->GetTitle());
328  TString hname = Form("dndeta%s%s",name,emp?"Emp":"");
329 
330  Printf("Read settings\n"
331  " Empirical: %s (%p)\n"
332  " System: %s (%p)\n"
333  " Energy: %d (%s - %p)\n"
334  " Trigger: %s (%s - %p)\n"
335  " Centrality: %s (%p)\n"
336  " Axis: %p",
337  (emp ? "yes" : "no"), empH,
338  sys.Data(), osys,
339  sNN, (osNN ? osNN->GetTitle() : ""), osNN,
340  trg.Data(), otrg->GetTitle(), otrg, mth.Data(),
341  ocentM, centA);
342  if (centA) Printf(" %d bins between %f and %f",
343  centA->GetNbins(), centA->GetXmin(), centA->GetXmax());
344  SysErrorAdder* adder = SysErrorAdder::Create(trg,sys,sNN,mth);
345  if (!adder) {
346  Warning("ExtractGSEs", "Failed to make adder");
347  return 0;
348  }
349 
350 
351  TCanvas* c = new TCanvas("c", "C");
352  c->SetTopMargin(0.01);
353  c->SetRightMargin(0.20);
354  TLegend* l = new TLegend(0.8,0.1,0.99,0.99,
355  Form("%s @ %s (%s)",
356  sys.Data(), osNN->GetTitle(), trg.Data()));
357  l->SetFillColor(0);
358  l->SetFillStyle(0);
359  l->SetBorderSize(0);
360 
361 
362  TList* ret = new TList;
363  Bool_t first = true;
364  TH1* frame = 0;
365  Double_t min = 100000, max = -1000000;
366  if (!centA || centA->GetNbins() < 1 || mth.IsNull()) {
367  Info("ExtractGSEs", "Doing pp-like extraction"
368  " all bin, %s, rebin=%d, eff=%f, %p, %s",
369  hname.Data(), rebin, eff, adder,
370  (cutEdges ? "cut edges" : "edges stay"));
371  TObject* all = DoOne(results,"all",hname,rebin,eff,adder,cutEdges,l);
372  if (!all) {
373  Warning("ExtractGSEs", "Nothing returned from DoOne(\"all\"...)");
374  return 0;
375  }
376  ret->Add(all);
377  frame = DrawOne(all, "SUM QUAD AXIS", min, max);
378  }
379  else {
380  for (Int_t i = 1; i <= centA->GetNbins(); i++) {
381  Double_t low = centA->GetBinLowEdge(i);
382  Double_t high = centA->GetBinUpEdge(i);
383  TString dir;
384  dir.Form("cent%03dd%02d_%03dd%02d",
385  Int_t(low), Int_t(low *100)%100,
386  Int_t(high), Int_t(high*100)%100);
387  TObject* g = DoOne(results, dir, hname, rebin, eff, adder, cutEdges,
388  (first ? l : 0));
389  if (!g) continue;
390  ret->Add(g);
391 
392  GraphSysErr* gse = static_cast<GraphSysErr*>(g);
393  Color_t col = PbPbColor(low, high);
394  gse->SetMarkerColor(col);
395  gse->SetSumFillColor(kWhite);
396  gse->SetSumFillStyle(0);
397  gse->SetSumLineColor(col);
398  gse->SetCommonSumFillColor(col);
399  gse->SetCommonSumFillStyle(1001);
400  gse->SetCommonSumLineColor(col);
401 
402  if (first) frame = DrawOne(g, "SUM QUAD AXIS", min, max);
403  else DrawOne(g, "SUM QUAD", min, max);
404  TLegendEntry* e = l->AddEntry("", Form("%3d-%3d%%",
405  int(low), int(high)),"F");
406  e->SetFillColor(col);
407  e->SetFillStyle(1001);
408  e->SetLineColor(col);
409  e->SetMarkerColor(col);
410  first = false;
411  }
412  }
413  if (!frame) {
414  Warning("ExtractGSEs", "No frame given");
415  }
416  else {
417  frame->SetMinimum(0.9*min);
418  frame->SetMaximum(1.1*max);
419  }
420  l->Draw();
421 
422  TString outName;
423  if (mth.EqualTo("V0M", TString::kIgnoreCase) &&
424  sys.EqualTo("PbPb",TString::kIgnoreCase))
425  mth = "";
426  if (raw) mth = "RAW";
427  outName.Form("%s_%05d_%s%s.root",sys.Data(),sNN,trg.Data(),mth.Data());
428  TFile* out = TFile::Open(outName,"RECREATE");
429  ret->Write("container",TObject::kSingleKey);
430  out->Write();
431 
432  Info("ExtractGSEs", "Written to %s", outName.Data());
433  return ret;
434 }
435 
TH1 * DrawOne(TObject *o, Option_t *opt, Double_t &min, Double_t &max)
Definition: ExtractGSEs.C:269
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
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
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
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:287
Definition: External.C:220
TFile * file
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 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