]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/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
Doxygen documentation fixes
[u/mrichter/AliRoot.git] / PWG2 / 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 *
10 * @ingroup pwg2_forward_analysis_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 *
39 * @ingroup pwg2_forward_analysis_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 *
76 * @ingroup pwg2_forward_analysis_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 *
100 * @ingroup pwg2_forward_analysis_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 *
149 * @ingroup pwg2_forward_analysis_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 *
173 * @ingroup pwg2_forward_analysis_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 *
192 * @ingroup pwg2_forward_analysis_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 *
220 * @ingroup pwg2_forward_analysis_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 *
244 * @ingroup pwg2_forward_analysis_scripts
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 *
305 * @ingroup pwg2_forward_analysis_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 *
348 * @ingroup pwg2_forward_analysis_scripts
797161e8 349 */
350TH1D*
351GetUA5Data(UShort_t type, TH1*& p, TH1*& n,
352 Double_t& xlow, Double_t& xhigh)
353{
354 gROOT->SetMacroPath(Form(".:$ALICE_ROOT.trunk/PWG2/FORWARD/analysis2/:%s",
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 *
56199f2b 389 *
390 * @ingroup pwg2_forward_analysis_scripts
797161e8 391 */
392void
393DrawUA5Ratios(const char* fname="forward_dndeta.root", UShort_t rebin=5)
394{
395 TFile* forward_dndeta = TFile::Open(fname, "READ");
396 if (!forward_dndeta) {
397 Error("DrawUA5Ratios", "%s not found", fname);
398 return;
399 }
400
401 UShort_t type = 1;
402 TH1D* forward = GetHist(forward_dndeta, "Forward", rebin);
403 TH1D* central = GetHist(forward_dndeta, "Central", rebin);
404
405 TObject* sys = GetObject(forward_dndeta->Get("ForwardResults"), "sys");
406 TObject* sNN = GetObject(forward_dndeta->Get("ForwardResults"), "sNN");
407 TObject* trg = GetObject(forward_dndeta->Get("ForwardResults"), "trigger");
408 if (!sys || !sNN || !trg) return;
409 if (sys->GetUniqueID() != 1) {
410 Error("DrawUA5Ratios", "Comparison only valid for pp, not %s",
411 sys->GetTitle());
412 return;
413 }
414 if (sNN->GetUniqueID() != 900) {
415 Error("DrawUA5Ratios", "Comparison only valid for 900GeV, not %dGeV",
416 sNN->GetUniqueID());
417 return;
418 }
419 if (trg->GetUniqueID() != 1 &&
420 trg->GetUniqueID() != 4) {
421 Error("DrawUA5Ratios",
422 "Comparison only valid for INEL or NSD, not %s (%d)",
423 trg->GetTitle(), trg->GetUniqueID());
424 return;
425 }
426
427
428 Double_t alilow, alihigh;
429 TH1D* ali = Merge(central, forward, alilow, alihigh);
430 TF1* alifit = FitMerged(ali, alilow, alihigh);
431 ali->SetTitle("Forward+Central");
432 alifit->SetLineColor(kRed+1);
433 alifit->SetLineStyle(2);
434 alifit->SetName("alifit");
435 alifit->SetTitle("Fit to Forward+Central");
436
437 Double_t ua5low, ua5high;
438 TH1* ua5p, *ua5n;
439 TH1D* ua5 = GetUA5Data(trg->GetUniqueID(), ua5p, ua5n, ua5low, ua5high);
440 TF1* ua5fit = FitMerged(ua5, ua5low, ua5high);
441 ua5fit->SetLineColor(kBlue+1);
442 ua5fit->SetLineStyle(3);
443 ua5fit->SetName("ua5fit");
444 ua5fit->SetTitle("Fit to UA5");
445
446 gStyle->SetOptTitle(0);
447 TCanvas* c = new TCanvas("c", "C", 900, 900);
448 c->SetFillColor(0);
449 c->SetFillStyle(0);
450 c->SetBorderMode(0);
451 c->SetBorderSize(0);
452
453 Double_t yd = .3;
454
455 TPad* p1 = new TPad("p1", "p1", 0, yd, 1, 1, 0, 0, 0);
456 p1->SetBorderSize(0);
457 p1->SetBorderMode(0);
458 p1->SetFillColor(0);
459 p1->SetFillStyle(0);
460 p1->SetRightMargin(0.02);
461 p1->SetTopMargin(0.02);
462 p1->SetBottomMargin(0.00);
463 p1->SetGridx();
464 p1->Draw();
465 p1->cd();
466
467 THStack* all = new THStack("all", "All");
468 all->Add(ua5p);
469 all->Add(ua5n);
470 // all->Add(ua5);
471 all->Add(forward);
472 all->Add(central);
473 // all->Add(merged);
474 all->Draw("nostack");
475 all->SetMinimum(-.07);
476 all->GetXaxis()->SetTitleFont(132);
477 all->GetYaxis()->SetTitleFont(132);
478 all->GetXaxis()->SetLabelFont(132);
479 all->GetYaxis()->SetLabelFont(132);
480 all->GetYaxis()->SetDecimals();
481 p1->Clear();
482 all->Draw("nostack");
483 // ua5p->Draw("same p");
484 // ua5m->Draw("same p");
485 alifit->Draw("same");
486 ua5fit->Draw("same");
487
488 TLegend* l = new TLegend(.2, .1, .8, .5,
489 Form("pp @ #sqrt{s}=900GeV, %s",trg->GetTitle()));
490 l->AddEntry(ua5, Form("U: %s", ua5->GetTitle()), "lp");
491 l->AddEntry(forward, "F: Forward", "lp");
492 l->AddEntry(central, "C: Central", "lp");
493 l->AddEntry(alifit,
494 Form("f: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
495 "%4.2fe^{(x/%4.2f)^{2}}#right]",
496 alifit->GetTitle(),
497 alifit->GetParameter(0),
498 alifit->GetParameter(2),
499 alifit->GetParameter(1),
500 alifit->GetParameter(3)), "l");
501 l->AddEntry(ua5fit,
502 Form("u: %s: %4.2f#left[e^{(x/%4.2f)^{2}}-"
503 "%4.2fe^{(x/%4.2f)^{2}}#right]",
504 ua5fit->GetTitle(),
505 ua5fit->GetParameter(0),
506 ua5fit->GetParameter(2),
507 ua5fit->GetParameter(1),
508 ua5fit->GetParameter(3)), "l");
509 l->SetTextFont(132);
510 l->SetBorderSize(0);
511 l->SetFillColor(0);
512 l->SetFillStyle(0);
513 l->Draw();
514
515 c->cd();
516 TPad* p2 = new TPad("p2", "p2", 0, 0, 1, yd, 0, 0, 0);
517 p2->SetBorderSize(0);
518 p2->SetBorderMode(0);
519 p2->SetFillColor(0);
520 p2->SetFillStyle(0);
521 p2->SetRightMargin(0.02);
522 p2->SetTopMargin(0.00);
523 p2->SetBottomMargin(0.15);
524 p2->SetGridx();
525 p2->Draw();
526 p2->cd();
527
528 THStack* ratios = new THStack("Ratios", "Ratios");
529 TH1* ua5ali = Ratio(ua5, alifit, 0);
530 TH1* aliua5 = Ratio(ali, ua5fit, 0);
531 ratios->Add(ua5ali);
532 ratios->Add(aliua5);
533 ratios->Draw("nostack");
534 ratios->SetMinimum(0.4);
535 ratios->GetYaxis()->SetTitleSize(2*ratios->GetYaxis()->GetTitleSize());
536 ratios->GetYaxis()->SetLabelSize(2*ratios->GetYaxis()->GetLabelSize());
537 ratios->GetYaxis()->SetNdivisions(508);
538 ratios->GetXaxis()->SetTitleSize(2*ratios->GetXaxis()->GetTitleSize());
539 ratios->GetXaxis()->SetLabelSize(2*ratios->GetXaxis()->GetLabelSize());
540 ratios->GetXaxis()->SetNdivisions(510);
541 ratios->GetXaxis()->SetTitle("#eta");
542 ratios->GetXaxis()->SetTitleFont(132);
543 ratios->GetYaxis()->SetTitleFont(132);
544 ratios->GetXaxis()->SetLabelFont(132);
545 ratios->GetYaxis()->SetLabelFont(132);
546 ratios->GetYaxis()->SetDecimals();
547 p2->Clear();
548
549 TGraphErrors* sysErr = new TGraphErrors(2);
550 sysErr->SetPoint(0, all->GetHistogram()->GetXaxis()->GetXmin(),1);
551 sysErr->SetPoint(1, all->GetHistogram()->GetXaxis()->GetXmax(),1);
552 sysErr->SetPointError(0,0,.07);
553 sysErr->SetPointError(1,0,.07);
554 sysErr->SetTitle("Systematic error on Forward+Central");
555 sysErr->SetFillColor(kYellow+1);
556 sysErr->SetFillStyle(3001);
557 sysErr->SetHistogram(ratios->GetHistogram());
558 sysErr->DrawClone("ael3");
559 ratios->DrawClone("nostack same");
560
561 TF1* fitfit = new TF1("fitfit", myRatio, alilow, alihigh, 8);
562 fitfit->SetParameters(ua5fit->GetParameter(0),
563 ua5fit->GetParameter(1),
564 ua5fit->GetParameter(2),
565 ua5fit->GetParameter(3),
566 alifit->GetParameter(0),
567 alifit->GetParameter(1),
568 alifit->GetParameter(2),
569 alifit->GetParameter(3));
570 fitfit->Draw("same");
571
572 TLegend* ll = new TLegend(.3,.15, .7, .6);
573 ll->AddEntry(sysErr, sysErr->GetTitle(), "f");
574 ll->AddEntry(ua5ali, ua5ali->GetTitle(), "lp");
575 ll->AddEntry(aliua5, aliua5->GetTitle(), "lp");
576 ll->AddEntry(fitfit, "UA5 (fit) / Forward+Central (fit)", "lp");
577 ll->SetTextFont(132);
578 ll->SetBorderSize(0);
579 ll->SetFillColor(0);
580 ll->SetFillStyle(0);
581 ll->Draw();
582
583
584 c->SaveAs(Form("ua5_ratios_%s_r%02d.png", trg->GetTitle(), rebin));
585 gROOT->GetInterpreter()->UnloadFile("OtherData.C");
586
587}
588
589//____________________________________________________________________
590//
591// EOF
592//