1 //____________________________________________________________________
3 * Get an object with specified name from TCollection @a l
6 * @param name Name of object to retrieve
8 * @return Object, or null
10 * @ingroup pwglf_forward_scripts
13 GetObject(const TObject* l, const char* name)
15 if (!l->IsA()->InheritsFrom(TCollection::Class())) {
16 Error("GetObject", "passed parent %s not a TCollection, but %s",
17 l->GetName(), l->IsA()->GetName());
20 TObject* o = static_cast<const TCollection*>(l)->FindObject(name);
22 Error("GetObject", "No object '%s' found in '%s'", name, l->GetName());
28 //____________________________________________________________________
30 * Get histogram from @a directory/which/sub/hname
32 * @param dir Directory
33 * @param which Name of parent list
34 * @param sub Name of sub-list
35 * @param hname Name of histogram
37 * @return Pointer to histogram or null
39 * @ingroup pwglf_forward_scripts
42 GetHist(TDirectory* dir,
48 TList* parent = static_cast<TList*>(dir->Get(which));
50 Error("GetHist", "List '%s' not found in '%s'", which, dir->GetName());
53 TList* child = static_cast<TList*>(parent->FindObject(sub));
55 Error("GetHist", "List '%s' not found in '%s'", sub, parent->GetName());
58 TObject* o = GetObject(child,hname);
60 return static_cast<TH1D*>(o);
64 //____________________________________________________________________
67 * <i>dir/which</i>Results<i>/sub/</i>dndeta<i>which</i>[_rebin<i>rebin</i>]
69 * @param dir Directory
71 * @param rebin Optional rebinning
72 * @param sub Sub-list name
74 * @return Histogram pointer or null
76 * @ingroup pwglf_forward_scripts
79 GetHist(TDirectory* dir,
82 const char* sub="all")
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);
89 //____________________________________________________________________
91 * Merge two histograms into one
93 * @param cen Central part
94 * @param fwd Forward part
95 * @param xlow On return, lower eta bound
96 * @param xhigh On return, upper eta bound
98 * @return Newly allocated histogram or null
100 * @ingroup pwglf_forward_scripts
103 Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
105 TH1* tmp = static_cast<TH1*>(fwd->Clone("dndetaMerged"));
106 // tmp->SetMarkerStyle(28);
107 // tmp->SetMarkerColor(kBlack);
108 tmp->SetDirectory(0);
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);
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);
126 ne = TMath::Sqrt(ec*ec + ef*ef);
129 tmp->SetBinContent(i, nc);
130 tmp->SetBinError(i, ne);
135 //____________________________________________________________________
137 * Function to calculate
139 * g(x;A_1,A_2,\sigma_1,\sigma_2) =
140 * A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} -
141 * A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
144 * @param xp Pointer to x array
145 * @param pp Pointer to parameter array
147 * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
149 * @ingroup pwglf_forward_scripts
151 Double_t myFunc(Double_t* xp, Double_t* pp)
158 return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
160 //____________________________________________________________________
164 * r(x) = \frac{g(x;A_1,A_2,\sigma_1,\sigma_2)}{
165 * g(x;A_1',A_2',\sigma'_1,\sigma'_2)}
168 * @param xp Pointer to X array
169 * @param pp Pointer to parameter array (8 entries)
173 * @ingroup pwglf_forward_scripts
175 Double_t myRatio(Double_t* xp, Double_t* pp)
177 Double_t denom = myFunc(xp, &(pp[4]));
178 if (TMath::Abs(denom) < 1.e-6) return 0;
179 return myFunc(xp, pp) / denom;
182 //____________________________________________________________________
184 * Fit @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
186 * @param tmp Histogram
187 * @param xlow Lower x bound
188 * @param xhigh Upper x bound
190 * @return Fitted function
192 * @ingroup pwglf_forward_scripts
195 FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
197 TF1* tmpf = new TF1("tmpf", "gaus", xlow, xhigh);
198 tmp->Fit(tmpf, "N", "");
199 tmp->SetDirectory(0);
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),
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","");
213 //____________________________________________________________________
215 * Make band of systematic errors
217 * @param tmp Histogram
220 * @ingroup pwglf_forward_scripts
223 MakeSysError(TH1* tmp, TF1* fit)
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);
236 //____________________________________________________________________
238 * Transform a graph into a histogram
244 * @ingroup pwglf_forward_scripts
247 Graph2Hist(const TGraphAsymmErrors* g)
249 Int_t nBins = g->GetN();
250 TArrayF bins(nBins+1);
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;
269 y.fArray[i] = g->GetY()[i];
270 ey.fArray[i] = TMath::Max(g->GetEYlow()[i],g->GetEYhigh()[i]);
273 TString name(g->GetName());
274 TString title(g->GetTitle());
277 h = new TH1D(name.Data(), title.Data(), nBins,
278 bins[0]-dx/2, bins[nBins]+dx/2);
281 h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
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]);
287 h->SetMarkerStyle(g->GetMarkerStyle());
288 h->SetMarkerColor(g->GetMarkerColor());
289 h->SetMarkerSize(g->GetMarkerSize());
295 //____________________________________________________________________
297 * Calculate ratio of histogram to function
301 * @param title (Optional) title
303 * @return Ratio in a histogram
305 * @ingroup pwglf_forward_scripts
308 Ratio(TH1* h, TF1* f, const char* title)
310 TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_%s",
313 ret->SetDirectory(0);
314 if (title) ret->SetTitle(title);
315 else ret->SetTitle(Form("%s (data) / %s",
316 h->GetTitle(),f->GetTitle()));
318 for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
319 Double_t cc = h->GetBinContent(i);
321 ret->SetBinContent(i,0);
322 ret->SetBinError(i,0);
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);
336 //____________________________________________________________________
340 * @param type Trigger type (1: INEL, 4: NSD)
341 * @param p On return, positive part or null
342 * @param n On return, negative part or null
343 * @param xlow On return, lower X bound
344 * @param xhigh On return, upper X bound
346 * @return Merged histogram or null
348 * @ingroup pwglf_forward_scripts
351 GetUA5Data(UShort_t type, TH1*& p, TH1*& n,
352 Double_t& xlow, Double_t& xhigh)
354 gROOT->SetMacroPath(Form(".:$ALICE_ROOT.trunk/PWGLF/FORWARD/analysis2/:%s",
355 gROOT->GetMacroPath()));
356 gROOT->LoadMacro("OtherData.C");
360 TGraphAsymmErrors* gp = GetSingle(UA5, 1, 900, type, 0, 0);
361 TGraphAsymmErrors* gn = GetSingle(UA5+10, 1, 900, type, 0, 0);
362 if (!gp || !gn) return 0;
367 Int_t nn = p->GetNbinsX();
368 xlow = n->GetXaxis()->GetXmin();
369 xhigh = p->GetXaxis()->GetXmax();
370 TH1D* ret = new TH1D("ua5", "UA5", 2*nn, xlow, xhigh);
371 ret->SetDirectory(0);
372 ret->SetMarkerColor(p->GetMarkerColor());
373 ret->SetMarkerStyle(p->GetMarkerStyle());
375 for (Int_t i = 1; i <= nn; i++) {
376 ret->SetBinContent(nn+i, p->GetBinContent(i));
377 ret->SetBinContent( i, n->GetBinContent(i));
378 ret->SetBinError(nn+i, p->GetBinError(i));
379 ret->SetBinError( i, n->GetBinError(i));
385 //____________________________________________________________________
387 * Draw ratios to UA5 data
389 * @param fname Input file name
390 * @param rebin Re-binning factor
392 * @ingroup pwglf_forward_scripts
395 DrawUA5Ratios(const char* fname="forward_dndeta.root", UShort_t rebin=5)
397 TFile* forward_dndeta = TFile::Open(fname, "READ");
398 if (!forward_dndeta) {
399 Error("DrawUA5Ratios", "%s not found", fname);
404 TH1D* forward = GetHist(forward_dndeta, "Forward", rebin);
405 TH1D* central = GetHist(forward_dndeta, "Central", rebin);
407 TObject* sys = GetObject(forward_dndeta->Get("ForwardResults"), "sys");
408 TObject* sNN = GetObject(forward_dndeta->Get("ForwardResults"), "sNN");
409 TObject* trg = GetObject(forward_dndeta->Get("ForwardResults"), "trigger");
410 if (!sys || !sNN || !trg) return;
411 if (sys->GetUniqueID() != 1) {
412 Error("DrawUA5Ratios", "Comparison only valid for pp, not %s",
416 if (sNN->GetUniqueID() != 900) {
417 Error("DrawUA5Ratios", "Comparison only valid for 900GeV, not %dGeV",
421 if (trg->GetUniqueID() != 1 &&
422 trg->GetUniqueID() != 4) {
423 Error("DrawUA5Ratios",
424 "Comparison only valid for INEL or NSD, not %s (%d)",
425 trg->GetTitle(), trg->GetUniqueID());
430 Double_t alilow, alihigh;
431 TH1D* ali = Merge(central, forward, alilow, alihigh);
432 TF1* alifit = FitMerged(ali, alilow, alihigh);
433 ali->SetTitle("Forward+Central");
434 alifit->SetLineColor(kRed+1);
435 alifit->SetLineStyle(2);
436 alifit->SetName("alifit");
437 alifit->SetTitle("Fit to Forward+Central");
439 Double_t ua5low, ua5high;
441 TH1D* ua5 = GetUA5Data(trg->GetUniqueID(), ua5p, ua5n, ua5low, ua5high);
442 TF1* ua5fit = FitMerged(ua5, ua5low, ua5high);
443 ua5fit->SetLineColor(kBlue+1);
444 ua5fit->SetLineStyle(3);
445 ua5fit->SetName("ua5fit");
446 ua5fit->SetTitle("Fit to UA5");
448 gStyle->SetOptTitle(0);
449 TCanvas* c = new TCanvas("c", "C", 900, 900);
457 TPad* p1 = new TPad("p1", "p1", 0, yd, 1, 1, 0, 0, 0);
458 p1->SetBorderSize(0);
459 p1->SetBorderMode(0);
462 p1->SetRightMargin(0.02);
463 p1->SetTopMargin(0.02);
464 p1->SetBottomMargin(0.00);
469 THStack* all = new THStack("all", "All");
476 all->Draw("nostack");
477 all->SetMinimum(-.07);
478 all->GetXaxis()->SetTitleFont(132);
479 all->GetYaxis()->SetTitleFont(132);
480 all->GetXaxis()->SetLabelFont(132);
481 all->GetYaxis()->SetLabelFont(132);
482 all->GetYaxis()->SetDecimals();
484 all->Draw("nostack");
485 // ua5p->Draw("same p");
486 // ua5m->Draw("same p");
487 alifit->Draw("same");
488 ua5fit->Draw("same");
490 TLegend* l = new TLegend(.2, .1, .8, .5,
491 Form("pp @ #sqrt{s}=900GeV, %s",trg->GetTitle()));
492 l->AddEntry(ua5, Form("U: %s", ua5->GetTitle()), "lp");
493 l->AddEntry(forward, "F: Forward", "lp");
494 l->AddEntry(central, "C: Central", "lp");
496 Form("f: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
497 "%4.2fe^{(x/%4.2f)^{2}}#right]",
499 alifit->GetParameter(0),
500 alifit->GetParameter(2),
501 alifit->GetParameter(1),
502 alifit->GetParameter(3)), "l");
504 Form("u: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
505 "%4.2fe^{(x/%4.2f)^{2}}#right]",
507 ua5fit->GetParameter(0),
508 ua5fit->GetParameter(2),
509 ua5fit->GetParameter(1),
510 ua5fit->GetParameter(3)), "l");
518 TPad* p2 = new TPad("p2", "p2", 0, 0, 1, yd, 0, 0, 0);
519 p2->SetBorderSize(0);
520 p2->SetBorderMode(0);
523 p2->SetRightMargin(0.02);
524 p2->SetTopMargin(0.00);
525 p2->SetBottomMargin(0.15);
530 THStack* ratios = new THStack("Ratios", "Ratios");
531 TH1* ua5ali = Ratio(ua5, alifit, 0);
532 TH1* aliua5 = Ratio(ali, ua5fit, 0);
535 ratios->Draw("nostack");
536 ratios->SetMinimum(0.4);
537 ratios->GetYaxis()->SetTitleSize(2*ratios->GetYaxis()->GetTitleSize());
538 ratios->GetYaxis()->SetLabelSize(2*ratios->GetYaxis()->GetLabelSize());
539 ratios->GetYaxis()->SetNdivisions(508);
540 ratios->GetXaxis()->SetTitleSize(2*ratios->GetXaxis()->GetTitleSize());
541 ratios->GetXaxis()->SetLabelSize(2*ratios->GetXaxis()->GetLabelSize());
542 ratios->GetXaxis()->SetNdivisions(510);
543 ratios->GetXaxis()->SetTitle("#eta");
544 ratios->GetXaxis()->SetTitleFont(132);
545 ratios->GetYaxis()->SetTitleFont(132);
546 ratios->GetXaxis()->SetLabelFont(132);
547 ratios->GetYaxis()->SetLabelFont(132);
548 ratios->GetYaxis()->SetDecimals();
551 TGraphErrors* sysErr = new TGraphErrors(2);
552 sysErr->SetPoint(0, all->GetHistogram()->GetXaxis()->GetXmin(),1);
553 sysErr->SetPoint(1, all->GetHistogram()->GetXaxis()->GetXmax(),1);
554 sysErr->SetPointError(0,0,.07);
555 sysErr->SetPointError(1,0,.07);
556 sysErr->SetTitle("Systematic error on Forward+Central");
557 sysErr->SetFillColor(kYellow+1);
558 sysErr->SetFillStyle(3001);
559 sysErr->SetHistogram(ratios->GetHistogram());
560 sysErr->DrawClone("ael3");
561 ratios->DrawClone("nostack same");
563 TF1* fitfit = new TF1("fitfit", myRatio, alilow, alihigh, 8);
564 fitfit->SetParameters(ua5fit->GetParameter(0),
565 ua5fit->GetParameter(1),
566 ua5fit->GetParameter(2),
567 ua5fit->GetParameter(3),
568 alifit->GetParameter(0),
569 alifit->GetParameter(1),
570 alifit->GetParameter(2),
571 alifit->GetParameter(3));
572 fitfit->Draw("same");
574 TLegend* ll = new TLegend(.3,.15, .7, .6);
575 ll->AddEntry(sysErr, sysErr->GetTitle(), "f");
576 ll->AddEntry(ua5ali, ua5ali->GetTitle(), "lp");
577 ll->AddEntry(aliua5, aliua5->GetTitle(), "lp");
578 ll->AddEntry(fitfit, "UA5 (fit) / Forward+Central (fit)", "lp");
579 ll->SetTextFont(132);
580 ll->SetBorderSize(0);
586 c->SaveAs(Form("ua5_ratios_%s_r%02d.png", trg->GetTitle(), rebin));
587 gROOT->GetInterpreter()->UnloadFile("OtherData.C");
591 //____________________________________________________________________