PWG2/FORWARD/analysis2/scripts/DrawMCResult.C

   1 THStack*
   2 GetStack(const TList& forward,  const char* sub, const char* name)
   3 {
   4   TList* lsub = static_cast<TList*>(forward.FindObject(sub));
   5   if (!lsub) {
   6     Warning("GetStack", "Sub list %s not found in %s", sub, forward.GetName());
   7     return 0;
   8   }
   9   THStack* ret = static_cast<THStack*>(lsub->FindObject(name));
  10   if (!ret)
  11     Warning("GetStack" "Stack %s not found in %s", name, sub);
  12   return ret;
  13 }
  14
  15 TH1*
  16 Rebin(TH1* h, Int_t rebin)
  17 {
  18   if (rebin <= 1) return h;
  19   h->Rebin(rebin);
  20   h->Scale(1. / rebin);
  21   return h;
  22 }
  23
  24 TH1*
  25 Ratio(const TH1* h1, const TH1* h2)
  26 {
  27   if (!h1) return;
  28   if (!h2) return;
  29
  30   TH1* copy = static_cast<TH1*>(h2->Clone("tmp"));
  31   copy->SetName(Form("%s_%s", h2->GetName(), h1->GetName()));
  32   copy->SetTitle(Form("%s/%s", h2->GetTitle(), h1->GetTitle()));
  33   copy->SetDirectory(0);
  34   copy->Divide(h1);
  35
  36   return copy;
  37 }
  38
  39 Int_t
  40 Ratio(THStack* r, const THStack* h1, const THStack* h2)
  41 {
  42   if (!h1) return 0;
  43   if (!h2) return 0;
  44
  45   int n1 = h1->GetHists()->GetEntries();
  46   int n2 = h2->GetHists()->GetEntries();
  47   int nH = 0;
  48   for (int i = 0; i < n1 && i < n2; i++) {
  49     TH1* hh1 = static_cast<TH1*>(h1->GetHists()->At(i));
  50     TH1* hh2 = static_cast<TH1*>(h2->GetHists()->At(i));
  51     TH1* h   = Ratio(hh1, hh2);
  52     if (!h) continue;
  53     nH++;
  54     r->Add(h);
  55   }
  56   return nH;
  57 }
  58
  59 void
  60 DrawMCResult(const char* filename="forward.root", Int_t rebin=1,
  61              Bool_t ratios=true)
  62 {
  63   gStyle->SetPalette(1);
  64   gStyle->SetOptFit(0);
  65   gStyle->SetOptStat(0);
  66
  67   TFile* file = TFile::Open(filename, "READ");
  68   if (!file) {
  69     Error("DrawMCResult", "failed to open %s", filename);
  70     return;
  71   }
  72
  73   TList* forward = static_cast<TList*>(file->Get("Forward"));
  74   if (!forward) {
  75     Error("DrawMCResult", "List Forward not found in %s", filename);
  76     return;
  77   }
  78   THStack* res    = GetStack(*forward, "ringResults", "all");
  79   THStack* mcRes  = GetStack(*forward, "mcRingResults", "all");
  80   THStack* deltas = GetStack(*forward, "fmdSharingFilter", "sums");
  81   THStack* nchs   = GetStack(*forward, "fmdDensityCalculator", "sums");
  82   THStack* prims  = GetStack(*forward, "fmdCorrector", "sums");
  83
  84   TH1* sumEta = static_cast<TH1*>(forward->FindObject("mcSumEta"));
  85   if (!sumEta) {
  86     Warning("DrawMCResults", "mcSumEta not found in Forward");
  87   }
  88
  89
  90   gStyle->SetTitleBorderSize(0);
  91   gStyle->SetTitleFillColor(0);
  92   gStyle->SetTitleStyle(0);
  93   // gStyle->SetTitleColor(kBlack);
  94
  95
  96   TCanvas* c = new TCanvas("c", "C", 900, 700);
  97   c->SetFillColor(0);
  98   c->SetBorderSize(0);
  99   c->SetTopMargin(0.05);
 100   c->SetRightMargin(0.05);
 101
 102   Double_t y1 = (ratios ? .3 : 0);
 103   TPad* p1 = new TPad("p1", "p1", 0, y1, 1, 1, 0, 0, 0);
 104   p1->SetBottomMargin(ratios ? 0.01 : .10);
 105   p1->SetFillColor(0);
 106   p1->SetBorderSize(0);
 107   p1->SetTopMargin(0.05);
 108   p1->SetRightMargin(0.05);
 109   p1->Draw();
 110   p1->cd();
 111
 112   THStack* all = new THStack("all", "Analysis steps");
 113   if (res) {
 114     res->SetTitle("dN_{ch}/d#eta");
 115     TIter next(res->GetHists());
 116     TH1* h = 0;
 117     while ((h = static_cast<TH1*>(next()))) all->Add(Rebin(h,rebin));
 118   }
 119   if (mcRes) {
 120     mcRes->SetTitle("Track-Refs");
 121     TIter next(mcRes->GetHists());
 122     TH1* h = 0;
 123     while ((h = static_cast<TH1*>(next()))) all->Add(Rebin(h,rebin));
 124   }
 125   if (deltas) {
 126     deltas->SetTitle("#sum_{} #Delta/#Delta_{mip}");
 127     TIter next(deltas->GetHists());
 128     TH1* h = 0;
 129     while ((h = static_cast<TH1*>(next()))) {
 130       h->SetMarkerStyle(25);
 131       all->Add(Rebin(h,rebin));
 132     }
 133   }
 134   if (nchs) {
 135     nchs->SetTitle("#sum_{} N_{ch,incl}");
 136     TIter next(nchs->GetHists());
 137     TH1* h = 0;
 138     while ((h = static_cast<TH1*>(next()))) {
 139       h->SetMarkerStyle(21);
 140       all->Add(Rebin(h,rebin));
 141     }
 142   }
 143   if (prims) {
 144     prims->SetTitle("#sum_{} N_{ch,primary}");
 145     TIter next(prims->GetHists());
 146     TH1* h = 0;
 147     while ((h = static_cast<TH1*>(next()))) {
 148       h->SetMarkerStyle(22);
 149       all->Add(Rebin(h,rebin));
 150     }
 151   }
 152   if (sumEta) all->Add(sumEta);
 153   all->Draw("nostack");
 154   all->GetHistogram()->SetXTitle("#eta");
 155   all->GetHistogram()->SetYTitle("signal");
 156   all->GetHistogram()->GetXaxis()->SetLabelFont(132);
 157   all->GetHistogram()->GetXaxis()->SetTitleFont(132);
 158   all->GetHistogram()->GetYaxis()->SetLabelFont(132);
 159   all->GetHistogram()->GetYaxis()->SetTitleFont(132);
 160   c->SetGridx();
 161
 162   TLegend* l = new TLegend(.33, .2, .53, .9);
 163   TLegendEntry* e = 0;
 164   l->SetFillColor(0);
 165   l->SetFillStyle(0);
 166   l->SetBorderSize(0);
 167   l->SetNColumns(1);
 168   l->SetTextFont(132);
 169   Int_t i = 0;
 170   if (res) {
 171     TIter next(res->GetHists());
 172     TH1* h = 0;
 173     while ((h = static_cast<TH1*>(next()))) {
 174       e = l->AddEntry(Form("dummy%02d", i++),h->GetTitle(),"pl");
 175       e->SetMarkerStyle(20);
 176       e->SetMarkerColor(h->GetMarkerColor());
 177     }
 178     e = l->AddEntry(Form("dummy%02d", i++),res->GetTitle(),"pl");
 179     e->SetMarkerStyle(20);
 180     e->SetMarkerColor(kBlack);
 181   }
 182   if (deltas) {
 183     e = l->AddEntry(Form("dummy%02d", i++), deltas->GetTitle(),"pl");
 184     TH1* h = static_cast<TH1*>(deltas->GetHists()->At(0));
 185     e->SetMarkerStyle(h->GetMarkerStyle());
 186     e->SetMarkerColor(kBlack);
 187   }
 188   if (nchs) {
 189     e = l->AddEntry(Form("dummy%02d",i++),nchs->GetTitle(),"pl");
 190     TH1* h = static_cast<TH1*>(nchs->GetHists()->At(0));
 191     e->SetMarkerStyle(h->GetMarkerStyle());
 192     e->SetMarkerColor(kBlack);
 193   }
 194   if (prims) {
 195     e = l->AddEntry(Form("dummy%02d", i++), prims->GetTitle(),"pl");
 196     TH1* h = static_cast<TH1*>(prims->GetHists()->At(0));
 197     e->SetMarkerStyle(h->GetMarkerStyle());
 198     e->SetMarkerColor(kBlack);
 199   }
 200
 201   if (mcRes) {
 202     e = l->AddEntry(Form("dummy%02d", i++), mcRes->GetTitle(), "pl");
 203     TH1* h = static_cast<TH1*>(mcRes->GetHists()->At(0));
 204     e->SetMarkerStyle(h->GetMarkerStyle());
 205     e->SetMarkerColor(kBlack);
 206   }
 207
 208   if (sumEta) l->AddEntry(sumEta);
 209   l->Draw();
 210
 211
 212   if (!ratios) return;
 213
 214   c->cd();
 215   TPad* p2 = new TPad("p2", "p2", 0, 0, 1, y1, 0, 0, 0);
 216   p2->SetTopMargin(0);
 217   p2->SetFillColor(0);
 218   p2->SetBorderSize(0);
 219   p2->SetRightMargin(0.05);
 220   p2->Draw();
 221   p2->cd();
 222
 223   THStack* rs = new THStack("ratios", "Ratios");
 224   Int_t nDN = Ratio(rs, deltas, nchs);
 225   Int_t nNR = Ratio(rs, nchs, res);
 226   Int_t nRP = Ratio(rs, res, prims);
 227   rs->Draw("nostack");
 228
 229   TLegend* ll = new TLegend(.38, .2, .48, .9);
 230   ll->SetFillColor(0);
 231   ll->SetFillStyle(0);
 232   ll->SetBorderSize(0);
 233   ll->SetNColumns(1);
 234   ll->SetTextFont(132);
 235   if (nDN) {
 236     e = ll->AddEntry("d1", Form("#frac{%s}{%s}",
 237                                               nchs->GetTitle(),
 238                                               deltas->GetTitle()), "lp");
 239     e->SetMarkerStyle(21);
 240   }
 241   if (nNR) {
 242     e = ll->AddEntry("d2", Form("#frac{%s}{%s}", res->GetTitle(),
 243                                 nchs->GetTitle()), "lp");
 244     e->SetMarkerStyle(20);
 245   }
 246   if (nRP) {
 247     e = ll->AddEntry("d3", Form("#frac{%s}{%s}", prims->GetTitle(),
 248                             res->GetTitle()), "lp");
 249     e->SetMarkerStyle(22);
 250   }
 251   ll->Draw();
 252
 253
 254 }
 255
 256
 257
 258