PWGLF/FORWARD/analysis2/qa/DrawSteps.C

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