AliPhysics  608b256 (608b256)
DrawUA5Ratios.C
Go to the documentation of this file.
1 //____________________________________________________________________
12 TObject*
13 GetObject(const TObject* l, const char* name)
14 {
15  if (!l->IsA()->InheritsFrom(TCollection::Class())) {
16  Error("GetObject", "passed parent %s not a TCollection, but %s",
17  l->GetName(), l->IsA()->GetName());
18  return 0;
19  }
20  TObject* o = static_cast<const TCollection*>(l)->FindObject(name);
21  if (!o) {
22  Error("GetObject", "No object '%s' found in '%s'", name, l->GetName());
23  return 0;
24  }
25  return o;
26 }
27 
28 //____________________________________________________________________
41 TH1D*
42 GetHist(TDirectory* dir,
43  const char* which,
44  const char* sub,
45  const char* hname)
46 {
47  if (!dir) return 0;
48  TList* parent = static_cast<TList*>(dir->Get(which));
49  if (!parent) {
50  Error("GetHist", "List '%s' not found in '%s'", which, dir->GetName());
51  return 0;
52  }
53  TList* child = static_cast<TList*>(parent->FindObject(sub));
54  if (!child) {
55  Error("GetHist", "List '%s' not found in '%s'", sub, parent->GetName());
56  return 0;
57  }
58  TObject* o = GetObject(child,hname);
59  if (!o) return 0;
60  return static_cast<TH1D*>(o);
61 }
62 
63 
64 //____________________________________________________________________
78 TH1D*
79 GetHist(TDirectory* dir,
80  const char* which,
82  const char* sub="all")
83 {
84  TString name(Form("dndeta%s", which));
85  if (rebin > 1) name.Append(Form("_rebin%02d", rebin));
86  return GetHist(dir, Form("%sResults", which), sub, name);
87 }
88 
89 //____________________________________________________________________
102 TH1*
103 Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
104 {
105  TH1* tmp = static_cast<TH1*>(fwd->Clone("dndetaMerged"));
106  // tmp->SetMarkerStyle(28);
107  // tmp->SetMarkerColor(kBlack);
108  tmp->SetDirectory(0);
109  xlow = 100;
110  xhigh = -100;
111  for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
112  Double_t cc = cen->GetBinContent(i);
113  Double_t cf = fwd->GetBinContent(i);
114  Double_t ec = cen->GetBinError(i);
115  Double_t ef = fwd->GetBinError(i);
116  Double_t nc = cf;
117  Double_t ne = ef;
118  if (cc < 0.001 && cf < 0.01) continue;
119  xlow = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
120  xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
121  if (cc > 0.001) {
122  nc = cc;
123  ne = ec;
124  if (cf > 0.001) {
125  nc = (cf + cc) / 2;
126  ne = TMath::Sqrt(ec*ec + ef*ef);
127  }
128  }
129  tmp->SetBinContent(i, nc);
130  tmp->SetBinError(i, ne);
131  }
132  return tmp;
133 }
134 
135 //____________________________________________________________________
152 {
153  Double_t x = xp[0];
154  Double_t a1 = pp[0];
155  Double_t a2 = pp[1];
156  Double_t s1 = pp[2];
157  Double_t s2 = pp[3];
158  return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
159 }
160 //____________________________________________________________________
176 {
177  Double_t denom = myFunc(xp, &(pp[4]));
178  if (TMath::Abs(denom) < 1.e-6) return 0;
179  return myFunc(xp, pp) / denom;
180 }
181 
182 //____________________________________________________________________
194 TF1*
195 FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
196 {
197  TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
198  tmp->Fit(tmpf, "N", "");
199  tmp->SetDirectory(0);
200 
201  TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
202  fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
203  fit->SetParameters(tmpf->GetParameter(0),
204  .2,
205  tmpf->GetParameter(2),
206  tmpf->GetParameter(2)/4);
207  fit->SetParLimits(3, 0, 100);
208  fit->SetParLimits(4, 0, 100);
209  tmp->Fit(fit,"0W","");
210  return fit;
211 }
212 
213 //____________________________________________________________________
222 void
223 MakeSysError(TH1* tmp, TF1* fit)
224 {
225  for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
226  Double_t tc = tmp->GetBinContent(i);
227  if (tc < 0.01) continue;
228  Double_t x = tmp->GetXaxis()->GetBinCenter(i);
229  Double_t y = fit->Eval(x);
230  tmp->SetBinContent(i, y);
231  tmp->SetBinError(i,sysErr*y);
232  }
233  return tmp;
234 }
235 
236 //____________________________________________________________________
246 TH1*
248 {
249  Int_t nBins = g->GetN();
250  TArrayF bins(nBins+1);
251  TArrayF y(nBins);
252  TArrayF ey(nBins);
253  Double_t dx = 0;
254  Double_t xmin = 10000;
255  Double_t xmax = -10000;
256  for (Int_t i = 0; i < nBins; i++) {
257  Double_t x = g->GetX()[i];
258  Double_t exl = g->GetEXlow()[i];
259  Double_t exh = g->GetEXhigh()[i];
260  xmin = TMath::Min(x-exl, xmin);
261  xmax = TMath::Max(x+exh, xmax);
262  bins.fArray[i] = x-exl;
263  bins.fArray[i+1] = x+exh;
264  Double_t dxi = exh+exl;
265  if (dxi == 0 && i != 0) dxi = bins.fArray[i]-bins.fArray[i-1];
266  if (dx == 0) dx = dxi;
267  else if (dxi != dx) dx = 0;
268 
269  y.fArray[i] = g->GetY()[i];
270  ey.fArray[i] = TMath::Max(g->GetEYlow()[i],g->GetEYhigh()[i]);
271 
272  }
273  TString name(g->GetName());
274  TString title(g->GetTitle());
275  TH1D* h = 0;
276  if (dx != 0) {
277  h = new TH1D(name.Data(), title.Data(), nBins,
278  bins[0]-dx/2, bins[nBins]+dx/2);
279  }
280  else {
281  h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
282  }
283  for (Int_t i = 1; i <= nBins; i++) {
284  h->SetBinContent(i, y.fArray[i-1]);
285  h->SetBinError(i, ey.fArray[i-1]);
286  }
287  h->SetMarkerStyle(g->GetMarkerStyle());
288  h->SetMarkerColor(g->GetMarkerColor());
289  h->SetMarkerSize(g->GetMarkerSize());
290  h->SetDirectory(0);
291 
292  return h;
293 }
294 
295 //____________________________________________________________________
307 TH1*
308 Ratio(TH1* h, TF1* f, const char* title)
309 {
310  TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_%s",
311  h->GetName(),
312  f->GetName())));
313  ret->SetDirectory(0);
314  if (title) ret->SetTitle(title);
315  else ret->SetTitle(Form("%s (data) / %s",
316  h->GetTitle(),f->GetTitle()));
317  ret->Reset();
318  for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
319  Double_t cc = h->GetBinContent(i);
320  if (cc < 0.01) {
321  ret->SetBinContent(i,0);
322  ret->SetBinError(i,0);
323  continue;
324  }
325  Double_t xx = h->GetBinCenter(i);
326  Double_t ee = h->GetBinError(i);
327  Double_t ff = f->Eval(xx);
328  Double_t yy = cc / ff;
329  Double_t ey = ee / ff;
330  ret->SetBinContent(i, yy);
331  ret->SetBinError(i, ey);
332  }
333  return ret;
334 }
335 
336 //____________________________________________________________________
351 TH1D*
352 GetUA5Data(UShort_t type, TH1*& p, TH1*& n,
353  Double_t& xlow, Double_t& xhigh)
354 {
355  gROOT->SetMacroPath(Form(".:$ALICE_PHYSICS.trunk/PWGLF/FORWARD/analysis2/:%s",
356  gROOT->GetMacroPath()));
357  gROOT->LoadMacro("OtherData.C");
358 
359  p = 0;
360  n = 0;
361  TGraphAsymmErrors* gp = GetSingle(UA5, 1, 900, type, 0, 0);
362  TGraphAsymmErrors* gn = GetSingle(UA5+10, 1, 900, type, 0, 0);
363  if (!gp || !gn) return 0;
364 
365  p = Graph2Hist(gp);
366  n = Graph2Hist(gn);
367 
368  Int_t nn = p->GetNbinsX();
369  xlow = n->GetXaxis()->GetXmin();
370  xhigh = p->GetXaxis()->GetXmax();
371  TH1D* ret = new TH1D("ua5", "UA5", 2*nn, xlow, xhigh);
372  ret->SetDirectory(0);
373  ret->SetMarkerColor(p->GetMarkerColor());
374  ret->SetMarkerStyle(p->GetMarkerStyle());
375 
376  for (Int_t i = 1; i <= nn; i++) {
377  ret->SetBinContent(nn+i, p->GetBinContent(i));
378  ret->SetBinContent( i, n->GetBinContent(i));
379  ret->SetBinError(nn+i, p->GetBinError(i));
380  ret->SetBinError( i, n->GetBinError(i));
381  }
382  return ret;
383 }
384 
385 
386 //____________________________________________________________________
396 void
397 DrawUA5Ratios(const char* fname="forward_dndeta.root", UShort_t rebin=5)
398 {
399  TFile* forward_dndeta = TFile::Open(fname, "READ");
400  if (!forward_dndeta) {
401  Error("DrawUA5Ratios", "%s not found", fname);
402  return;
403  }
404 
405  UShort_t type = 1;
406  TH1D* forward = GetHist(forward_dndeta, "Forward", rebin);
407  TH1D* central = GetHist(forward_dndeta, "Central", rebin);
408 
409  TObject* sys = GetObject(forward_dndeta->Get("ForwardResults"), "sys");
410  TObject* sNN = GetObject(forward_dndeta->Get("ForwardResults"), "sNN");
411  TObject* trg = GetObject(forward_dndeta->Get("ForwardResults"), "trigger");
412  if (!sys || !sNN || !trg) return;
413  if (sys->GetUniqueID() != 1) {
414  Error("DrawUA5Ratios", "Comparison only valid for pp, not %s",
415  sys->GetTitle());
416  return;
417  }
418  if (sNN->GetUniqueID() != 900) {
419  Error("DrawUA5Ratios", "Comparison only valid for 900GeV, not %dGeV",
420  sNN->GetUniqueID());
421  return;
422  }
423  if (trg->GetUniqueID() != 1 &&
424  trg->GetUniqueID() != 4) {
425  Error("DrawUA5Ratios",
426  "Comparison only valid for INEL or NSD, not %s (%d)",
427  trg->GetTitle(), trg->GetUniqueID());
428  return;
429  }
430 
431 
432  Double_t alilow, alihigh;
433  TH1D* ali = Merge(central, forward, alilow, alihigh);
434  TF1* alifit = FitMerged(ali, alilow, alihigh);
435  ali->SetTitle("Forward+Central");
436  alifit->SetLineColor(kRed+1);
437  alifit->SetLineStyle(2);
438  alifit->SetName("alifit");
439  alifit->SetTitle("Fit to Forward+Central");
440 
441  Double_t ua5low, ua5high;
442  TH1* ua5p, *ua5n;
443  TH1D* ua5 = GetUA5Data(trg->GetUniqueID(), ua5p, ua5n, ua5low, ua5high);
444  TF1* ua5fit = FitMerged(ua5, ua5low, ua5high);
445  ua5fit->SetLineColor(kBlue+1);
446  ua5fit->SetLineStyle(3);
447  ua5fit->SetName("ua5fit");
448  ua5fit->SetTitle("Fit to UA5");
449 
450  gStyle->SetOptTitle(0);
451  TCanvas* c = new TCanvas("c", "C", 900, 900);
452  c->SetFillColor(0);
453  c->SetFillStyle(0);
454  c->SetBorderMode(0);
455  c->SetBorderSize(0);
456 
457  Double_t yd = .3;
458 
459  TPad* p1 = new TPad("p1", "p1", 0, yd, 1, 1, 0, 0, 0);
460  p1->SetBorderSize(0);
461  p1->SetBorderMode(0);
462  p1->SetFillColor(0);
463  p1->SetFillStyle(0);
464  p1->SetRightMargin(0.02);
465  p1->SetTopMargin(0.02);
466  p1->SetBottomMargin(0.00);
467  p1->SetGridx();
468  p1->Draw();
469  p1->cd();
470 
471  THStack* all = new THStack("all", "All");
472  all->Add(ua5p);
473  all->Add(ua5n);
474  // all->Add(ua5);
475  all->Add(forward);
476  all->Add(central);
477  // all->Add(merged);
478  all->Draw("nostack");
479  all->SetMinimum(-.07);
480  all->GetXaxis()->SetTitleFont(132);
481  all->GetYaxis()->SetTitleFont(132);
482  all->GetXaxis()->SetLabelFont(132);
483  all->GetYaxis()->SetLabelFont(132);
484  all->GetYaxis()->SetDecimals();
485  p1->Clear();
486  all->Draw("nostack");
487  // ua5p->Draw("same p");
488  // ua5m->Draw("same p");
489  alifit->Draw("same");
490  ua5fit->Draw("same");
491 
492  TLegend* l = new TLegend(.2, .1, .8, .5,
493  Form("pp @ #sqrt{s}=900GeV, %s",trg->GetTitle()));
494  l->AddEntry(ua5, Form("U: %s", ua5->GetTitle()), "lp");
495  l->AddEntry(forward, "F: Forward", "lp");
496  l->AddEntry(central, "C: Central", "lp");
497  l->AddEntry(alifit,
498  Form("f: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
499  "%4.2fe^{(x/%4.2f)^{2}}#right]",
500  alifit->GetTitle(),
501  alifit->GetParameter(0),
502  alifit->GetParameter(2),
503  alifit->GetParameter(1),
504  alifit->GetParameter(3)), "l");
505  l->AddEntry(ua5fit,
506  Form("u: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
507  "%4.2fe^{(x/%4.2f)^{2}}#right]",
508  ua5fit->GetTitle(),
509  ua5fit->GetParameter(0),
510  ua5fit->GetParameter(2),
511  ua5fit->GetParameter(1),
512  ua5fit->GetParameter(3)), "l");
513  l->SetTextFont(132);
514  l->SetBorderSize(0);
515  l->SetFillColor(0);
516  l->SetFillStyle(0);
517  l->Draw();
518 
519  c->cd();
520  TPad* p2 = new TPad("p2", "p2", 0, 0, 1, yd, 0, 0, 0);
521  p2->SetBorderSize(0);
522  p2->SetBorderMode(0);
523  p2->SetFillColor(0);
524  p2->SetFillStyle(0);
525  p2->SetRightMargin(0.02);
526  p2->SetTopMargin(0.00);
527  p2->SetBottomMargin(0.15);
528  p2->SetGridx();
529  p2->Draw();
530  p2->cd();
531 
532  THStack* ratios = new THStack("Ratios", "Ratios");
533  TH1* ua5ali = Ratio(ua5, alifit, 0);
534  TH1* aliua5 = Ratio(ali, ua5fit, 0);
535  ratios->Add(ua5ali);
536  ratios->Add(aliua5);
537  ratios->Draw("nostack");
538  ratios->SetMinimum(0.4);
539  ratios->GetYaxis()->SetTitleSize(2*ratios->GetYaxis()->GetTitleSize());
540  ratios->GetYaxis()->SetLabelSize(2*ratios->GetYaxis()->GetLabelSize());
541  ratios->GetYaxis()->SetNdivisions(508);
542  ratios->GetXaxis()->SetTitleSize(2*ratios->GetXaxis()->GetTitleSize());
543  ratios->GetXaxis()->SetLabelSize(2*ratios->GetXaxis()->GetLabelSize());
544  ratios->GetXaxis()->SetNdivisions(510);
545  ratios->GetXaxis()->SetTitle("#eta");
546  ratios->GetXaxis()->SetTitleFont(132);
547  ratios->GetYaxis()->SetTitleFont(132);
548  ratios->GetXaxis()->SetLabelFont(132);
549  ratios->GetYaxis()->SetLabelFont(132);
550  ratios->GetYaxis()->SetDecimals();
551  p2->Clear();
552 
553  TGraphErrors* sysErr = new TGraphErrors(2);
554  sysErr->SetPoint(0, all->GetHistogram()->GetXaxis()->GetXmin(),1);
555  sysErr->SetPoint(1, all->GetHistogram()->GetXaxis()->GetXmax(),1);
556  sysErr->SetPointError(0,0,.07);
557  sysErr->SetPointError(1,0,.07);
558  sysErr->SetTitle("Systematic error on Forward+Central");
559  sysErr->SetFillColor(kYellow+1);
560  sysErr->SetFillStyle(3001);
561  sysErr->SetHistogram(ratios->GetHistogram());
562  sysErr->DrawClone("ael3");
563  ratios->DrawClone("nostack same");
564 
565  TF1* fitfit = new TF1("fitfit", myRatio, alilow, alihigh, 8);
566  fitfit->SetParameters(ua5fit->GetParameter(0),
567  ua5fit->GetParameter(1),
568  ua5fit->GetParameter(2),
569  ua5fit->GetParameter(3),
570  alifit->GetParameter(0),
571  alifit->GetParameter(1),
572  alifit->GetParameter(2),
573  alifit->GetParameter(3));
574  fitfit->Draw("same");
575 
576  TLegend* ll = new TLegend(.3,.15, .7, .6);
577  ll->AddEntry(sysErr, sysErr->GetTitle(), "f");
578  ll->AddEntry(ua5ali, ua5ali->GetTitle(), "lp");
579  ll->AddEntry(aliua5, aliua5->GetTitle(), "lp");
580  ll->AddEntry(fitfit, "UA5 (fit) / Forward+Central (fit)", "lp");
581  ll->SetTextFont(132);
582  ll->SetBorderSize(0);
583  ll->SetFillColor(0);
584  ll->SetFillStyle(0);
585  ll->Draw();
586 
587 
588  c->SaveAs(Form("ua5_ratios_%s_r%02d.png", trg->GetTitle(), rebin));
589  gROOT->GetInterpreter()->UnloadFile("OtherData.C");
590 
591 }
592 
593 //____________________________________________________________________
594 //
595 // EOF
596 //
Double_t myRatio(Double_t *xp, Double_t *pp)
const Color_t cc[]
Definition: DrawKs.C:1
TH1D * GetUA5Data(UShort_t type, TH1 *&p, TH1 *&n, Double_t &xlow, Double_t &xhigh)
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
TCanvas * c
Definition: TestFitELoss.C:172
TH1 * Graph2Hist(const TGraphAsymmErrors *g)
int Int_t
Definition: External.C:63
TObject * GetObject(const TObject *l, const char *name)
Definition: DrawUA5Ratios.C:13
TH1 * Ratio(TH1 *h, TF1 *f, const char *title)
void MakeSysError(TH1 *tmp, TF1 *fit)
Definition: External.C:212
TH1 * Merge(const TH1 *cen, const TH1 *fwd, Double_t &xlow, Double_t &xhigh)
TObject * FindObject(int bin, const char *nameH, const TList *lst, Bool_t normPerEvent=kTRUE)
Double_t myFunc(Double_t *xp, Double_t *pp)
TF1 * FitMerged(TH1 *tmp, Double_t xlow, Double_t xhigh)
const char * fwd
Int_t rebin
TH1D * GetHist(TDirectory *dir, const char *which, const char *sub, const char *hname)
Definition: DrawUA5Ratios.C:42
unsigned short UShort_t
Definition: External.C:28
TList * ef
Definition: TestFitELoss.C:145
void DrawUA5Ratios(const char *fname="forward_dndeta.root", UShort_t rebin=5)
Definition: External.C:196
TDirectoryFile * dir