PWG2/FORWARD/analysis2/scripts/SimpledNdeta.C

   1 #ifndef __CINT__
   2 # include <TStyle.h>
   3 # include <TFile.h>
   4 # include <TList.h>
   5 # include <TH1.h>
   6 # include <THStack.h>
   7 # include <TLegend.h>
   8 # include <TLegendEntry.h>
   9 # include <TLatex.h>
  10 # include <TLine.h>
  11 # include <TString.h>
  12 # include <TCanvas.h>
  13 # include <TError.h>
  14 # include <TColor.h>
  15 #else
  16 class THStack;
  17 class TAxis;
  18 #endif
  19
  20 /**
  21  * A simple script to draw results from MakedNdeta.C (or similar)
  22  *
  23  */
  24 /**
  25  * Get a stack from the passed list
  26  *
  27  * @param list   List to get the stack from
  28  * @param name   Name of stack
  29  * @param rebin  Optional rebinning - must exists in list
  30  *
  31  * @return Stack or null
  32  */
  33 THStack*
  34 GetStack(const TList* list, const char* name, Int_t rebin)
  35 {
  36   if (!list) {
  37     Warning("GetStack", "List is null");
  38     return 0;
  39   }
  40
  41   TString n(name);
  42   if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
  43   TObject* o = list->FindObject(n);
  44   if (!o) {
  45     Warning("GetStack", "No %s object found in %s", n.Data(), list->GetName());
  46     return 0;
  47   }
  48   return static_cast<THStack*>(o);
  49 }
  50
  51 TH1*
  52 GetHist(const TList* list, const char* name, Int_t rebin)
  53 {
  54   if (!list) {
  55     Warning("GetStack", "List is null");
  56     return 0;
  57   }
  58   TList* all = static_cast<TList*>(list->FindObject("all"));
  59   if (!all) {
  60     Warning("GetHist", "List all not found in %s", list->GetName());
  61     return 0;
  62   }
  63
  64   TString n(name);
  65   if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
  66   TObject* o = all->FindObject(n);
  67   if (!o) {
  68     Warning("GetStack", "No %s object found in %s", n.Data(), all->GetName());
  69     return 0;
  70   }
  71   return static_cast<TH1*>(o);
  72 }
  73
  74 /**
  75  * Add histograms from one stack to another
  76  *
  77  * @param p      Parent stack to add to
  78  * @param list   List to look for the stack in
  79  * @param name   Name of stack to add
  80  * @param rebin  Optional rebinning - must exists in list
  81  *
  82  * @return Added stack or null
  83  */
  84 THStack*
  85 AddStack(THStack* p, const TList* list, const char* name, Int_t rebin)
  86 {
  87   THStack* s = GetStack(list, name, rebin);
  88   if (!s) return 0;
  89
  90   TIter next(s->GetHists());
  91   TH1*  hist = 0;
  92   while ((hist = static_cast<TH1*>(next())))
  93     p->Add(hist);
  94   return s;
  95 }
  96
  97 /**
  98  * Build up a centrality legend
  99  *
 100  * @param c Centrality axis
 101  */
 102 void
 103 BuildCentLegend(const TAxis* c)
 104 {
 105   if (!c) return;
 106   TLegend* l = new TLegend(1.2*gPad->GetLeftMargin(),
 107                            .55, .35, 1-gPad->GetTopMargin());
 108   l->SetFillColor(0);
 109   l->SetFillStyle(0);
 110   l->SetBorderSize(0);
 111   l->SetTextFont(22);
 112   l->SetHeader("Centralities");
 113   l->SetTextFont(132);
 114
 115   Int_t    nCol     = gStyle->GetNumberOfColors();
 116   for (Int_t i = 0; i < c->GetNbins(); i++) {
 117     UShort_t low    = c->GetBinLowEdge(i+1);
 118     UShort_t high   = c->GetBinUpEdge(i+1);
 119     Float_t  fc     = (low+double(high-low)/2) / 100;
 120     Int_t    icol   = TMath::Min(nCol-1,int(fc * nCol + .5));
 121     Int_t    col    = gStyle->GetColorPalette(icol);
 122     TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
 123                                   Form("%3d%% - %3d%%", low, high), "p");
 124     e->SetMarkerColor(col);
 125
 126   }
 127   l->Draw();
 128 }
 129
 130 /**
 131  * Build a legend.  Histograms a filtered for the same title
 132  *
 133  * @param stack Stack of histograms
 134  * @param c     Centrality axis.  If present, markers are black
 135  */
 136 void
 137 BuildLegend(const THStack* stack, const TAxis* c)
 138 {
 139   Double_t x1 = .75, x2 = 1-gPad->GetRightMargin();
 140   Double_t y1 = .8,  y2 = 1-gPad->GetTopMargin();
 141   if (!c) {
 142     // PP
 143     x1 = .4;
 144     y1 = .4;
 145     x2 = .8;
 146     y2 = .6;
 147   }
 148   TLegend* l = new TLegend(x1, y1, x2, y2, "");
 149   l->SetFillColor(0);
 150   l->SetFillStyle(0);
 151   l->SetBorderSize(0);
 152   l->SetTextFont(132);
 153
 154   TObjArray seen;
 155   seen.SetOwner();
 156   TIter next(stack->GetHists());
 157   TH1* hist = 0;
 158   while ((hist = static_cast<TH1*>(next()))) {
 159     TString n(hist->GetTitle());
 160     if (n.Contains("mirrored")) continue;
 161     if (seen.FindObject(n.Data())) continue;
 162     TObjString* ns = new TObjString(n);
 163     ns->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
 164                     ((hist->GetMarkerColor() & 0xFFFF) <<  0));
 165     seen.Add(ns);
 166   }
 167   seen.Sort();
 168   TIter nextu(&seen);
 169   TObject* s = 0;
 170   Int_t i = 0;
 171   while ((s = nextu())) {
 172     TLegendEntry* dd = l->AddEntry(Form("data%2d", i++),
 173                                    s->GetName(), "p");
 174     Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
 175     Int_t color = (s->GetUniqueID() >>  0) & 0xFFFF;
 176     dd->SetLineColor(kBlack);
 177     if (c) dd->SetMarkerColor(kBlack);
 178     else   dd->SetMarkerColor(color);
 179     dd->SetMarkerStyle(style);
 180     Int_t st = dd->GetMarkerStyle();
 181     if (st == 27 || st == 28 || st == 29 || st == 30 || st == 33 || st == 34)
 182       dd->SetMarkerSize(1.4*dd->GetMarkerSize());
 183   }
 184   // TLegendEntry* sep = l->AddEntry("s", "_", "h");
 185   // sep->SetTextSize(0.01);
 186   TLine* sep = new TLine(0,0,1,1);
 187   sep->SetLineWidth(1);
 188   sep->DrawLineNDC(x1+.02, y1-.005, x2-.03, y1-.005);
 189   // sep->Draw();
 190
 191   TLegend* l2 = new TLegend(x1, y1-.005, x2, y1-.15, "");
 192   l2->SetFillColor(0);
 193   l2->SetFillStyle(0);
 194   l2->SetBorderSize(0);
 195   l2->SetTextFont(132);
 196
 197   TLegendEntry* d1 = l2->AddEntry("d1", "Data", "p");
 198   d1->SetLineColor(kBlack);
 199   d1->SetMarkerColor(kBlack);
 200   d1->SetMarkerStyle(20);
 201   TLegendEntry* d2 = l2->AddEntry("d2", "Mirrored data", "p");
 202   d2->SetLineColor(kBlack);
 203   d2->SetMarkerColor(kBlack);
 204   d2->SetMarkerStyle(24);
 205
 206   l->Draw();
 207   l2->Draw();
 208 }
 209
 210 void
 211 AddInformation(TList* forward, bool prelim=true)
 212 {
 213   Double_t x  = .12;
 214   Double_t y  = .95;
 215   Double_t ts = 0.05;
 216   if (prelim) {
 217     TLatex* wip = new TLatex(x, y, "Work in progress");
 218     wip->SetNDC();
 219     wip->SetTextColor(TColor::GetColor(234,26,46));
 220     wip->SetTextAlign(13);
 221     wip->SetTextFont(132);
 222     wip->SetTextSize(ts);
 223     wip->Draw();
 224     y -= ts;
 225   }
 226
 227   TObject* sNN = forward->FindObject("sNN");
 228   TObject* sys = forward->FindObject("sys");
 229   TObject* trg = forward->FindObject("trigger");
 230   TObject* vtx = forward->FindObject("vtxAxis");
 231   TObject* sch = forward->FindObject("scheme");
 232
 233   TString t(sys->GetTitle());
 234   Bool_t isPP = t == "pp";
 235
 236
 237   TString s = Form("%s @ #sqrt{s%s}=",
 238                    sys->GetTitle(),
 239                    (isPP ? "" : "_{NN}"));
 240   Int_t cms = sNN->GetUniqueID();
 241   if (cms > 1000) s.Append(Form("%5.2fTeV", float(cms)/1000));
 242   else            s.Append(Form("%3dGeV", cms));
 243   s.Append(Form(", %s", trg->GetTitle()));
 244
 245   // if (isPP) { x = .3; y = .3; }
 246   if (isPP) {
 247     x  = 1-gPad->GetRightMargin()-.01;
 248     y  = 1-gPad->GetTopMargin()-.01;
 249     ts = .04;
 250   }
 251   TLatex* ltx = new TLatex(x, y, s.Data());
 252   ltx->SetNDC();
 253   ltx->SetTextColor(TColor::GetColor(41,73,156));
 254   ltx->SetTextAlign((isPP ? 33 : 13));
 255   ltx->SetTextFont(132);
 256   ltx->SetTextSize(ts);
 257   ltx->Draw();
 258   y -= ts;
 259   ltx->DrawLatex(x, y, vtx->GetTitle());
 260   y -= ts;
 261   ltx->DrawLatex(x, y, sch->GetTitle());
 262 }
 263
 264 Double_t myFunc(Double_t* xp, Double_t* pp)
 265 {
 266   Double_t x  = xp[0];
 267   Double_t a1 = pp[0];
 268   Double_t a2 = pp[1];
 269   Double_t s1 = pp[2];
 270   Double_t s2 = pp[3];
 271   return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
 272 }
 273
 274 TH1*
 275 MakeSysError(const TH1* cen, const TH1* fwd, Double_t sysErr=0.7)
 276 {
 277   TH1* tmp = static_cast<TH1*>(fwd->Clone("dndetaFitted"));
 278   tmp->SetMarkerStyle(0);
 279   tmp->SetFillColor(kGray);
 280   tmp->SetFillStyle(3001);
 281   tmp->SetDirectory(0);
 282   Double_t xlow  = 100;
 283   Double_t xhigh = 0;
 284   for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
 285     Double_t cc = cen->GetBinContent(i);
 286     Double_t cf = fwd->GetBinContent(i);
 287     Double_t ec = fwd->GetBinError(i);
 288     Double_t ef = fwd->GetBinError(i);
 289     Double_t nc = cf;
 290     Double_t ne = ef;
 291     if (cc < 0.001 && cf < 0.01) continue;
 292     xlow  = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
 293     xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
 294     if (cc > 0.001) {
 295       nc = cc;
 296       ne = ec;
 297       if (cf > 0.001) {
 298         nc  = (cf + cc) / 2;
 299         ne  = TMath::Sqrt(ec*ec + ef*ef);
 300       }
 301     }
 302     tmp->SetBinContent(i, nc);
 303     tmp->SetBinError(i, ne);
 304   }
 305   TF1* tmpf  = new TF1("tmpf", "gaus", xlow, xhigh);
 306   tmp->Fit(tmpf, "NQ", "");
 307
 308   TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
 309   fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
 310   fit->SetParameters(tmpf->GetParameter(0),
 311                      .2,
 312                      tmpf->GetParameter(2),
 313                      tmpf->GetParameter(2)/4);
 314   tmp->Fit(fit,"0WQ","");
 315   for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
 316     Double_t tc = tmp->GetBinContent(i);
 317     if (tc < 0.01) continue;
 318     Double_t x = tmp->GetXaxis()->GetBinCenter(i);
 319     Double_t y = fit->Eval(x);
 320     tmp->SetBinContent(i, y);
 321     tmp->SetBinError(i,sysErr*y);
 322   }
 323   return tmp;
 324 }
 325
 326 /**
 327  * Function to draw the results from forward_dndeta.root file
 328  *
 329  * @param rebin    Rebinnig.  Note, the data must be present in the file
 330  * @param filename File to open and draw stuff from >
 331  */
 332 void
 333 SimpledNdeta(Int_t what=0x5,
 334              Int_t rebin=5, const char* filename="forward_dndeta.root")
 335 {
 336   TFile* file = TFile::Open(filename, "READ");
 337   if (!file) {
 338     Error("SimpledNdeta", "File %s not found", filename);
 339     return;
 340   }
 341
 342   TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
 343   TList* central = static_cast<TList*>(file->Get("CentralResults"));
 344   TList* mctruth = static_cast<TList*>(file->Get("MCTruthResults"));
 345
 346   THStack* all = new THStack("all", "All");
 347   if (what & 0x1) AddStack(all, forward, "dndeta",      rebin);
 348   if (what & 0x2) AddStack(all, forward, "dndetaMC",    rebin);
 349   if (what & 0x1) AddStack(all, central, "dndeta",      rebin);
 350   if (what & 0x2) AddStack(all, central, "dndetaMC",    rebin);
 351   if (what & 0x4) AddStack(all, mctruth, "dndeta",      rebin);
 352
 353   TH1* tmp = 0;
 354   if (what & 0x1) {
 355     Double_t sysErr = 0.07;
 356     TH1* fwd = GetHist(forward, "dndetaForward", rebin);
 357     TH1* cen = GetHist(central, "dndetaCentral", rebin);
 358     tmp = MakeSysError(cen, fwd, sysErr);
 359     // fit = tmpf;
 360   }
 361
 362   TAxis* centAxis = static_cast<TAxis*>(forward->FindObject("centAxis"));
 363   Bool_t isPP = centAxis == 0;
 364
 365   gStyle->SetPalette(1);
 366   gStyle->SetOptTitle(0);
 367   TCanvas* c = new TCanvas("dndeta", "dN/deta results", 900, 600);
 368   c->SetFillColor(0);
 369   c->SetFillStyle(0);
 370   c->SetBorderSize(0);
 371   c->SetBorderMode(0);
 372   c->SetTopMargin(0.03);
 373   c->SetRightMargin(0.03);
 374
 375   all->Draw("nostack");
 376   TAxis* xa = all->GetHistogram()->GetXaxis();
 377   xa->SetTitleFont(132);
 378   xa->SetLabelFont(132);
 379   xa->SetTitle("#eta");
 380   TAxis* ya = all->GetHistogram()->GetYaxis();
 381   ya->SetTitleFont(132);
 382   ya->SetLabelFont(132);
 383   ya->SetTitle("dN_{ch}/d#eta");
 384
 385   if (tmp) {
 386     tmp->Draw("e5 same");
 387     all->Draw("same nostack");
 388   }
 389
 390
 391   // if (fit) fit->Draw("same");
 392   // if (tmp) tmp->Draw("same");
 393
 394   BuildCentLegend(centAxis);
 395   BuildLegend(all, centAxis);
 396
 397   AddInformation(forward);
 398
 399   c->SaveAs("dndeta_simple.C");
 400   c->SaveAs("dndeta_simple.png");
 401   c->SaveAs("dndeta_simple.root");
 402 }
 403 //
 404 // EOF
 405 //