]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/scripts/DrawUA5Ratios.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'TPCdev' of https://git.cern.ch/reps/AliRoot into TPCdev
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / DrawUA5Ratios.C
CommitLineData
797161e8 1//____________________________________________________________________
2/**
3 * Get an object with specified name from TCollection @a l
4 *
5 * @param l Collection
6 * @param name Name of object to retrieve
7 *
8 * @return Object, or null
56199f2b 9 *
bd6f5206 10 * @ingroup pwglf_forward_scripts
797161e8 11 */
12TObject*
13GetObject(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//____________________________________________________________________
29/**
30 * Get histogram from @a directory/which/sub/hname
31 *
32 * @param dir Directory
33 * @param which Name of parent list
34 * @param sub Name of sub-list
35 * @param hname Name of histogram
36 *
37 * @return Pointer to histogram or null
56199f2b 38 *
bd6f5206 39 * @ingroup pwglf_forward_scripts
797161e8 40 */
41TH1D*
42GetHist(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//____________________________________________________________________
65/**
66 * Get histogram from
67 * <i>dir/which</i>Results<i>/sub/</i>dndeta<i>which</i>[_rebin<i>rebin</i>]
68 *
69 * @param dir Directory
70 * @param which Name
71 * @param rebin Optional rebinning
72 * @param sub Sub-list name
73 *
74 * @return Histogram pointer or null
56199f2b 75 *
bd6f5206 76 * @ingroup pwglf_forward_scripts
797161e8 77 */
78TH1D*
79GetHist(TDirectory* dir,
80 const char* which,
81 UShort_t rebin,
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//____________________________________________________________________
90/**
91 * Merge two histograms into one
92 *
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
97 *
98 * @return Newly allocated histogram or null
56199f2b 99 *
bd6f5206 100 * @ingroup pwglf_forward_scripts
797161e8 101 */
102TH1*
103Merge(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//____________________________________________________________________
136/**
137 * Function to calculate
138 * @f[
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)
142 * @f]
143 *
144 * @param xp Pointer to x array
145 * @param pp Pointer to parameter array
146 *
147 * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
56199f2b 148 *
bd6f5206 149 * @ingroup pwglf_forward_scripts
797161e8 150 */
151Double_t myFunc(Double_t* xp, Double_t* pp)
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//____________________________________________________________________
161/**
162 * Calculate
163 * @f[
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)}
166 * @f]
167 *
168 * @param xp Pointer to X array
169 * @param pp Pointer to parameter array (8 entries)
170 *
171 * @return @f$r(x)@f$
56199f2b 172 *
bd6f5206 173 * @ingroup pwglf_forward_scripts
797161e8 174 */
175Double_t myRatio(Double_t* xp, Double_t* pp)
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//____________________________________________________________________
183/**
184 * Fit @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
185 *
186 * @param tmp Histogram
187 * @param xlow Lower x bound
188 * @param xhigh Upper x bound
189 *
190 * @return Fitted function
56199f2b 191 *
bd6f5206 192 * @ingroup pwglf_forward_scripts
797161e8 193 */
194TF1*
195FitMerged(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//____________________________________________________________________
214/**
215 * Make band of systematic errors
216 *
217 * @param tmp Histogram
218 * @param fit Fit
56199f2b 219 *
bd6f5206 220 * @ingroup pwglf_forward_scripts
797161e8 221 */
222void
223MakeSysError(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//____________________________________________________________________
237/**
238 * Transform a graph into a histogram
239 *
240 * @param g
241 *
242 * @return
56199f2b 243 *
bd6f5206 244 * @ingroup pwglf_forward_scripts
56199f2b 245 */
797161e8 246TH1*
247Graph2Hist(const TGraphAsymmErrors* g)
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//____________________________________________________________________
296/**
297 * Calculate ratio of histogram to function
298 *
299 * @param h Histogram
300 * @param f Function
301 * @param title (Optional) title
302 *
303 * @return Ratio in a histogram
56199f2b 304 *
bd6f5206 305 * @ingroup pwglf_forward_scripts
797161e8 306 */
307TH1*
308Ratio(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//____________________________________________________________________
337/**
338 * Get the UA5 data
339 *
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
345 *
346 * @return Merged histogram or null
56199f2b 347 *
bd6f5206 348 * @ingroup pwglf_forward_scripts
797161e8 349 */
350TH1D*
351GetUA5Data(UShort_t type, TH1*& p, TH1*& n,
352 Double_t& xlow, Double_t& xhigh)
353{
bd6f5206 354 gROOT->SetMacroPath(Form(".:$ALICE_ROOT.trunk/PWGLF/FORWARD/analysis2/:%s",
797161e8 355 gROOT->GetMacroPath()));
356 gROOT->LoadMacro("OtherData.C");
357
358 p = 0;
359 n = 0;
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;
363
364 p = Graph2Hist(gp);
365 n = Graph2Hist(gn);
366
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());
374
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));
380 }
381 return ret;
382}
383
384
385//____________________________________________________________________
386/**
56199f2b 387 * Draw ratios to UA5 data
797161e8 388 *
95c09e1c 389 * @param fname Input file name
390 * @param rebin Re-binning factor
56199f2b 391 *
bd6f5206 392 * @ingroup pwglf_forward_scripts
797161e8 393 */
394void
395DrawUA5Ratios(const char* fname="forward_dndeta.root", UShort_t rebin=5)
396{
397 TFile* forward_dndeta = TFile::Open(fname, "READ");
398 if (!forward_dndeta) {
399 Error("DrawUA5Ratios", "%s not found", fname);
400 return;
401 }
402
403 UShort_t type = 1;
404 TH1D* forward = GetHist(forward_dndeta, "Forward", rebin);
405 TH1D* central = GetHist(forward_dndeta, "Central", rebin);
406
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",
413 sys->GetTitle());
414 return;
415 }
416 if (sNN->GetUniqueID() != 900) {
417 Error("DrawUA5Ratios", "Comparison only valid for 900GeV, not %dGeV",
418 sNN->GetUniqueID());
419 return;
420 }
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());
426 return;
427 }
428
429
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");
438
439 Double_t ua5low, ua5high;
440 TH1* ua5p, *ua5n;
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");
447
448 gStyle->SetOptTitle(0);
449 TCanvas* c = new TCanvas("c", "C", 900, 900);
450 c->SetFillColor(0);
451 c->SetFillStyle(0);
452 c->SetBorderMode(0);
453 c->SetBorderSize(0);
454
455 Double_t yd = .3;
456
457 TPad* p1 = new TPad("p1", "p1", 0, yd, 1, 1, 0, 0, 0);
458 p1->SetBorderSize(0);
459 p1->SetBorderMode(0);
460 p1->SetFillColor(0);
461 p1->SetFillStyle(0);
462 p1->SetRightMargin(0.02);
463 p1->SetTopMargin(0.02);
464 p1->SetBottomMargin(0.00);
465 p1->SetGridx();
466 p1->Draw();
467 p1->cd();
468
469 THStack* all = new THStack("all", "All");
470 all->Add(ua5p);
471 all->Add(ua5n);
472 // all->Add(ua5);
473 all->Add(forward);
474 all->Add(central);
475 // all->Add(merged);
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();
483 p1->Clear();
484 all->Draw("nostack");
485 // ua5p->Draw("same p");
486 // ua5m->Draw("same p");
487 alifit->Draw("same");
488 ua5fit->Draw("same");
489
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");
495 l->AddEntry(alifit,
496 Form("f: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
497 "%4.2fe^{(x/%4.2f)^{2}}#right]",
498 alifit->GetTitle(),
499 alifit->GetParameter(0),
500 alifit->GetParameter(2),
501 alifit->GetParameter(1),
502 alifit->GetParameter(3)), "l");
503 l->AddEntry(ua5fit,
504 Form("u: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
505 "%4.2fe^{(x/%4.2f)^{2}}#right]",
506 ua5fit->GetTitle(),
507 ua5fit->GetParameter(0),
508 ua5fit->GetParameter(2),
509 ua5fit->GetParameter(1),
510 ua5fit->GetParameter(3)), "l");
511 l->SetTextFont(132);
512 l->SetBorderSize(0);
513 l->SetFillColor(0);
514 l->SetFillStyle(0);
515 l->Draw();
516
517 c->cd();
518 TPad* p2 = new TPad("p2", "p2", 0, 0, 1, yd, 0, 0, 0);
519 p2->SetBorderSize(0);
520 p2->SetBorderMode(0);
521 p2->SetFillColor(0);
522 p2->SetFillStyle(0);
523 p2->SetRightMargin(0.02);
524 p2->SetTopMargin(0.00);
525 p2->SetBottomMargin(0.15);
526 p2->SetGridx();
527 p2->Draw();
528 p2->cd();
529
530 THStack* ratios = new THStack("Ratios", "Ratios");
531 TH1* ua5ali = Ratio(ua5, alifit, 0);
532 TH1* aliua5 = Ratio(ali, ua5fit, 0);
533 ratios->Add(ua5ali);
534 ratios->Add(aliua5);
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();
549 p2->Clear();
550
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");
562
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");
573
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);
581 ll->SetFillColor(0);
582 ll->SetFillStyle(0);
583 ll->Draw();
584
585
586 c->SaveAs(Form("ua5_ratios_%s_r%02d.png", trg->GetTitle(), rebin));
587 gROOT->GetInterpreter()->UnloadFile("OtherData.C");
588
589}
590
591//____________________________________________________________________
592//
593// EOF
594//