PWG2/FORWARD/analysis2/scripts/DrawUA5Ratios.C

   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
   9  *
  10  * @ingroup pwg2_forward_scripts
  11  */
  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 //____________________________________________________________________
  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
  38  *
  39  * @ingroup pwg2_forward_scripts
  40  */
  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 //____________________________________________________________________
  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
  75  *
  76  * @ingroup pwg2_forward_scripts
  77  */
  78 TH1D*
  79 GetHist(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
  99  *
 100  * @ingroup pwg2_forward_scripts
 101  */
 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 //____________________________________________________________________
 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$
 148  *
 149  * @ingroup pwg2_forward_scripts
 150  */
 151 Double_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$
 172  *
 173  * @ingroup pwg2_forward_scripts
 174  */
 175 Double_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
 191  *
 192  * @ingroup pwg2_forward_scripts
 193  */
 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 //____________________________________________________________________
 214 /**
 215  * Make band of systematic errors
 216  *
 217  * @param tmp Histogram
 218  * @param fit Fit
 219  *
 220  * @ingroup pwg2_forward_scripts
 221  */
 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 //____________________________________________________________________
 237 /**
 238  * Transform a graph into a histogram
 239  *
 240  * @param g
 241  *
 242  * @return
 243  *
 244  * @ingroup pwg2_forward_scripts
 245  */
 246 TH1*
 247 Graph2Hist(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
 304  *
 305  * @ingroup pwg2_forward_scripts
 306  */
 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 //____________________________________________________________________
 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
 347  *
 348  * @ingroup pwg2_forward_scripts
 349  */
 350 TH1D*
 351 GetUA5Data(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 /**
 387  * Draw ratios to UA5 data
 388  *
 389  *
 390  * @ingroup pwg2_forward_scripts
 391  */
 392 void
 393 DrawUA5Ratios(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 //