PWGLF/FORWARD/analysis2/scripts/ExtractData.C

   1 // --- Get a single histogram ----------------------------------------
   2 TH1*
   3 GetOne(const TList* list, const char* which, bool mirror)
   4 {
   5   if (!list) {
   6     Error("GetOne", "No list passed");
   7     return 0;
   8   }
   9   TString n(Form("dndeta%s_rebin05", which));
  10   if (mirror) n.Append("_mirror");
  11
  12   TObject* o = list->FindObject(n);
  13   if (!o) {
  14     Error("GetOne", "Object %s not found in %s", n.Data(), list->GetName());
  15     return 0;
  16   }
  17   TH1* ret = static_cast<TH1*>(o);
  18   ret->SetLineColor(ret->GetMarkerColor());
  19   ret->SetDirectory(0);
  20   return ret;
  21 }
  22
  23 // --- Make point-to-point systematic errors -------------------------
  24 TH1*
  25 MakeSysError(TH1* h, Double_t sysErr)
  26 {
  27   TString n(h->GetName());
  28   n.Append("_syserr");
  29   TH1* ret = static_cast<TH1*>(h->Clone(n));
  30   ret->SetMarkerStyle(1);
  31   ret->SetMarkerSize(0);
  32   ret->SetMarkerColor(kBlue-10);
  33   ret->SetFillColor(kBlue-10);
  34   ret->SetFillStyle(1001);
  35   ret->SetLineColor(kBlue-10);
  36   ret->SetLineWidth(0);
  37   ret->SetDirectory(0);
  38
  39   for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
  40     Double_t c = ret->GetBinContent(i);
  41     if (c < 0.001) {
  42       ret->SetBinContent(i,0);
  43       ret->SetBinError(i,0);
  44     }
  45     Double_t e = c * sysErr/100;
  46     ret->SetBinError(i,e);
  47   }
  48
  49   return ret;
  50 }
  51
  52 // --- Turn a TGraphAsymmErrors into a histogram ---------------------
  53 TH1* Graph2Hist(const TGraphAsymmErrors* g)
  54 {
  55   Int_t    nBins = g->GetN();
  56   TArrayF  bins(nBins+1);
  57   TArrayF  y(nBins);
  58   TArrayF  ey(nBins);
  59   Double_t dx = 0;
  60   Double_t xmin = 10000;
  61   Double_t xmax = -10000;
  62   for (Int_t i = 0; i < nBins; i++) {
  63     Double_t x   = g->GetX()[i];
  64     Double_t exl = g->GetEXlow()[i];
  65     Double_t exh = g->GetEXhigh()[i];
  66     xmin             = TMath::Min(x-exl, xmin);
  67     xmax             = TMath::Max(x+exh, xmax);
  68     bins.fArray[i]   = x-exl;
  69     bins.fArray[i+1] = x+exh;
  70     Double_t dxi = exh+exl;
  71     if (dxi == 0 && i != 0) dxi = bins.fArray[i]-bins.fArray[i-1];
  72     if (dx == 0) dx  = dxi;
  73     else if (dxi != dx) dx = 0;
  74
  75     y.fArray[i]  = g->GetY()[i];
  76     ey.fArray[i] = TMath::Max(g->GetEYlow()[i],g->GetEYhigh()[i]);
  77
  78   }
  79   TString name(g->GetName());
  80   TString title(g->GetTitle());
  81   TH1D* h = 0;
  82   if (dx != 0) {
  83     h = new TH1D(name.Data(), title.Data(), nBins,
  84                  bins[0]-dx/2, bins[nBins]+dx/2);
  85   }
  86   else {
  87     h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
  88   }
  89   for (Int_t i = 1; i <= nBins; i++) {
  90     h->SetBinContent(i, y.fArray[i-1]);
  91     h->SetBinError(i, ey.fArray[i-1]);
  92   }
  93   h->SetMarkerStyle(g->GetMarkerStyle());
  94   h->SetMarkerColor(g->GetMarkerColor());
  95   h->SetMarkerSize(g->GetMarkerSize());
  96   h->SetDirectory(0);
  97
  98   return h;
  99 }
 100
 101 // --- Set Histogram properties --------------------------------------
 102 void SetAttributes(TH1* h, UShort_t sNN, Int_t which, Bool_t mirror)
 103 {
 104   Int_t marker = (sNN == 900 ? 20     : 21) + (mirror ? 4 : 0);
 105   Int_t color  = (which == 0 ? kRed+2  :
 106                   which == 1 ? kMagenta+2 :
 107                   kBlue + 2);
 108   h->SetMarkerColor(color);
 109   h->SetLineColor(color);
 110   h->SetMarkerStyle(marker);
 111 }
 112
 113
 114 // --- Get publlished data via script --------------------------------
 115 TH1* GetPublished(UShort_t sNN, Bool_t isNSD)
 116 {
 117   Long_t ptr = 0;
 118   Int_t  typ = (isNSD ? 4 : 1);
 119
 120   ptr = gROOT->ProcessLine(Form("GetSingle(2,1,%d,%d)", sNN, typ));
 121   if (!ptr) return 0;
 122
 123   TGraphAsymmErrors* g = reinterpret_cast<TGraphAsymmErrors*>(ptr);
 124   TH1* h = Graph2Hist(g);
 125   h->SetLineColor(h->GetMarkerColor());
 126
 127   return h;
 128 }
 129
 130 // --- Add points to stack -------------------------------------------
 131 Double_t
 132 AddToStack(THStack* stack, UShort_t sNN, Bool_t isNSD,
 133            TList* forward, TList* central,
 134            Double_t fwdSysErr, Double_t cenSysErr,
 135            Double_t strangeCorr)
 136 {
 137   TH1* fd1 = GetOne(forward, "Forward", false);
 138   TH1* fd2 = GetOne(forward, "Forward", true);
 139   fd1->Scale(1-strangeCorr/100);
 140   fd2->Scale(1-strangeCorr/100);
 141   TH1* fs1 = MakeSysError(fd1, fwdSysErr);
 142   TH1* fs2 = MakeSysError(fd2, fwdSysErr);
 143
 144   TH1* cd  = 0;
 145   TH1* cs  = 0;
 146   if (central) {
 147     cd = GetOne(central, "Central", false);
 148     cs = MakeSysError(cd, cenSysErr);
 149   }
 150   else {
 151     cd = GetPublished(sNN, isNSD);
 152   }
 153   SetAttributes(fd1, sNN, 0, false);
 154   SetAttributes(fd2, sNN, 0, true);
 155   if (cd) SetAttributes(cd,  sNN, central ? 1 : 2, false);
 156
 157   if (cs) stack->Add(cs,  "e2");
 158   stack->Add(fs1, "e2");
 159   stack->Add(fs2, "e2");
 160   if (cd) stack->Add(cd,  "ep");
 161   stack->Add(fd1, "ep");
 162   stack->Add(fd2, "ep");
 163
 164   Double_t mcd = (cd ? cd->GetMaximum() : 0);
 165   Double_t mfs = fs1->GetMaximum();
 166
 167   return TMath::Max(mcd, mfs);
 168 }
 169
 170 // --- Find a list in directory --------------------------------------
 171 TList*
 172 GetList(const TDirectory* d, const char* what)
 173 {
 174   if (!d) {
 175     Error("GetList", "No diretory passed");
 176     return 0;
 177   }
 178   TList* p = static_cast<TList*>(d->Get(Form("%sResults", what)));
 179   if (!p) {
 180     Error("GetList", "%sResults not found in %s", what, d->GetName());
 181     return 0;
 182   }
 183   TList* r = static_cast<TList*>(p->FindObject("all"));
 184   if (!r) {
 185     Error("GetList", "all not found in %s", p->GetName());
 186     return 0;
 187   }
 188   return r;
 189 }
 190
 191 void
 192 AdjustAxis(TAxis* a)
 193 {
 194   a->SetTitleFont(132);
 195   a->SetLabelFont(132);
 196   a->SetTitleSize(0.08);
 197   a->SetLabelSize(0.08);
 198   a->SetTitleOffset(0.5);
 199   a->SetNdivisions(10);
 200 }
 201
 202 // --- Make a stack for a single plot --------------------------------
 203 THStack*
 204 MakeStack(Bool_t isNSD,  Bool_t showClusters,
 205           Double_t strangeCorr, Double_t maxFactor=1.1)
 206 {
 207   const char* trg = isNSD ? "nsd" : "inel";
 208   TFile* f0900 =TFile::Open(Form("forward_dndeta_%s%04d.root",trg,900, "READ"));
 209   TFile* f7000 =TFile::Open(Form("forward_dndeta_%s%04d.root",trg,7000,"READ"));
 210   if (!f0900 || !f7000) {
 211     Error("MakeStack", "Failed to open one or more files (%p, %p)",
 212           f0900, f7000);
 213     return 0;
 214   }
 215   TList* forward0900 = GetList(f0900, "Forward");
 216   TList* forward7000 = GetList(f7000, "Forward");
 217   TList* central0900 = 0;
 218   TList* central7000 = showClusters ? GetList(f7000, "Central") : 0;
 219
 220   THStack* stack = new THStack("stack", "Stack");
 221   Double_t sysDen = 5;
 222   Double_t sysPt  = 2;
 223   Double_t sysMix = 2;
 224   Double_t sysNch = 5;
 225   Double_t sysNor = 2;
 226   Double_t fwdSys = TMath::Sqrt(sysDen * sysDen +
 227                                 sysPt  * sysPt  +
 228                                 sysMix * sysMix +
 229                                 sysNch * sysNch +
 230                                 sysNor * sysNor);
 231   Info("", "Forward systematic error: %4.1f", fwdSys);
 232   Double_t m0900 =
 233     AddToStack(stack, 900,  isNSD, forward0900, central0900,
 234                fwdSys, 5, strangeCorr);
 235   Double_t m7000 =
 236     AddToStack(stack, 7000, isNSD, forward7000, central7000,
 237                fwdSys, 5, strangeCorr);
 238
 239   stack->SetMaximum(maxFactor * TMath::Max(m0900, m7000));
 240
 241   f0900->Close();
 242   f7000->Close();
 243
 244   return stack;
 245 }
 246
 247 // --- Draw one stack ------------------------------------------------
 248 Double_t
 249 DrawStack(Bool_t isNSD, Bool_t showClusters,
 250           Double_t strangeCorr, Double_t maxFactor=1.1)
 251 {
 252   THStack* stack = MakeStack(isNSD, showClusters, strangeCorr, maxFactor);
 253   stack->Draw("nostack ep");
 254   stack->GetHistogram()->SetXTitle("#eta");
 255   if (!isNSD)
 256     stack->GetHistogram()->SetYTitle("#frac{1}{N} #frac{dN_{ch}}{d#eta}");
 257   if (!isNSD)
 258     stack->SetMinimum(.0001);
 259   AdjustAxis(stack->GetHistogram()->GetXaxis());
 260   AdjustAxis(stack->GetHistogram()->GetYaxis());
 261   stack->GetHistogram()->GetXaxis()->SetTitleOffset(1);
 262   gPad->Clear();
 263   stack->Draw("nostack ep");
 264
 265   TLatex* title = new TLatex(1-gPad->GetRightMargin()-.01,
 266                              1-gPad->GetTopMargin()-.01,
 267                              (isNSD ? "NSD" : "INEL"));
 268   title->SetNDC();
 269   title->SetTextSize(0.1);
 270   title->SetTextFont(132);
 271   title->SetTextAlign(33);
 272   title->SetTextColor(kBlue+2);
 273   title->Draw();
 274
 275   return stack->GetMaximum();
 276 }
 277
 278 // --- Draw one stack ------------------------------------------------
 279 void
 280 ExtractData(Bool_t showClusters = true)
 281 {
 282   gROOT->LoadMacro("OtherData.C");
 283   gStyle->SetOptTitle(0);
 284   gStyle->SetGridColor(kGray);
 285   TCanvas* c = new TCanvas("c", "C", 1200, 850);
 286   c->SetRightMargin(0.01);
 287   c->SetLeftMargin(0.1);
 288   c->SetTopMargin(0.02);
 289   c->SetBottomMargin(0.15);
 290   c->SetFillColor(0);
 291   c->SetBorderSize(0);
 292   c->SetBorderMode(0);
 293   c->cd();
 294
 295   c->Divide(1, 2, 0, 0);
 296
 297   TVirtualPad* p = c->cd(1);
 298   p->SetFillColor(0);
 299   p->SetRightMargin(0.02);
 300   p->SetGridx();
 301   p->SetGridy();
 302   DrawStack(false, showClusters, 1.5);
 303
 304   Double_t ms = 1.2;
 305
 306   TLegend* l = new TLegend(.33, .01, .53, .35);
 307   l->SetBorderSize(0);
 308   l->SetFillColor(0);
 309   l->SetTextFont(132);
 310   l->SetFillStyle(0);
 311   TLegendEntry* e = l->AddEntry("", "900GeV", "p");
 312   e->SetMarkerStyle(20);
 313   e->SetMarkerSize(ms+.1);
 314   e = l->AddEntry("", "    7TeV", "p");
 315   e->SetMarkerStyle(21);
 316   e->SetMarkerSize(ms);
 317   e = l->AddEntry("", "Mirrored data", "p");
 318   e->SetMarkerStyle(24);
 319   e->SetMarkerSize(ms);
 320   l->Draw();
 321
 322   TLegend* l2 = new TLegend(.5, .01, .8, .35);
 323   l2->SetBorderSize(0);
 324   l2->SetFillColor(0);
 325   l2->SetTextFont(132);
 326   l2->SetFillStyle(0);
 327   l2->SetColumnSeperation(-.05);
 328   e = l2->AddEntry("", "Forward", "p");
 329   e->SetMarkerStyle(20);
 330   e->SetMarkerColor(kRed+2);
 331   e->SetLineColor(kRed+2);
 332   e->SetMarkerSize(ms);
 333   if (showClusters) {
 334     e = l2->AddEntry("", "Central", "p");
 335     e->SetMarkerStyle(20);
 336     e->SetMarkerColor(kMagenta+2);
 337     e->SetLineColor(kMagenta+2);
 338     e->SetMarkerSize(ms);
 339   }
 340   e = l2->AddEntry("", "#splitline{Eur.Phys.J.#font[22]{C68}:89-108}"
 341                    "{Eur.Phys.J.#font[22]{C68}:345--354}", "p");
 342   e->SetMarkerStyle(20);
 343   e->SetMarkerColor(kBlue+2);
 344   e->SetLineColor(kBlue+2);
 345   e->SetMarkerSize(ms);
 346   l2->Draw();
 347
 348
 349   p = c->cd(2);
 350   p->SetGridx();
 351   p->SetGridy();
 352   p->SetFillColor(0);
 353   p->SetRightMargin(0.02);
 354   DrawStack(true, showClusters, 1.5);
 355
 356   // TPad, x1, y1, x2, y2, ts, t1, t2, prel
 357   gROOT->LoadMacro("AddLogo.C");
 358   AddLogo((TPad*)p, .4, p->GetBottomMargin()+.05,
 359           .5, .44, 0.08, "", "pp data", true);
 360
 361   c->cd();
 362   TString out("dndeta_pp_forward");
 363   if (!showClusters) out.Append("_noclusters");
 364   c->SaveAs(Form("%s.png", out.Data()));
 365   c->SaveAs(Form("%s.eps", out.Data()));
 366 }
 367
 368
 369
 370