AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GenerateEmpirical.C
Go to the documentation of this file.
1 #ifndef __CINT__
2 #include <TFile.h>
3 #include <TSystem.h>
4 #include <TROOT.h>
5 #include <TCollection.h>
6 #include <TString.h>
7 #include <TH1.h>
8 #include <TH2.h>
9 #include <TGraph.h>
10 #include <TGraphAsymmErrors.h>
11 #include <TGraphErrors.h>
12 #include <TClass.h>
13 #include <TMultiGraph.h>
14 #include <THStack.h>
15 #include <TError.h>
16 #else
17 class TCollection;
18 class TDirectory;
19 class TFile;
20 class TH1;
21 class TH2;
22 class TGraph;
23 class THStack;
24 class TAxis;
25 class TMultiGraph;
26 #endif
27 
37  static TFile* OpenFile(const TString& fname, Bool_t rw=false)
38  {
39  TFile* file = TFile::Open(fname, rw ? "RECREATE" : "READ");
40  if (!file) {
41  ::Error("OpenFile", "Couldn't open %s for %s",
42  fname.Data(), (rw ? "read/write" : "read"));
43  return 0;
44  }
45  return file;
46  }
56  static Bool_t CheckType(TObject* o, TClass* cl, Bool_t quiet=false)
57  {
58  if (!o || !cl) {
59  if (!quiet) ::Error("CheckType", "No object and/or class given");
60  return false;
61  }
62  if (!o->IsA()->InheritsFrom(cl)) {
63  if (!quiet) ::Error("CheckType", "Object %s is not a %s, but a %s",
64  o->GetName(), cl->GetName(), o->ClassName());
65  return false;
66  }
67  return true;
68  }
79  static TObject* GetObject(TDirectory* d, const TString& name,
80  TClass* cl=0, Bool_t quiet=false)
81  {
82  if (!d) {
83  if (!quiet) ::Error("GetObject", "No parent directory");
84  return 0;
85  }
86  TObject* o = d->Get(name.Data());
87  if (!o) {
88  if (!quiet) ::Error("GetObject", "didn't find the object %s in %s",
89  name.Data(), d->GetName());
90  if (!quiet) d->ls();
91  return 0;
92  }
93  if (cl && !CheckType(o, cl, quiet)) return 0;
94 
95  return o;
96  }
107  static TObject* GetObject(TCollection* d, const TString& name,
108  TClass* cl=0, Bool_t quiet=false)
109  {
110  if (!d) {
111  if (!quiet) ::Error("GetObject", "No parent directory");
112  return 0;
113  }
114  TObject* o = d->FindObject(name.Data());
115  if (!o) {
116  if (!quiet) ::Error("GetObject", "didn't find the object %s in %s",
117  name.Data(), d->GetName());
118  if (!quiet) d->ls();
119  return 0;
120  }
121  if (cl && !CheckType(o, cl, quiet)) return 0;
122 
123  return o;
124  }
134  static TCollection* GetCollection(TDirectory* d, const TString& name,
135  Bool_t quiet=false)
136  {
137  return static_cast<TCollection*>(GetObject(d, name, TCollection::Class(),
138  quiet));
139  }
149  static TCollection* GetCollection(TCollection* d, const TString& name,
150  Bool_t quiet=false)
151  {
152  return static_cast<TCollection*>(GetObject(d, name, TCollection::Class(),
153  quiet));
154  }
164  static TH1* GetH1(TDirectory* d, const TString& name, Bool_t quiet=false)
165  {
166  return static_cast<TH1*>(GetObject(d, name, TH1::Class(), quiet));
167  }
177  static TH1* GetH1(TCollection* d, const TString& name, Bool_t quiet=false)
178  {
179  return static_cast<TH1*>(GetObject(d, name, TH1::Class(), quiet));
180  }
190  static TGraph* GetGraph(TCollection* d, const TString& name, Bool_t quiet=false)
191  {
192  return static_cast<TGraph*>(GetObject(d, name, TGraph::Class(), quiet));
193  }
203  static TAxis* GetAxis(TCollection* d, const TString& name,
204  Bool_t quiet=false)
205  {
206  return static_cast<TAxis*>(GetObject(d, name, TAxis::Class(), quiet));
207  }
217  static TH1* RatioHG(const TH1* h, const TGraph* g)
218  {
219  if (!h || !g) return 0;
220 
221  TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
222  ret->Reset();
223  ret->SetDirectory(0);
224  Double_t xlow = g->GetX()[0];
225  Double_t xhigh = g->GetX()[g->GetN()-1];
226  if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
227 
228  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
229  Double_t c = h->GetBinContent(i);
230  if (c <= 0) continue;
231 
232  Double_t x = h->GetBinCenter(i);
233  if (x < xlow || x > xhigh) continue;
234 
235  Double_t f = g->Eval(x);
236  if (f <= 0) continue;
237 
238  ret->SetBinContent(i, c / f);
239  ret->SetBinError(i, h->GetBinError(i) / f);
240  }
241  return ret;
242  }
243  //==================================================================
248  : fRefLoaded(false)
249  {}
255  {
256 #if 0
257  if (fRefLoaded) return;
258  // --- Set the macro pathand load other data script --------------
259  TString savPath(gROOT->GetMacroPath());
260  gROOT->SetMacroPath(Form("%s:$(ALICE_PHYSICS)/PWGLF/FORWARD/analysis2",
261  gROOT->GetMacroPath()));
262  // Always recompile
263  if (!gROOT->GetClass("RefData"))
264  gROOT->LoadMacro("OtherData.C+");
265  gROOT->SetMacroPath(savPath);
266  fRefLoaded = true;
267 #endif
268  }
277  TMultiGraph* GetReference(UShort_t c1, UShort_t c2)
278  {
279  return 0;
280 #if 0
281 
282  if (!fRefLoaded) LoadReferences();
283 
284  TMultiGraph* other = 0;
285  UShort_t sys = 2; // (fSysString ? fSysString->GetUniqueID() : 0);
286  UShort_t trg = 0; // (fTrigString ? fTrigString->GetUniqueID() : 0);
287  UShort_t snn = 2760; // (fSNNString ? fSNNString->GetUniqueID() : 0);
288  Long_t ret =
289  gROOT->ProcessLine(Form("RefData::GetData(%d,%d,%d,%d,%d,%d);",
290  sys,snn,trg,c1,c2,0xF));
291  if (!ret) {
292  Warning("GetReference",
293  "No other data for sys=%d trg=%d sNN=%d c=%3d%%-%3d%%",
294  sys, trg, snn, c1, c2);
295  return 0;
296  }
297 
298  other = reinterpret_cast<TMultiGraph*>(ret);
299  return other;
300 #endif
301  }
307  void Run(const TString& fileName)
308  {
309  TFile* file = OpenFile(fileName, false);
310  if (!file) return;
311 
312  TFile* out = OpenFile("empirical.root", true);
313  if (!out) return;
314 
315  ProcessComponent(file, out, "Forward");
316  ProcessComponent(file, out, "Central");
317 
318  out->ls();
319  out->Close();
320  file->Close();
321  }
329  void ProcessComponent(TDirectory* d, TDirectory* out, const TString& name)
330  {
331  TDirectory* store = out->mkdir(name);
332  store->cd();
333 
334  TString resName(Form("%sdNdetaResults", name.Data()));
335  TCollection* results = GetCollection(d, resName);
336  if (!results) return;
337 
338  TAxis* centAxis = GetAxis(results, "centAxis");
339  if (!centAxis) return;
340 
341  THStack* corrs = new THStack("empirical", "Empirical correction");
342  // TMultiGraph* corrs = new TMultiGraph("empirical");
343  corrs->SetTitle(Form("Empirical corrections for %s", name.Data()));
344  for (Int_t i=1; i<=centAxis->GetNbins(); i++) {
345  UShort_t c1 = UShort_t(centAxis->GetBinLowEdge(i));
346  UShort_t c2 = UShort_t(centAxis->GetBinUpEdge(i));
347 
348  ProcessCent(c1, c2, results, name, corrs, store);
349  }
350 
351  if (!corrs->GetHists() ||
352  corrs->GetHists()->GetEntries() <= 0) return;
353 
354  store->cd();
355 
356  TH1* sum = static_cast<TH1*>(corrs->GetHists()->At(0)->Clone("mean"));
357  sum->SetTitle("mean");
358  sum->Reset();
359  sum->SetMarkerColor(kBlack);
360  sum->SetDirectory(0);
361  Int_t nHist = corrs->GetHists()->GetEntries();
362  for (Int_t i = 0; i < nHist; i++)
363  sum->Add(static_cast<TH1*>(corrs->GetHists()->At(i)));
364  sum->Scale(1. / nHist);
365  sum->Write("default");
366 
367  THStack* ratios = new THStack("ratios", "Ratio to mean");
368  for (Int_t i = 0; i < nHist; i++) {
369  TH1* h = static_cast<TH1*>(corrs->GetHists()->At(i));
370  TH1* r = static_cast<TH1*>(h->Clone());
371  r->SetDirectory(0);
372  r->Divide(sum);
373  ratios->Add(r);
374  }
375 
376 
377  TCanvas* c = new TCanvas(name);
378  c->SetRightMargin(0.05);
379  c->SetTopMargin(0.05);
380  c->Divide(1,2,0,0);
381  TVirtualPad* p = c->cd(1);
382  corrs->Add(sum);
383  corrs->DrawClone("nostack");
384  p->BuildLegend(0.2, 0.05, 0.4, 0.4);
385  corrs->Write();
386 
387  p = c->cd(2);
388  THStack* cpy = static_cast<THStack*>(ratios->DrawClone("nostack"));
389  cpy->SetMinimum(0.95);
390  cpy->SetMaximum(1.05);
391  p->BuildLegend(0.35,0.15,0.55,0.4);
392  ratios->Write();
393 
394  c->Print(Form("%sEmpirical.pdf", name.Data()));
395  }
407  const TString& name, THStack* s, TDirectory* d)
408  {
409  // @todo reimplement to get from GSE's
410  TMultiGraph* others = GetReference(c1, c2);
411  if (!others) return;
412 
413  TGraph* alice = GetGraph(others->GetListOfGraphs(),
414  Form("alice_pbpb2760"));
415  if (!alice) return;
416 
417  TString folderName(Form("cent%03d_%03d", c1, c2));
418  TCollection* centBin = GetCollection(c, folderName);
419  if (!centBin) return;
420 
421  TH1* ana = GetH1(centBin, Form("dndeta%s", name.Data()));
422  if (!ana) return;
423 
424  TH1* ratio = RatioHG(ana, alice);
425  if (!ratio) return;
426 
427  ::Info("", "Adding %s/%s", name.Data(), folderName.Data());
428  ratio->SetName(folderName);
429  ratio->SetTitle(Form("%3d%% - %3d%%", c1, c2));
430 
431  // if (c1 == 0) ratio->Write("default");
432 
433  TDirectory* store = d->mkdir(folderName);
434  store->cd();
435  ana->Write();
436  alice->Write();
437  d->cd();
438 
439  s->Add(ratio);
440  }
441 
443 };
444 
445 void GenerateEmpirical(const Char_t* fileName="PbPb_2760_dndeta_nosec_CENT_20140513_1349/forward_dndeta.root")
446 {
447  EmpiricalMaker m;
448  m.Run(fileName);
449 }
void ana(Int_t mode=mGRID)
Definition: ana.C:87
double Double_t
Definition: External.C:58
static TH1 * RatioHG(const TH1 *h, const TGraph *g)
TString fileName
char Char_t
Definition: External.C:18
static Bool_t CheckType(TObject *o, TClass *cl, Bool_t quiet=false)
void Run(const TString &fileName)
static TAxis * GetAxis(TCollection *d, const TString &name, Bool_t quiet=false)
void GenerateEmpirical(const Char_t *fileName="PbPb_2760_dndeta_nosec_CENT_20140513_1349/forward_dndeta.root")
TCanvas * c
Definition: TestFitELoss.C:172
TMultiGraph * GetReference(UShort_t c1, UShort_t c2)
static TH1 * GetH1(TDirectory *d, const TString &name, Bool_t quiet=false)
static TObject * GetObject(TCollection *d, const TString &name, TClass *cl=0, Bool_t quiet=false)
static TFile * OpenFile(const TString &fname, Bool_t rw=false)
int Int_t
Definition: External.C:63
void ProcessComponent(TDirectory *d, TDirectory *out, const TString &name)
void ProcessCent(UShort_t c1, UShort_t c2, TCollection *c, const TString &name, THStack *s, TDirectory *d)
static TCollection * GetCollection(TDirectory *d, const TString &name, Bool_t quiet=false)
static TObject * GetObject(TDirectory *d, const TString &name, TClass *cl=0, Bool_t quiet=false)
static TH1 * GetH1(TCollection *d, const TString &name, Bool_t quiet=false)
static TGraph * GetGraph(TCollection *d, const TString &name, Bool_t quiet=false)
Definition: External.C:220
TFile * file
TList with histograms for a given trigger.
unsigned short UShort_t
Definition: External.C:28
bool Bool_t
Definition: External.C:53
static TCollection * GetCollection(TCollection *d, const TString &name, Bool_t quiet=false)
Definition: External.C:196