PWGLF/FORWARD/analysis2/scan/Trend.C

   1 #ifndef __TREND_C_
   2 #define __TREND_C_
   3 #include "SummaryDrawer.C"
   4 #include <TH1.h>
   5 #include <THStack.h>
   6 #include <TMultiGraph.h>
   7 #include <TGraphAsymmErrors.h>
   8 #include <TTree.h>
   9 #include <TList.h>
  10 #include <TLegend.h>
  11
  12 /**
  13  * Structure to make trendind plots from cut scans
  14  *
  15  */
  16 struct Trend : public SummaryDrawer
  17 {
  18   //__________________________________________________________________
  19   /**
  20    * Constructor
  21    */
  22   Trend()
  23   {
  24     fSLCuts.SetName("sl");
  25     fSHCuts.SetName("sh");
  26     fDCCuts.SetName("dc");
  27     fOrder[0] = &fSLCuts;
  28     fOrder[1] = &fSHCuts;
  29     fOrder[2] = &fDCCuts;
  30   }
  31   //=== Customization member function ================================
  32   /**
  33    * Add a sharing filter low cut
  34    *
  35    * @param name    Name of the cut type (fix,mpv,sig,xi,prob)
  36    * @param values  Space-separated cut parameters
  37    */
  38   void AddSLCut(const TString& name, const TString& values)
  39   {
  40     fSLCuts.Add(new TNamed(name, values));
  41   }
  42   //__________________________________________________________________
  43   /**
  44    * Add a sharing filter high cut
  45    *
  46    * @param name    Name of the cut type (fix,mpv,sig,xi,prob)
  47    * @param values  Space-separated cut parameters
  48    */
  49   void AddSHCut(const TString& name, const TString& values)
  50   {
  51     fSHCuts.Add(new TNamed(name, values));
  52   }
  53   //__________________________________________________________________
  54   /**
  55    * Add a density calculator threshold
  56    *
  57    * @param name    Name of the cut type (fix,mpv,sig,xi,prob)
  58    * @param values  Space-separated cut parameters
  59    */
  60   void AddDCCut(const TString& name, const TString& values)
  61   {
  62     fDCCuts.Add(new TNamed(name, values));
  63   }
  64   //__________________________________________________________________
  65   /**
  66    * Add a run to analyse
  67    *
  68    * @param name Run number
  69    */
  70   void AddRun(const TString& name)
  71   {
  72     fRuns.Add(new TNamed(name, name));
  73   }
  74   //__________________________________________________________________
  75   /**
  76    * Add a centrality bin to analyse
  77    *
  78    * @param l  Lower bound
  79    * @param h  Upper bound
  80    */
  81   void AddCentrality(UShort_t l, UShort_t h)
  82   {
  83     TObjString* n = 0;
  84     if (l <= 0 && h >= 100) n = new TObjString("all");
  85     else                    n = new TObjString(Form("cent%03d_%03d", l, h));
  86     n->SetUniqueID((h & 0xFFFF) << 16 | (l & 0xFFFF));
  87     fCents.Add(n);
  88   }
  89   //__________________________________________________________________
  90   /**
  91    * Set the order in which to do the comparison.  Most important
  92    * should be last.
  93    *
  94    * @param order Permutation of sl, sh, and dc, space separated
  95    */
  96   void SetOrder(const TString& order)
  97   {
  98     TObjArray* a = order.Tokenize(" ");
  99     if (a->GetEntries() < 3) {
 100       Error("SetOrder", "Order must contain a pertubation of "
 101             "\"sl\", \"sh\", and \"dc\", separated by spaces");
 102       fOrder[0] = fOrder[1] = fOrder[2] = 0;
 103       return;
 104     }
 105     for (Int_t i = 0; i < 3; i++) {
 106       const TString& n = static_cast<TObjString*>(a->At(i))->String();
 107       TList* l = 0;
 108       if      (n.EqualTo("sl", TString::kIgnoreCase)) l = &fSLCuts;
 109       else if (n.EqualTo("sh", TString::kIgnoreCase)) l = &fSHCuts;
 110       else if (n.EqualTo("dc", TString::kIgnoreCase)) l = &fDCCuts;
 111       if (!l) {
 112         Error("SetOrder", "Unknown cut \"%s\"", n.Data());
 113         fOrder[0] = fOrder[1] = fOrder[2] = 0;
 114         return;
 115       }
 116       fOrder[i] = l;
 117     }
 118     a->Delete();
 119   }
 120   //=== Processing member functions ==================================
 121   /**
 122    * Run it
 123    *
 124    * @param sys
 125    * @param sNN
 126    * @param trg
 127    */
 128   void Run(const char* output="trending.root",
 129            UShort_t sys=2,
 130            UShort_t sNN=2760,
 131            UShort_t trg=1)
 132   {
 133     TString outName(output);
 134     outName.ReplaceAll(".pdf", ".root");
 135     if (!outName.EndsWith(".root")) outName.Append(".root");
 136     TString pdfName(outName);
 137     pdfName.ReplaceAll(".root", ".pdf");
 138
 139     CreateCanvas(pdfName, true);
 140
 141     TFile*      out    = TFile::Open(outName, "RECREATE");
 142     TDirectory* refDir = out->mkdir("reference");
 143
 144     TIter iCent(&fCents);
 145     TObject* pCent = 0;
 146     while ((pCent = iCent())) {
 147       UShort_t low   = (pCent->GetUniqueID() & 0xFFFF);
 148       UShort_t high  = (pCent->GetUniqueID() >> 16) & 0xFFFF;
 149       TGraph*  graph = GetOther(sys, sNN, trg, low, high);
 150       if (!graph) break;
 151
 152       TString name(graph->GetName());
 153       name.Append(Form("_%s", pCent->GetName()));
 154       graph->SetName(name);
 155
 156       fOthers.Add(graph);
 157       refDir->Add(graph);
 158     }
 159
 160     TIter    iRun(&fRuns);
 161     TObject* pRun = 0;
 162     while ((pRun = iRun())) {
 163       TString     sRun = pRun->GetName();
 164       TDirectory* dRun = out->mkdir(sRun);
 165
 166       MakeChapter(sRun);
 167
 168       for (Int_t i = 0; i < 1/*2*/; i++) {
 169         TString     three(Form("%dstrip", i+2));
 170         TDirectory* dStrip = dRun->mkdir(three);
 171
 172         MakeChapter(Form("_%s", three.Data()));
 173
 174         three.Append("_slX_shX_dcX");
 175
 176         LoopCuts(sRun, three, fOrder[0], fOrder[1], fOrder[2], dStrip);
 177       }
 178     }
 179     out->Write();
 180     out->Close();
 181     CloseCanvas();
 182   }
 183   //__________________________________________________________________
 184   /**
 185    * Loop over cuts.  Recursively calls self.
 186    *
 187    * @param run      Run number
 188    * @param cur      Current path
 189    * @param cuts1    Cuts
 190    * @param cuts2    Cuts (or null)
 191    * @param cuts3    Cuts (or null)
 192    * @param out      Output diretory
 193    * @param upBin    Parents bin number
 194    * @param upMean   Parents mean
 195    * @param upRatios Parents ratios
 196    */
 197   void LoopCuts(const TString& run,
 198                 const TString& cur,
 199                 const TList*   cuts1,
 200                 const TList*   cuts2,
 201                 const TList*   cuts3,
 202                 TDirectory*    out,
 203                 Int_t          upBin=0,
 204                 THStack*       upMean=0,
 205                 THStack*       upRatios=0)
 206   {
 207     TString pre(Form("%s%s",(cuts3 ? "" : "_"), (cuts2 ? "" : "_")));
 208     Printf("%s%s", pre.Data(), cur.Data());
 209
 210     TDirectory* dCut = out->mkdir(cuts1->GetName());
 211     TIter       iCut(cuts1);
 212     TObject*    pCut = 0;
 213     dCut->cd();
 214     while ((pCut = iCut())) {
 215       TString     sMethod = pCut->GetName();
 216       TDirectory* dMethod = dCut->mkdir(sMethod);
 217       TString     sValues = pCut->GetTitle();
 218       TObjArray*  aValues = sValues.Tokenize(" ");
 219       TIter       iValue(aValues);
 220       TObjString* pValue = 0;
 221       TString     templ(Form("%s%s", cuts1->GetName(), sMethod.Data()));
 222       TString     base(cur);
 223       base.ReplaceAll(Form("%sX", cuts1->GetName()), templ);
 224       aValues->SetName(cuts1->GetName());
 225
 226       THStack*     mean   = 0;
 227       THStack*     ratios = 0;
 228       dMethod->cd();
 229       MakeStacks(base, aValues, mean, ratios);
 230       dMethod->Add(mean);
 231       dMethod->Add(ratios);
 232
 233       Int_t xbin = 1;
 234       while ((pValue = static_cast<TObjString*>(iValue()))) {
 235         TString     sValue(Form("%05.2f", pValue->String().Atof()));
 236         sValue.ReplaceAll(".", "d");
 237         TDirectory* dValue = dMethod->mkdir(sValue);
 238         TString now(base);
 239         TString sub(Form("%s%s", templ.Data(), sValue.Data()));
 240         now.ReplaceAll(templ.Data(), sub.Data());
 241         // now.Append(sValue);
 242         dValue->cd();
 243
 244         if (cuts2) {
 245           // Loop over sub-cut.
 246           LoopCuts(run, now, cuts2, cuts3, 0,
 247                    dValue, xbin, mean, ratios);
 248         }
 249         else {
 250           // Process for a given cut with the current value of that cut.
 251           if (!NextFile(run, now, dValue, xbin,
 252                         mean, ratios)) {
 253             // xbin++;
 254             dMethod->cd();
 255             continue;
 256           }
 257           // After this, we have points for the current cut value
 258           // stored in the stacks mean and wSpread, and added ratios
 259           // for the current cut to the stack ratios.  Since we're not
 260           // done with the cut values just yet, we should wait to
 261           // update the parent stacks.
 262         }
 263
 264         if (upMean /*&& upWSpread*/) {
 265           // For each centrality extract
 266           // - mean of mean of ...
 267           // - var of var of ...
 268           // - spread of spread of ...
 269           TIter    iCent(&fCents);
 270           TObject* pCent = 0;
 271           Int_t    jCent = 0;
 272           while ((pCent = iCent())) {
 273             TH1* h = static_cast<TH1*>(mean->GetHists()->At(jCent));
 274             // Info("", "Updating parent %s with %s @ %d",
 275             //      upMean->GetTitle(), h->GetTitle(),upBin);
 276             UpdateStacks(h, jCent, upBin, upMean);
 277             jCent++;
 278           }
 279         } // if ups
 280
 281         xbin++;
 282         dMethod->cd();
 283       } // for values
 284
 285       if (upRatios) {
 286         TIter nextHist(ratios->GetHists());
 287         TH1*  hist = 0;
 288         while ((hist = static_cast<TH1*>(nextHist())))
 289           upRatios->Add(hist);
 290       }
 291       dCut->cd();
 292
 293       FixMinMax(ratios);
 294       FixMinMax(mean);
 295
 296       if (mean && mean->GetHists()->GetEntries() > 0 &&
 297           ratios && ratios->GetHists() &&
 298           ratios->GetHists()->GetEntries() > 0) {
 299         fBody->Divide(2,1);
 300         DrawStacks(fBody->GetPad(1), mean, kNorth|kWest);
 301         DrawStacks(fBody->GetPad(2), ratios, kNorth|kCenter, kSouth|kCenter);
 302         PrintCanvas(Form("%s_%s", pre.Data(), base.Data()), 0.5);
 303       }
 304     } // for methods
 305     out->cd();
 306   }
 307   //__________________________________________________________________
 308   /**
 309    * Process the next file
 310    *
 311    * @param dir        Directory
 312    * @param filename   File name
 313    * @param out        Output directory
 314    * @param binx       Bin number
 315    * @param mean       Graphs of mean +/- variance
 316    * @param wspread    Graphs of mean +/- max/min
 317    * @param ratios     Ratios stack
 318    *
 319    * @return true on success
 320    */
 321   Bool_t NextFile(const TString& run, const TString& now,
 322                   TDirectory* out, Int_t binx,
 323                   THStack* mean, THStack* ratios)
 324   {
 325     TString dir(Form("%s_dndeta_%s", run.Data(), now.Data()));
 326     TString path(Form("%s/forward_dndeta.root", dir.Data()));
 327     if (gSystem->AccessPathName(path.Data())) {
 328       Warning("NextFile", "%s not found", path.Data());
 329       return false;
 330     }
 331     TFile* file = TFile::Open(path, "READ");
 332     if (!file) {
 333       Warning("NextFile", "Failed to open %s", path.Data());
 334       return false;
 335     }
 336     // Info("NextFile", "Opened %s", path.Data());
 337
 338     TCollection* results = GetCollection(file, "ForwarddNdetaResults");
 339     if (!results) return false;
 340
 341     THStack* all = new THStack("all",    "All");
 342     THStack* rat = new THStack("ratios", "Ratios");
 343
 344     TIter    iCent(&fCents);
 345     TObject* pCent = 0;
 346     Int_t    jCent = 0;
 347     while ((pCent = iCent())) {
 348       TGraph*  graph = static_cast<TGraph*>(fOthers.At(jCent));
 349       if (!graph) break;
 350
 351       TString folderName(pCent->GetName());
 352       TCollection* centFolder = GetCollection(results, folderName);
 353       if (!centFolder) break;
 354
 355       TH1* dNdeta = GetH1(centFolder, Form("dndetaForward%s",
 356                                            fRebinned ? "_rebin05" : ""));
 357       if (!dNdeta) {
 358         Warning("", "Didn't get histogram for jCent=%d in %s",
 359                 jCent, path.Data());
 360         break;
 361       }
 362       dNdeta->SetDirectory(out);
 363       dNdeta->SetName(folderName);
 364       dNdeta->SetMarkerColor(jCent+1);
 365       dNdeta->SetLineColor(jCent+1);
 366
 367       TH1* other = G2H(graph, *(dNdeta->GetXaxis()));
 368       other->SetMarkerColor(dNdeta->GetMarkerColor());
 369       other->SetMarkerSize(dNdeta->GetMarkerSize());
 370       other->SetLineColor(dNdeta->GetLineColor());
 371       other->SetDirectory(out);
 372       folderName.ReplaceAll("cent", "other");
 373       folderName.ReplaceAll("all",  "other");
 374       other->SetName(folderName);
 375
 376       all->Add(other);
 377       all->Add(dNdeta);
 378
 379       folderName.ReplaceAll("other", "ratio");
 380       TH1* ratio = static_cast<TH1*>(dNdeta->Clone(folderName));
 381       ratio->SetTitle(Form("%s %s", dir.Data(), pCent->GetName()));
 382       ratio->SetDirectory(out);
 383       ratio->Divide(other);
 384       ratio->SetMarkerStyle(20+binx-1);
 385       ratio->SetMarkerColor(jCent+1);
 386       ratio->SetLineColor(kGray);
 387       ratio->SetLineWidth(0);
 388       ratios->Add(ratio);
 389
 390       rat->Add(ratio);
 391
 392       UpdateStacks(ratio, jCent, binx, mean);
 393       jCent++;
 394     }
 395     out->Add(all);
 396     out->Add(rat);
 397
 398     FixMinMax(all);
 399     FixMinMax(rat);
 400
 401     fBody->Divide(1,2,0,0);
 402     DrawStacks(fBody->GetPad(1), all/*, 21*/);
 403     DrawStacks(fBody->GetPad(2), rat, kNorth|kCenter);
 404     PrintCanvas(Form("____%s", now.Data()), .4);
 405
 406     file->Close();
 407     if (jCent <= 0) return false;
 408
 409     return true;
 410   }
 411   //=== Stack functions ==============================================
 412   /**
 413    * Make stacks
 414    *
 415    * @param run      Run number
 416    * @param values   Cut parameters
 417    * @param mean     Stack of means
 418    * @param ratios   Stack of ratios
 419    */
 420   void MakeStacks(const TString&   run,
 421                   const TObjArray* values,
 422                   THStack*&        mean,
 423                   THStack*&        ratios)
 424   {
 425     mean    = new THStack("mean", run);
 426     ratios  = new THStack("ratios", run);
 427
 428     // --- Create histograms and graphs ------------------------
 429     Int_t    nValues = values->GetEntriesFast();
 430     TIter    nextCent(&fCents);
 431     TObject* pcent = 0;
 432     Int_t    col   = 1;
 433     while ((pcent = nextCent())) {
 434       UShort_t low   = (pcent->GetUniqueID() & 0xFFFF);
 435       UShort_t high  = (pcent->GetUniqueID() >> 16) & 0xFFFF;
 436
 437       TH1* hMean = new TH1D(pcent->GetName(),
 438                             Form("%s %d%%-%d%% central", run.Data(),
 439                                  low, high),
 440                             nValues, .5, nValues+.5);
 441       hMean->SetMarkerColor(col);
 442       hMean->SetMarkerStyle(20);
 443       hMean->SetLineColor(col);
 444       hMean->SetXTitle(Form("Cut parameter X_{%s}", values->GetName()));
 445       hMean->SetYTitle("Average, spread, and min/max of ratio");
 446       hMean->SetDirectory(0);
 447       TIter nextV(values);
 448       TObjString* pvalue = 0;
 449       Int_t xbin = 1;
 450       while ((pvalue = static_cast<TObjString*>(nextV()))) {
 451         TString& value = pvalue->String();
 452         hMean->GetXaxis()->SetBinLabel(xbin, value);
 453         xbin++;
 454       }
 455
 456       TGraphAsymmErrors* gWSpread = new TGraphAsymmErrors(nValues);
 457       gWSpread->SetName("minmax");
 458       gWSpread->SetTitle(hMean->GetTitle());
 459       gWSpread->SetMarkerColor(col);
 460       gWSpread->SetLineColor(col);
 461       gWSpread->SetFillColor(0);
 462       gWSpread->SetFillStyle(0);
 463
 464       hMean->GetListOfFunctions()->Add(gWSpread, "[]pl same");
 465       mean->Add(hMean, "x0 e1");
 466       col++;
 467     }
 468   }
 469   //__________________________________________________________________
 470   /**
 471    * Update stacks
 472    *
 473    * @param h       Histogram
 474    * @param i       Index
 475    * @param binx    Bin number
 476    * @param stack   Stack to update
 477    */
 478   void UpdateStacks(TH1*         h,
 479                     Int_t        i,
 480                     Int_t        binx,
 481                     THStack*     stack)
 482   {
 483     Double_t avg = 0;
 484     Double_t var = 0;
 485     Double_t min = +100000;
 486     Double_t max = -100000;
 487     HistStatistics(h, avg, var, min, max);
 488     h->SetMinimum(0.95*min);
 489     h->SetMaximum(1.1*max);
 490
 491     TH1* hMean = static_cast<TH1*>(stack->GetHists()->At(i));
 492     hMean->SetBinContent(binx, avg);
 493     hMean->SetBinError(binx,var);
 494
 495     TObject* pG = hMean->GetListOfFunctions()->FindObject("minmax");
 496     if (!pG) {
 497       hMean->GetListOfFunctions()->ls();
 498       return;
 499     }
 500     TGraphAsymmErrors* gWSpread = static_cast<TGraphAsymmErrors*>(pG);
 501     gWSpread->SetPoint(binx-1, binx, avg);
 502     gWSpread->SetPointError(binx-1,0, 0, /*0.5,0.5,*/
 503                             TMath::Abs(avg-min), TMath::Abs(max-avg));
 504   }
 505   //=== Graphics member functions ====================================
 506   //__________________________________________________________________
 507   /**
 508    * Build a centrality legend
 509    *
 510    * @param p     Pad
 511    * @param where See above
 512    * @param stack Stack to make legend for
 513    */
 514   void BuildCentLegend(TVirtualPad* p, UInt_t where, THStack* stack=0)
 515   {
 516     TLegend* l     = MakeLegend(p, where, false);
 517     TIter    iCent(&fCents);
 518     TObject* pCent = 0;
 519     Int_t    col   = 1;
 520     while ((pCent = iCent())) {
 521       UShort_t low   = (pCent->GetUniqueID() & 0xFFFF);
 522       UShort_t high  = (pCent->GetUniqueID() >> 16) & 0xFFFF;
 523
 524       TLegendEntry* e = l->AddEntry("dummy",
 525                                     Form("%2d%% - %2d%% central", low, high),
 526                                     "pl");
 527       e->SetFillColor(col);
 528       e->SetMarkerColor(col);
 529       e->SetLineColor(col);
 530       e->SetMarkerStyle(20);
 531       col++;
 532     }
 533     if (stack && stack->GetHistogram()) {
 534       stack->GetHistogram()->GetListOfFunctions()->Add(l);
 535     }
 536     else
 537       l->Draw();
 538
 539   }
 540   //__________________________________________________________________
 541   /**
 542    * Make a cut legend.
 543    *
 544    * @param p      Pad
 545    * @param where  See above
 546    * @param stack  Stack to make legend for
 547    */
 548   void BuildCutLegend(TVirtualPad* p, UInt_t where, THStack* stack)
 549   {
 550     TLegend* l     = MakeLegend(p, where, false);
 551     l->SetX1(p->GetLeftMargin());
 552     l->SetX2(1-p->GetRightMargin());
 553     // l->SetBorderSize(1);
 554     l->SetNColumns(1);
 555
 556     TList    seen;
 557     TIter    iHist(stack->GetHists());
 558     TH1*     pHist = 0;
 559     while ((pHist = static_cast<TH1*>(iHist()))) {
 560       TString title(pHist->GetTitle());
 561
 562       Int_t   idx = title.Index(" cent");
 563       if (idx != kNPOS) title.Remove(idx, 5+3+3+1);
 564
 565       idx = title.Index("_dndeta");
 566       if (idx != kNPOS) title.Remove(0, idx+6+1+1);
 567
 568       TObject* before = seen.FindObject(title);
 569       if (before) continue;
 570
 571       seen.Add(new TObjString(title));
 572
 573       TLegendEntry* e = l->AddEntry("dummy", title, "p");
 574       e->SetMarkerColor(kBlack);
 575       e->SetMarkerStyle(pHist->GetMarkerStyle());
 576     }
 577     seen.IsOwner();
 578     if (stack->GetHistogram()) {
 579       stack->GetHistogram()->GetListOfFunctions()->Add(l);
 580     }
 581     else
 582       l->Draw();
 583
 584   }
 585   //__________________________________________________________________
 586   /**
 587    * Draw stacks
 588    *
 589    * @param p      Pad
 590    * @param stack  Stacks
 591    * @param cent   Centrality legend location
 592    * @param cuts   Cut legend location
 593    */
 594   void DrawStacks(TVirtualPad* p,
 595                   THStack* stack,
 596                   UInt_t cent=0,
 597                   UInt_t cuts=0)
 598   {
 599     if (!stack) {
 600       Warning("DrawStacks", "Stack is missing!");
 601       return;
 602     }
 603     TH1*    first = static_cast<TH1*>(stack->GetHists()->At(0));
 604     TString xT    = first->GetXaxis()->GetTitle();
 605
 606     p->cd();
 607     stack->Draw("nostack");
 608     // DrawInPad(p, 0, stack, "nostack", flags);
 609     stack->GetXaxis()->SetTitle(xT);
 610     FixMinMax(stack);
 611
 612     if (cent > 0) BuildCentLegend(p, cent, stack);
 613     if (cuts > 0) BuildCutLegend(p, cuts, stack);
 614
 615     p->Modified();
 616     p->Update();
 617     p->cd();
 618   }
 619   //=== Utility static member functions ==============================
 620   /**
 621    * Update statistics
 622    *
 623    * @param y    Current observation
 624    * @param cnt  On entry, last count, on return current count
 625    * @param mean On entry, last mean, on return current mean
 626    * @param var  On entry, last variance, on return current variance
 627    */
 628   static void Statistics(Double_t  y,
 629                          Int_t&    cnt,
 630                          Double_t& mean,
 631                          Double_t& var)
 632   {
 633     cnt        += 1;
 634     mean       += (y - mean) / cnt;
 635     var        += (cnt > 1 ? (TMath::Power(y-mean,2)/(cnt-1)-var/cnt) : 0);
 636   }
 637   //__________________________________________________________________
 638   /**
 639    * Calculate histogram statistics
 640    *
 641    * @param h     Histogram
 642    * @param mean  On return, the Y-mean
 643    * @param var   On return, the Y variance
 644    * @param min   On return, the least Y
 645    * @param max   On return, the largest Y
 646    */
 647   static void HistStatistics(const TH1* h,
 648                              Double_t& mean,
 649                              Double_t& var,
 650                              Double_t& min,
 651                              Double_t& max)
 652   {
 653     mean      = 0;
 654     var       = 0;
 655     min       = +100000;
 656     max       = -100000;
 657     Int_t cnt = 0;
 658     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 659       Double_t y = h->GetBinContent(i);
 660       if (TMath::Abs(y) <= 1e-9) continue;
 661       min        =  TMath::Min(min, y);
 662       max        =  TMath::Max(max, y);
 663       Statistics(y, cnt, mean, var);
 664     }
 665     // Info("", "Stats for %s:  mean=%f +/- %f [%f,%f]",
 666     //      h->GetTitle(), mean, var, min, max);
 667   }
 668   //__________________________________________________________________
 669   /**
 670    * Fix the min/max of a stack
 671    *
 672    * @param stack
 673    */
 674   static void FixMinMax(THStack* stack)
 675   {
 676     TIter iHist(stack->GetHists());
 677     TH1*  pHist = 0;
 678     Double_t m1 = 10000000;
 679     Double_t m2 = -10000000;
 680     while((pHist = static_cast<TH1*>(iHist()))) {
 681       m1 = TMath::Min(m1, pHist->GetMinimum(1e-6));
 682       m2 = TMath::Max(m2, pHist->GetMaximum());
 683     }
 684     // Double_t m1 = stack->GetMinimum("nostack e");
 685     // Double_t m2 = stack->GetMaximum("nostack e");
 686     // Printf("Stack %s minimum: %f", stack->GetTitle(), m1);
 687     stack->SetMinimum((m1 < 0 ? 1.05 : 0.95)  * m1);
 688     stack->SetMaximum((m2 > 0 ? 1.05 : 0.95)  * m2);
 689   }
 690   //__________________________________________________________________
 691   /**
 692    * Turn a graph into a histogram
 693    *
 694    * @param g    Graph
 695    * @param axis Axis to use
 696    *
 697    * @return Newly allocated histogram
 698    */
 699   static TH1* G2H(const TGraph* g, const TAxis& axis)
 700   {
 701     TH1* h = 0;
 702     if (axis.GetXbins()->GetArray())
 703       h = new TH1D(g->GetName(), g->GetTitle(),
 704                    axis.GetNbins(), axis.GetXbins()->GetArray());
 705     else
 706       h = new TH1D(g->GetName(), g->GetTitle(),
 707                    axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
 708     h->SetMarkerColor(g->GetMarkerColor());
 709     h->SetMarkerStyle(g->GetMarkerStyle());
 710     h->SetMarkerSize(g->GetMarkerSize());
 711     h->SetLineColor(g->GetLineColor());
 712     h->SetLineStyle(g->GetLineStyle());
 713     h->SetLineWidth(g->GetLineWidth());
 714     h->SetFillColor(g->GetFillColor());
 715     h->SetFillStyle(g->GetFillStyle());
 716     h->SetDirectory(0);
 717
 718     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 719       Double_t x = h->GetXaxis()->GetBinCenter(i);
 720       Double_t y = g->Eval(x); // , 0, "S");
 721       h->SetBinContent(i, y);
 722     }
 723     return h;
 724   }
 725   //__________________________________________________________________
 726   /**
 727    * Get other data
 728    *
 729    * @param sys    Collision system
 730    * @param sNN    Collision energy
 731    * @param trg    Trigger
 732    * @param lowC   Low centrality
 733    * @param highC  High centrality
 734    *
 735    * @return Newly allocated histogram
 736    */
 737   static TGraph* GetOther(UShort_t      sys,
 738                           UShort_t      sNN,
 739                           UShort_t      trg,
 740                           UShort_t      lowC,
 741                           UShort_t      highC)
 742   {
 743     // --- Set the macro pathand load other data script --------------
 744     // Always recompile
 745     if (!gROOT->GetClass("RefData")) {
 746       TString savPath(gROOT->GetMacroPath());
 747       gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2",
 748                                gROOT->GetMacroPath()));
 749       gROOT->LoadMacro("OtherData.C++");
 750       gROOT->SetMacroPath(savPath);
 751     }
 752     Long_t ret = gROOT->ProcessLine(Form("RefData::GetData(%d,%d,%d,%d,%d,0x4)",
 753                                          sys, sNN, trg, lowC, highC));
 754     if (!ret) {
 755       Warning("GetOther", "No other data for %d %d %d %d%%-%d%% central",
 756               sys, sNN, trg, lowC, highC);
 757       return 0;
 758     }
 759     TMultiGraph* others = reinterpret_cast<TMultiGraph*>(ret);
 760     TGraph*      other  =static_cast<TGraph*>(others->GetListOfGraphs()->At(0));
 761     if (!other) {
 762       Warning("GetOther", "No ALICE data for %d %d %d %d%%-%d%% central",
 763               sys, sNN, trg, lowC, highC);
 764       return 0;
 765     }
 766
 767     // Info("", "Got graph %s/%s", other->GetName(), other->GetTitle());
 768
 769     return other;
 770   }
 771
 772   //__________________________________________________________________
 773   TList     fSLCuts;
 774   TList     fSHCuts;
 775   TList     fDCCuts;
 776   TList     fRuns;
 777   TList     fCents;
 778   TList     fOthers;
 779   TList*    fOrder[3];
 780   Bool_t    fRebinned;
 781 };
 782
 783 #endif