PWG2/FORWARD/analysis2/qa/DrawSteps.C

   1 /**
   2  *
   3  * @defgroup pwg2_forward_scripts_qa Quality Assurance scripts
   4  *
   5  *
   6  * @ingroup pwg2_forward_scripts
   7  *
   8  */
   9
  10 /**
  11  * Get a stack
  12  *
  13  * @param forward   Input list
  14  * @param sub       Sub-list
  15  * @param name      Name of stack
  16  *
  17  * @return A stack or null
  18  *
  19  * @ingroup pwg2_forward_scripts_qa
  20  */
  21 THStack*
  22 GetStack(const TList& forward,  const char* sub, const char* name)
  23 {
  24   TList* lsub = &forward;
  25
  26   if (sub && sub[0] != '\0')
  27     lsub = static_cast<TList*>(forward.FindObject(sub));
  28
  29   if (!lsub) {
  30     Warning("GetStack", "Sub list %s not found in %s", sub, forward.GetName());
  31     return 0;
  32   }
  33   THStack* ret = static_cast<THStack*>(lsub->FindObject(name));
  34   if (!ret)
  35     Warning("GetStack" "Stack %s not found in %s", name, sub);
  36   return ret;
  37 }
  38
  39 /**
  40  * Rebin a histogram
  41  *
  42  * @param h      Histogram
  43  * @param rebin  Rebinning factor
  44  *
  45  * @return Histogram
  46  *
  47  * @ingroup pwg2_forward_scripts_qa
  48  */
  49 TH1*
  50 Rebin(TH1* h, Int_t rebin)
  51 {
  52   if (rebin <= 1) return h;
  53   h->Rebin(rebin);
  54   h->Scale(1. / rebin);
  55   return h;
  56 }
  57
  58 /**
  59  * Ratio of two histograms
  60  *
  61  * @param h1 numerator
  62  * @param h2 denominator
  63  *
  64  * @return Ratio
  65  *
  66  * @ingroup pwg2_forward_scripts_qa
  67  */
  68 TH1*
  69 Ratio(const TH1* h1, const TH1* h2)
  70 {
  71   if (!h1) return;
  72   if (!h2) return;
  73
  74   TH1* copy = static_cast<TH1*>(h2->Clone("tmp"));
  75   copy->SetName(Form("%s_%s", h2->GetName(), h1->GetName()));
  76   copy->SetTitle(Form("%s/%s", h2->GetTitle(), h1->GetTitle()));
  77   copy->SetDirectory(0);
  78   copy->Divide(h1);
  79
  80   return copy;
  81 }
  82
  83 /**
  84  * Ratio all histograms in stacks
  85  *
  86  * @param r  Result
  87  * @param h1 Numerators
  88  * @param h2 Denominators
  89  *
  90  * @return Number of histograms
  91   *
  92  * @ingroup pwg2_forward_scripts_qa
  93 */
  94 Int_t
  95 Ratio(THStack* r, const THStack* h1, const THStack* h2)
  96 {
  97   if (!h1) return 0;
  98   if (!h2) return 0;
  99
 100   int n1 = h1->GetHists()->GetEntries();
 101   int n2 = h2->GetHists()->GetEntries();
 102   int nH = 0;
 103   for (int i = 0; i < n1 && i < n2; i++) {
 104     TH1* hh1 = static_cast<TH1*>(h1->GetHists()->At(i));
 105     TH1* hh2 = static_cast<TH1*>(h2->GetHists()->At(i));
 106     TH1* h   = Ratio(hh1, hh2);
 107     if (!h) continue;
 108     nH++;
 109     r->Add(h);
 110   }
 111   return nH;
 112 }
 113
 114 /**
 115  * Add a histogram to the all stack
 116  *
 117  * @param all         Stack
 118  * @param h           Histogram
 119  * @param singleStep  Showing individual steps?
 120  *
 121  * @ingroup pwg2_forward_scripts_qa
 122 */
 123 void
 124 AddToAll(THStack* all, const TH1* h, Bool_t singleStep)
 125 {
 126   TH1* copy = static_cast<TH1*>(h->Clone(Form("%s_copy", h->GetName())));
 127   copy->SetDirectory(0);
 128   if (singleStep) {
 129     copy->SetMarkerColor(kGray);
 130     copy->SetLineColor(kGray);
 131   }
 132   all->Add(copy);
 133 }
 134
 135 /**
 136  * Dim an entry
 137  *
 138  * @param thisId  This step
 139  * @param step    Current step
 140  * @param e       Entry in legend
 141  *
 142  * @ingroup pwg2_forward_scripts_qa
 143  */
 144 void
 145 DimEntry(Int_t thisId, Int_t step, TLegendEntry* e)
 146 {
 147
 148   Int_t col = (thisId == step || step <= 0) ? kBlack : kGray;
 149   e->SetMarkerColor(col);
 150   e->SetLineColor(col);
 151   e->SetTextColor(col);
 152 }
 153
 154 /**
 155  * Draw a step
 156  *
 157  * @param deltas   From energy loss
 158  * @param nchs     After 2nd correction
 159  * @param prims    Primaries
 160  * @param dndeta   Result
 161  * @param step     Step number
 162  *
 163  * @ingroup pwg2_forward_scripts_qa
 164  */
 165 void
 166 DrawStep(THStack* deltas, THStack* nchs, THStack* prims,
 167          TH1*     dndeta, Int_t step)
 168 {
 169   THStack* all = new THStack("all", "Analysis steps");
 170   if (step > 0) all->SetTitle(Form("Step %d", step));
 171
 172   if (deltas) {
 173     deltas->SetTitle("#sum_{} #Delta/#Delta_{mip}");
 174     TIter next(deltas->GetHists());
 175     TH1* h = 0;
 176     while ((h = static_cast<TH1*>(next()))) {
 177       h->SetMarkerStyle(25);
 178       // Info("DrawStep", "Adding %s", h->GetName());
 179       AddToAll(all, h, step>0);
 180     }
 181   }
 182   if (nchs) {
 183     nchs->SetTitle("#sum_{} N_{ch,incl}");
 184     TIter next(nchs->GetHists());
 185     TH1* h = 0;
 186     while ((h = static_cast<TH1*>(next()))) {
 187       h->SetMarkerStyle(21);
 188       // Info("DrawStep", "Adding %s", h->GetName());
 189       AddToAll(all, h, step>0);
 190     }
 191   }
 192   if (prims) {
 193     prims->SetTitle("#sum_{} N_{ch,primary}");
 194     TIter next(prims->GetHists());
 195     TH1* h = 0;
 196     while ((h = static_cast<TH1*>(next()))) {
 197       h->SetMarkerStyle(22);
 198       // Info("DrawStep", "Adding %s", h->GetName());
 199       AddToAll(all, h, step>0);
 200     }
 201   }
 202   if (dndeta) {
 203     dndeta->SetTitle("1/N dN_{ch}/d#eta");
 204     dndeta->SetMarkerStyle(20);
 205     dndeta->SetMarkerColor(kBlack);
 206     // Info("DrawStep", "Adding %s", dndeta->GetName());
 207     AddToAll(all, dndeta, step>0);
 208   }
 209
 210   all->Draw("nostack");
 211   all->GetHistogram()->SetXTitle("#eta");
 212   all->GetHistogram()->SetYTitle("signal");
 213   all->GetHistogram()->GetXaxis()->SetLabelFont(132);
 214   all->GetHistogram()->GetXaxis()->SetTitleFont(132);
 215   all->GetHistogram()->GetYaxis()->SetLabelFont(132);
 216   all->GetHistogram()->GetYaxis()->SetTitleFont(132);
 217   c->SetGridx();
 218
 219   TLegend* l = new TLegend(.33, .2, .53, .9);
 220   TLegendEntry* e = 0;
 221   l->SetFillColor(0);
 222   l->SetFillStyle(0);
 223   l->SetBorderSize(0);
 224   l->SetNColumns(1);
 225   l->SetTextFont(132);
 226   Int_t i = 0;
 227   if (deltas) {
 228     TIter next(deltas->GetHists());
 229     TH1*  h = 0;
 230     while ((h = static_cast<TH1*>(next()))) {
 231       e = l->AddEntry(Form("dummy%02d", i++),h->GetTitle(),"pl");
 232       e->SetMarkerStyle(20);
 233       e->SetMarkerColor(h->GetMarkerColor());
 234     }
 235     e = l->AddEntry(Form("dummy%02d", i++), deltas->GetTitle(),"pl");
 236     TH1* h = static_cast<TH1*>(deltas->GetHists()->At(0));
 237     e->SetMarkerStyle(h->GetMarkerStyle());
 238     DimEntry(1, step, e);
 239   }
 240   if (nchs) {
 241     e = l->AddEntry(Form("dummy%02d",i++),nchs->GetTitle(),"pl");
 242     TH1* h = static_cast<TH1*>(nchs->GetHists()->At(0));
 243     e->SetMarkerStyle(h->GetMarkerStyle());
 244     DimEntry(2, step, e);
 245   }
 246   if (prims) {
 247     e = l->AddEntry(Form("dummy%02d", i++), prims->GetTitle(),"pl");
 248     TH1* h = static_cast<TH1*>(prims->GetHists()->At(0));
 249     e->SetMarkerStyle(h->GetMarkerStyle());
 250     DimEntry(3, step, e);
 251   }
 252   if (dndeta) {
 253     e = l->AddEntry(Form("dummy%02d", i++), dndeta->GetTitle(),"pl");
 254     e->SetMarkerStyle(dndeta->GetMarkerStyle());
 255     DimEntry(4, step, e);
 256   }
 257   l->Draw();
 258
 259   TString what;
 260   if (step > 0) {
 261     switch (step) {
 262     case 1:
 263       deltas->Draw("same nostack");
 264       what = "After merging";
 265       break;
 266     case 2:
 267       nchs->Draw("same nostack");
 268       what = "After particle counting";
 269       break;
 270     case 3:
 271       prims->Draw("same nostack");
 272       what = "After corrections";
 273       break;
 274     case 4:
 275       dndeta->Draw("same");
 276       what = "After normalisation";
 277       break;
 278     default:
 279       Error("DrawSteps", "Unknown step: %d (must be in 1-4)");
 280       break;
 281     }
 282   }
 283   TLatex* ltx = new TLatex(.95, .85, what);
 284   ltx->SetNDC();
 285   ltx->SetTextSize(.07);
 286   ltx->SetTextAlign(33);
 287   ltx->SetTextFont(132);
 288   ltx->Draw();
 289 }
 290
 291 /**
 292  * Draw steps
 293  *
 294  * @param filename Input file
 295  * @param single   Whether to show individial steps
 296  *
 297  * @ingroup pwg2_forward_scripts_qa
 298  */
 299 void DrawSteps(const char* filename="forward.root", Bool_t single=true)
 300 {
 301   gStyle->SetPalette(1);
 302   gStyle->SetOptFit(0);
 303   gStyle->SetOptStat(0);
 304
 305   TFile* file = TFile::Open(filename, "READ");
 306   if (!file) {
 307     Error("DrawMCResult", "failed to open %s", filename);
 308     return;
 309   }
 310   const char* fname2 = "forward_dndeta.root";
 311   TFile* file2 = TFile::Open(fname2, "READ");
 312   if (!file2) {
 313     Error("DrawSteps", "File %s not found", fname2);
 314   }
 315
 316   TList* forward = static_cast<TList*>(file->Get("Forward"));
 317   if (!forward) {
 318     Error("DrawMCResult", "List Forward not found in %s", filename);
 319     return;
 320   }
 321   TList* forwardRes = (file2 ?
 322                        static_cast<TList*>(file2->Get("ForwardResults")) :
 323                        0);
 324   TList* forwardAll = (forwardRes ?
 325                        static_cast<TList*>(forwardRes->FindObject("all")) :
 326                        0);
 327
 328
 329   // THStack* res    = GetStack(*forward, "ringResults", "all");
 330   // THStack* mcRes  = GetStack(*forward, "mcRingResults", "all");
 331   THStack* deltas = GetStack(*forward, "fmdSharingFilter", "sums");
 332   THStack* nchs   = GetStack(*forward, "fmdDensityCalculator", "sums");
 333   THStack* prims  = GetStack(*forward, "fmdCorrector", "sums");
 334   TH1*     dndeta = (forwardAll ?
 335                      static_cast<TH1*>(forwardAll->FindObject("dndetaForward")):
 336                      0);
 337
 338   Info("DrawSteps", "Got steps deltas=%p, nchs=%p, prims=%p, dndeta=%p",
 339        deltas, nchs, prims, dndeta);
 340
 341
 342   gStyle->SetTitleBorderSize(0);
 343   gStyle->SetTitleFillColor(0);
 344   gStyle->SetTitleStyle(0);
 345   gStyle->SetTitleX(.7);
 346   gStyle->SetTitleY(.95);
 347   gStyle->SetTitleH(.1);
 348   gStyle->SetTitleW(.25);
 349   gStyle->SetOptTitle(1);
 350   // gStyle->SetTitleColor(kBlack);
 351
 352
 353
 354   if (!single) {
 355     TCanvas* c = new TCanvas("c", "C", 900, 700);
 356     c->SetFillColor(0);
 357     c->SetBorderSize(0);
 358     c->SetTopMargin(0.05);
 359     c->SetRightMargin(0.05);
 360
 361     DrawStep(deltas, nchs, prims, dndeta, 0);
 362     c->SaveAs("steps_all.png");
 363     return;
 364   }
 365   Int_t nSteps = 0;
 366   if (deltas) nSteps++;
 367   if (nchs)   nSteps++;
 368   if (prims)  nSteps++;
 369   if (dndeta) nSteps++;
 370
 371   Int_t w = (nSteps >= 4 ? 1100 :  700);
 372   Int_t h = (nSteps >= 4 ? 800  : 1100);
 373
 374   TCanvas* c = new TCanvas("c", "C", w, h);
 375   c->SetFillColor(0);
 376   c->SetBorderSize(0);
 377   c->SetTopMargin(0.05);
 378   c->SetRightMargin(0.05);
 379
 380   if (nSteps >= 4)
 381     c->Divide(2,(nSteps+1)/2,0,0);
 382   else
 383     c->Divide(1,nSteps,0,0);
 384
 385   for (Int_t i=1; i<=nSteps; i++) {
 386     TVirtualPad* p = c->cd(i);
 387     p->SetFillColor(0);
 388     p->SetFillStyle(0);
 389     p->SetBorderSize(0);
 390     p->SetGridx();
 391     p->SetGridy();
 392
 393     DrawStep(deltas, nchs, prims, dndeta, i);
 394   }
 395   c->SaveAs("steps_comic.png");
 396 }
 397 //
 398 // EOF
 399 //
 400
 401