PWG2/FORWARD/analysis2/DrawdNdeta.C

   1 /**
   2  * Script to visualise the dN/deta
   3  *
   4  * This script is independent of any AliROOT code - and should stay that way.
   5  */
   6 #include <TH1.h>
   7 #include <THStack.h>
   8 #include <TGraphErrors.h>
   9 #include <TGraphAsymmErrors.h>
  10 #include <TMultiGraph.h>
  11 #include <TFile.h>
  12 #include <TList.h>
  13 #include <TString.h>
  14 #include <TError.h>
  15 #include <TSystem.h>
  16 #include <TROOT.h>
  17 #include <TMath.h>
  18 #include <TCanvas.h>
  19 #include <TPad.h>
  20 #include <TStyle.h>
  21 #include <TLegend.h>
  22 #include <TLegendEntry.h>
  23 #include <TLatex.h>
  24 #include <TImage.h>
  25
  26 /**
  27  * Class to draw dN/deta results
  28  *
  29  */
  30 struct dNdetaDrawer
  31 {
  32   /**
  33    * POD of data for range zooming
  34    */
  35   struct RangeParam
  36   {
  37     TAxis*       fMasterAxis; // Master axis
  38     TAxis*       fSlave1Axis; // First slave axis
  39     TVirtualPad* fSlave1Pad;  // First slave pad
  40     TAxis*       fSlave2Axis; // Second slave axis
  41     TVirtualPad* fSlave2Pad;  // Second slave pad
  42   };
  43   //__________________________________________________________________
  44   /**
  45    * Constructor
  46    *
  47    */
  48   dNdetaDrawer()
  49     : fShowOthers(false),       // Bool_t
  50       fShowAlice(false),        // Bool_t
  51       fShowRatios(false),       // Bool_t
  52       fShowLeftRight(false),    // Bool_t
  53       fRebin(5),                // UShort_t
  54       fCutEdges(false),         // Bool_t
  55       fTitle(""),               // TString
  56       fTrigString(0),           // TNamed*
  57       fSNNString(0),            // TNamed*
  58       fSysString(0),            // TNamed*
  59       fVtxAxis(0),              // TAxis*
  60       // fForward(0),           // TH1*
  61       // fForwardMC(0),         // TH1*
  62       // fForwardHHD(0),        // TH1*
  63       // fTruth(0),             // TH1*
  64       // fCentral(0),           // TH1*
  65       // fForwardSym(0),        // TH1*
  66       // fForwardMCSym(0),      // TH1*
  67       // fForwardHHDSym(0),     // TH1*
  68       fTriggers(0),             // TH1*
  69       fRangeParam(0)
  70
  71   {
  72     fRangeParam = new RangeParam;
  73     fRangeParam->fMasterAxis = 0;
  74     fRangeParam->fSlave1Axis = 0;
  75     fRangeParam->fSlave1Pad  = 0;
  76     fRangeParam->fSlave2Axis = 0;
  77     fRangeParam->fSlave2Pad  = 0;
  78   }
  79   //==================================================================
  80   /**
  81    * @{
  82    * @name Set parameters
  83    */
  84   /**
  85    * Show other (UA5, CMS, ...) data
  86    *
  87    * @param x Whether to show or not
  88    */
  89   void SetShowOthers(Bool_t x)    { fShowOthers = x; }
  90   //__________________________________________________________________
  91   /**
  92    * Show ALICE published data
  93    *
  94    * @param x Wheter to show or not
  95    */
  96   void SetShowAlice(Bool_t x)     { fShowAlice = x; }
  97   //__________________________________________________________________
  98   /**
  99    * Whether to show ratios or not.  If there's nothing to compare to,
 100    * the ratio panel will be implicitly disabled
 101    *
 102    * @param x Whether to show or not
 103    */
 104   void SetShowRatios(Bool_t x)    { fShowRatios = x; }
 105   //__________________________________________________________________
 106   /**
 107    *
 108    * Whether to show the left/right asymmetry
 109    *
 110    * @param x To show or not
 111    */
 112   void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
 113   //__________________________________________________________________
 114   /**
 115    * Set the rebinning factor
 116    *
 117    * @param x Rebinning factor (must be a divisor in the number of bins)
 118    */
 119   void SetRebin(UShort_t x)       { fRebin = x; }
 120   //__________________________________________________________________
 121   /**
 122    * Wheter to cut away the edges
 123    *
 124    * @param x Whether or not to cut away edges
 125    */
 126   void SetCutEdges(Bool_t x)      { fCutEdges = x; }
 127   //__________________________________________________________________
 128   /**
 129    * Set the title of the plot
 130    *
 131    * @param x Title
 132    */
 133   void SetTitle(TString x)        { fTitle = x; }
 134   /* @} */
 135   //==================================================================
 136   /**
 137    * @{
 138    * @name Override settings from input
 139    */
 140   /**
 141    * Override setting from file
 142    *
 143    * @param sNN Center of mass energy per nucleon pair (GeV)
 144    */
 145   void SetSNN(UShort_t sNN)
 146   {
 147     fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
 148     fSNNString->SetUniqueID(sNN);
 149   }
 150   //__________________________________________________________________
 151   /**
 152    * Set the collision system
 153    * - 1: pp
 154    * - 2: PbPb
 155    *
 156    * @param sys collision system
 157    */
 158   void SetSys(UShort_t sys)
 159   {
 160     fSysString = new TNamed("sys", (sys == 1 ? "pp" :
 161                                     sys == 2 ? "PbPb" : "unknown"));
 162     fSysString->SetUniqueID(sys);
 163   }
 164   //__________________________________________________________________
 165   /**
 166    * Set the vertex range in centimeters
 167    *
 168    * @param vzMin Min @f$ v_z@f$
 169    * @param vzMax Max @f$ v_z@f$
 170    */
 171   void SetVertexRange(Double_t vzMin, Double_t vzMax)
 172   {
 173     fVtxAxis = new TAxis(10, vzMin, vzMax);
 174     fVtxAxis->SetName("vtxAxis");
 175     fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax));
 176   }
 177   //__________________________________________________________________
 178   void SetTrigger(UShort_t trig)
 179   {
 180     fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" :
 181                                             trig & 0x2 ? "INEL>0" :
 182                                             trig & 0x4 ? "NSD" :
 183                                             "unknown"));
 184     fTrigString->SetUniqueID(trig);
 185   }
 186
 187
 188   //==================================================================
 189   /**
 190    * @{
 191    * @name Main steering functions
 192    */
 193   /**
 194    * Run the code to produce the final result.
 195    *
 196    * @param filename  File containing the data
 197    */
 198   void Run(const char* filename="forward_dndeta.C")
 199   {
 200     Double_t max = 0;
 201
 202     if (!Open(filename, max)) return;
 203     Info("Run", "Got data from %s with maximum %f", filename, max);
 204
 205     // Create our stack of results
 206     THStack* results = fResults; // StackResults(max);
 207
 208     // Create our stack of other results
 209     TMultiGraph* other = 0;
 210     if (fShowOthers || fShowAlice) other = StackOther(max);
 211
 212     Double_t smax = 0;
 213     THStack* ratios = 0;
 214     if (fShowRatios) ratios = StackRatios(other, smax);
 215
 216     Double_t amax = 0;
 217     THStack* leftright = 0;
 218     if (fShowLeftRight) leftright = StackLeftRight(amax);
 219
 220     Plot(results, other, max, ratios, smax, leftright, amax);
 221   }
 222
 223   //__________________________________________________________________
 224   /**
 225    * Fetch the information on the run from the results list
 226    *
 227    * @param results  Results list
 228    */
 229   void FetchInformation(const TList* results)
 230   {
 231     if (!fTrigString)
 232       fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
 233     if (!fSNNString)
 234       fSNNString  = static_cast<TNamed*>(results->FindObject("sNN"));
 235     if (!fSysString)
 236       fSysString  = static_cast<TNamed*>(results->FindObject("sys"));
 237     if (!fVtxAxis)
 238       fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
 239     if (!fTrigString) fTrigString = new TNamed("trigString", "unknown");
 240     if (!fSNNString)  fSNNString  = new TNamed("sNN", "unknown");
 241     if (!fSysString)  fSysString  = new TNamed("sys", "unknown");
 242     if (!fVtxAxis) {
 243       fVtxAxis    = new TAxis(1,0,0);
 244       fVtxAxis->SetName("vtxAxis");
 245       fVtxAxis->SetTitle("v_{z} range unspecified");
 246     }
 247
 248     Info("FetchInformation",
 249          "Initialized for\n"
 250          "   Trigger:    %s  (%d)\n"
 251          "   sqrt(sNN):  %s  (%dGeV)\n"
 252          "   System:     %s  (%d)\n"
 253          "   Vz range:   %s  (%f,%f)",
 254          fTrigString->GetTitle(), fTrigString->GetUniqueID(),
 255          fSNNString->GetTitle(),  fSNNString->GetUniqueID(),
 256          fSysString->GetTitle(),  fSysString->GetUniqueID(),
 257          fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
 258   }
 259
 260   //__________________________________________________________________
 261   /**
 262    * Open input file, and find data
 263    *
 264    * @param filename File name
 265    *
 266    * @return true on success
 267    */
 268   Bool_t Open(const char* filename, Double_t& max)
 269   {
 270     TFile* file = TFile::Open(filename, "READ");
 271     if (!file) {
 272       Error("Open", "Cannot open %s", filename);
 273       return false;
 274     }
 275
 276     TList* results = static_cast<TList*>(file->Get("ForwardResults"));
 277     if (!results) {
 278       Error("Open", "Couldn't find list ForwardResults");
 279       return false;
 280     }
 281
 282     // Get information on the run
 283     FetchInformation(results);
 284
 285     TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
 286     if (!clusters) Warning("Open", "Couldn't find list CentralResults");
 287
 288     fResults   = new THStack("results", "Results");
 289     FetchResults(fResults, results,  "Forward", max);
 290     FetchResults(fResults, clusters, "Central", max);
 291
 292     TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
 293     if (sums) fTriggers = FetchResult(sums, "triggers");
 294
 295     file->Close();
 296
 297     return true;
 298   }
 299   //__________________________________________________________________
 300   void FetchResults(THStack* stack, const TList* list, const char* name,
 301                     Double_t& max)
 302   {
 303     TH1* dndeta      = FetchResult(list, Form("dndeta%s", name));
 304     TH1* dndetaMC    = FetchResult(list, Form("dndeta%sMC", name));
 305     TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
 306     TH1* dndetaSym   = 0;
 307     TH1* dndetaMCSym = 0;
 308
 309     max = TMath::Max(max, AddHistogram(stack, dndetaTruth, "e5 p"));
 310     max = TMath::Max(max, AddHistogram(stack, dndetaMC,    "", dndetaMCSym));
 311     max = TMath::Max(max, AddHistogram(stack, dndeta,      "", dndetaSym));
 312
 313     Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
 314          dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
 315
 316     if (dndetaMCSym) fNumerators.Add(dndetaMCSym);
 317     if (dndetaMC)    fNumerators.Add(dndetaMC);
 318     if (dndetaSym)   fNumerators.Add(dndetaSym);
 319     if (dndeta)      fNumerators.Add(dndeta);
 320     if (dndetaTruth) fDenominators.Add(dndetaTruth);
 321   }
 322   //__________________________________________________________________
 323   /**
 324    * Make a histogram stack of results
 325    *
 326    * @param max On return, the maximum value in the stack
 327    *
 328    * @return Newly allocated stack
 329    */
 330   TMultiGraph* StackOther(Double_t& max)
 331   {
 332     gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C");
 333     Int_t    error = 0;
 334     Bool_t   onlya = (fShowOthers ? false : true);
 335     Int_t    trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
 336     UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
 337     Long_t   ret   = gROOT->ProcessLine(Form("GetData(%d,%d,%d);",
 338                                              snn,trg,onlya));
 339     if (error) {
 340       Error("StackOther", "Failed to execute GetData(%d,%d,%d)",
 341             snn, trg, onlya);
 342       return 0;
 343     }
 344     if (!ret) {
 345       Warning("StackOther", "No other data found for sNN=%d, trigger=%d",
 346               snn, trg);
 347       return 0;
 348     }
 349     TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret);
 350
 351     TGraphAsymmErrors* o      = 0;
 352     TIter              next(other->GetListOfGraphs());
 353     while ((o = static_cast<TGraphAsymmErrors*>(next()))) {
 354       max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY()));
 355       fDenominators.Add(o);
 356     }
 357
 358     return other;
 359   }
 360   //__________________________________________________________________
 361   /**
 362    * Make a histogram stack of ratios of results to other data
 363    *
 364    * @param max On return, the maximum value in the stack
 365    *
 366    * @return Newly allocated stack
 367    */
 368   THStack* StackRatios(TMultiGraph* others, Double_t& max)
 369   {
 370     THStack* ratios = new THStack("ratios", "Ratios");
 371
 372     if (others) {
 373       TObject* ua5_1  = 0;
 374       TObject* ua5_2  = 0;
 375       TObject* alice  = 0;
 376       TObject* cms    = 0;
 377       TObject* denom  = 0;
 378       TIter    nextD(&fDenominators);
 379       while ((denom = nextD())) {
 380         TIter nextN(&fNumerators);
 381         TObject* numer = 0;
 382         while ((numer = nextN())) {
 383           ratios->Add(Ratio(numer, denom, max));
 384         }
 385         TString oName(denom->GetName());
 386         oName.ToLower();
 387         if (oName.Contains("ua5"))  {
 388           if (ua5_1) ua5_2 = denom;
 389           else       ua5_1 = denom;
 390         }
 391         if (oName.Contains("alice")) alice = denom;
 392         if (oName.Contains("cms"))   cms   = denom;
 393       }
 394       if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max));
 395       if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max));
 396       if (cms   && alice) ratios->Add(Ratio(alice, cms,   max));
 397     }
 398     // Check if we have ratios
 399     if (!ratios->GetHists() ||
 400         (ratios->GetHists()->GetEntries() <= 0)) {
 401       delete ratios;
 402       ratios = 0;
 403     }
 404     return ratios;
 405   }
 406   //__________________________________________________________________
 407   /**
 408    * Make a histogram stack of the left-right asymmetry
 409    *
 410    * @param max On return, the maximum value in the stack
 411    *
 412    * @return Newly allocated stack
 413    */
 414   THStack* StackLeftRight(Double_t& max)
 415   {
 416     THStack* ret = new THStack("leftright", "Left-right asymmetry");
 417     // ret->Add(Asymmetry(fForward,    max));
 418     // ret->Add(Asymmetry(fForwardMC,  max));
 419
 420     if (!ret->GetHists() ||
 421         (ret->GetHists()->GetEntries() <= 0)) {
 422       delete ret;
 423       ret = 0;
 424     }
 425     return ret;
 426   }
 427   //__________________________________________________________________
 428   /**
 429    * Plot the results
 430    *
 431    * @param results    Results
 432    * @param others     Other data
 433    * @param max        Max value
 434    * @param ratios     Stack of ratios (optional)
 435    * @param rmax       Maximum diviation from 1 of ratios
 436    * @param leftright  Stack of left-right asymmetry (optional)
 437    * @param amax       Maximum diviation from 1 of asymmetries
 438    */
 439   void Plot(THStack*     results,
 440             TMultiGraph* others,
 441             Double_t     max,
 442             THStack*     ratios,
 443             Double_t     rmax,
 444             THStack*     leftright,
 445             Double_t     amax)
 446   {
 447     gStyle->SetOptTitle(0);
 448     gStyle->SetTitleFont(132, "xyz");
 449     gStyle->SetLabelFont(132, "xyz");
 450
 451     Int_t    h = 800;
 452     Int_t    w = 800; // h / TMath::Sqrt(2);
 453     if (!ratios) w *= 1.4;
 454     if (!leftright) w *= 1.4;
 455     TCanvas* c = new TCanvas("Results", "Results", w, h);
 456     c->SetFillColor(0);
 457     c->SetBorderSize(0);
 458     c->SetBorderMode(0);
 459
 460     Double_t y1 = 0;
 461     Double_t y2 = 0;
 462     Double_t y3 = 0;
 463     if (ratios)    y1 = 0.3;
 464     if (leftright) {
 465       if (y1 > 0.0001) {
 466         y2 = 0.2;
 467         y1 = 0.4;
 468       }
 469       else {
 470         y1 = 0.2;
 471         y2 = y1;
 472       }
 473     }
 474     PlotResults(results, others, max, y1);
 475     c->cd();
 476
 477     PlotRatios(ratios, rmax, y2, y1);
 478     c->cd( );
 479
 480     PlotLeftRight(leftright, amax, y3, y2);
 481     c->cd();
 482
 483
 484     Int_t   vMin = fVtxAxis->GetXmin();
 485     Int_t   vMax = fVtxAxis->GetXmax();
 486     TString trg(fTrigString->GetTitle());
 487     Int_t   nev  = fTriggers->GetBinContent(fTriggers->GetNbinsX());
 488     trg          = trg.Strip(TString::kBoth);
 489     TString base(Form("dndeta_%s_%s_%s_%c%02d%c%02dcm_%09dev",
 490                       fSysString->GetTitle(),
 491                       fSNNString->GetTitle(),
 492                       trg.Data(),
 493                       vMin < 0 ? 'm' : 'p',  TMath::Abs(vMin),
 494                       vMax < 0 ? 'm' : 'p',  TMath::Abs(vMax),
 495                       nev));
 496     c->SaveAs(Form("%s.png",  base.Data()));
 497     c->SaveAs(Form("%s.root", base.Data()));
 498     c->SaveAs(Form("%s.C",    base.Data()));
 499   }
 500   //__________________________________________________________________
 501   /**
 502    * Plot the results
 503    *
 504    * @param results   Results
 505    * @param others    Other data
 506    * @param max       Maximum
 507    * @param yd        Bottom position of pad
 508    */
 509   void PlotResults(THStack* results, TMultiGraph* others,
 510                    Double_t max, Double_t yd)
 511   {
 512     // Make a sub-pad for the result itself
 513     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
 514     p1->SetTopMargin(0.05);
 515     p1->SetBorderSize(0);
 516     p1->SetBorderMode(0);
 517     p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
 518     p1->SetRightMargin(0.05);
 519     p1->SetGridx();
 520     p1->SetTicks(1,1);
 521     p1->SetNumber(1);
 522     p1->Draw();
 523     p1->cd();
 524
 525     Info("PlotResults", "Plotting results with max=%f", max);
 526     results->ls();
 527     results->SetMaximum(1.15*max);
 528     results->SetMinimum(yd > 0.00001 ? -0.1 : 0);
 529
 530     FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
 531
 532     p1->Clear();
 533     results->DrawClone("nostack e1");
 534
 535     fRangeParam->fSlave1Axis = results->GetXaxis();
 536     fRangeParam->fSlave1Pad  = p1;
 537
 538     // Draw other data
 539     if (others) {
 540       TGraphAsymmErrors* o      = 0;
 541       TIter              nextG(others->GetListOfGraphs());
 542       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
 543         o->DrawClone("same p");
 544     }
 545
 546     // Make a legend in the main result pad
 547     TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
 548     l->SetNColumns(2);
 549     l->SetFillColor(0);
 550     l->SetFillStyle(0);
 551     l->SetBorderSize(0);
 552     l->SetTextFont(132);
 553
 554     // Put a title on top
 555     TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
 556     tit->SetNDC();
 557     tit->SetTextFont(132);
 558     tit->SetTextSize(0.05);
 559     tit->Draw();
 560
 561     // Put a nice label in the plot
 562     TString     eS;
 563     UShort_t    snn = fSNNString->GetUniqueID();
 564     const char* sys = fSysString->GetTitle();
 565     if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
 566     else            eS = Form("%3dGeV", snn);
 567     TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s",
 568                                            sys,
 569                                            eS.Data(),
 570                                            fTrigString->GetTitle()));
 571     tt->SetNDC();
 572     tt->SetTextFont(132);
 573     tt->SetTextAlign(33);
 574     tt->Draw();
 575
 576     // Put number of accepted events on the plot
 577     Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
 578     TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
 579     et->SetNDC();
 580     et->SetTextFont(132);
 581     et->SetTextAlign(33);
 582     et->Draw();
 583
 584     // Put number of accepted events on the plot
 585     if (fVtxAxis) {
 586       TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
 587       vt->SetNDC();
 588       vt->SetTextFont(132);
 589       vt->SetTextAlign(33);
 590       vt->Draw();
 591     }
 592     // results->Draw("nostack e1 same");
 593
 594     fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName());
 595     fRangeParam->fSlave1Pad  = p1;
 596
 597
 598     // Mark the plot as preliminary
 599     TLatex* pt = new TLatex(.12, .93, "Preliminary");
 600     pt->SetNDC();
 601     pt->SetTextFont(22);
 602     pt->SetTextSize(0.07);
 603     pt->SetTextColor(kRed+1);
 604     pt->SetTextAlign(13);
 605     pt->Draw();
 606
 607     if (!gSystem->AccessPathName("ALICE.png")) {
 608       TPad* logo = new TPad("logo", "logo", .12, .65, .25, .85, 0, 0, 0);
 609       logo->SetFillStyle(0);
 610       logo->Draw();
 611       logo->cd();
 612       TImage* i = TImage::Create();
 613       i->ReadImage("ALICE.png");
 614       i->Draw();
 615     }
 616     p1->cd();
 617   }
 618   //__________________________________________________________________
 619   /**
 620    * Plot the ratios
 621    *
 622    * @param ratios  Ratios to plot (if any)
 623    * @param max     Maximum diviation from 1
 624    * @param y1      Lower y coordinate of pad
 625    * @param y2      Upper y coordinate of pad
 626    */
 627   void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2)
 628   {
 629     if (!ratios) return;
 630     bool isBottom = (y1 < 0.0001);
 631     Double_t yd = y2 - y1;
 632     // Make a sub-pad for the result itself
 633     TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
 634     p2->SetTopMargin(0.001);
 635     p2->SetRightMargin(0.05);
 636     p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
 637     p2->SetGridx();
 638     p2->SetTicks(1,1);
 639     p2->SetNumber(2);
 640     p2->Draw();
 641     p2->cd();
 642
 643     // Fix up axis
 644     FixAxis(ratios, 1/yd/1.7, "Ratios", 7);
 645
 646     ratios->SetMaximum(1+TMath::Max(.22,1.05*max));
 647     ratios->SetMinimum(1-TMath::Max(.32,1.05*max));
 648     p2->Clear();
 649     ratios->DrawClone("nostack e1");
 650
 651
 652     // Make a legend
 653     TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
 654                                   isBottom ? .6 : .4);
 655     l2->SetNColumns(2);
 656     l2->SetFillColor(0);
 657     l2->SetFillStyle(0);
 658     l2->SetBorderSize(0);
 659     l2->SetTextFont(132);
 660
 661     // Make a nice band from 0.9 to 1.1
 662     TGraphErrors* band = new TGraphErrors(2);
 663     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
 664     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
 665     band->SetPointError(0, 0, .1);
 666     band->SetPointError(1, 0, .1);
 667     band->SetFillColor(kYellow+2);
 668     band->SetFillStyle(3002);
 669     band->SetLineStyle(2);
 670     band->SetLineWidth(1);
 671     band->Draw("3 same");
 672     band->DrawClone("X L same");
 673
 674     // Replot the ratios on top
 675     ratios->DrawClone("nostack e1 same");
 676
 677     if (isBottom) {
 678       fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName());
 679       p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
 680                                 fRangeParam));
 681     }
 682     else {
 683       fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName());
 684       fRangeParam->fSlave2Pad  = p2;
 685     }
 686   }
 687   //__________________________________________________________________
 688   /**
 689    * Plot the asymmetries
 690    *
 691    * @param ratios  Asymmetries to plot (if any)
 692    * @param max     Maximum diviation from 1
 693    * @param y1      Lower y coordinate of pad
 694    * @param y2      Upper y coordinate of pad
 695    */
 696   void PlotLeftRight(THStack* leftright, Double_t max,
 697                      Double_t y1, Double_t y2)
 698   {
 699     if (!leftright) return;
 700     bool isBottom = (y1 < 0.0001);
 701     Double_t yd = y2 - y1;
 702     // Make a sub-pad for the result itself
 703     TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
 704     p3->SetTopMargin(0.001);
 705     p3->SetRightMargin(0.05);
 706     p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
 707     p3->SetGridx();
 708     p3->SetTicks(1,1);
 709     p3->SetNumber(2);
 710     p3->Draw();
 711     p3->cd();
 712
 713     TH1* dummy = 0;
 714     if (leftright->GetHists()->GetEntries() == 1) {
 715       // Add dummy histogram
 716       dummy = new TH1F("dummy","", 10, -6, 6);
 717       dummy->SetLineColor(0);
 718       dummy->SetFillColor(0);
 719       dummy->SetMarkerColor(0);
 720       leftright->Add(dummy);
 721     }
 722
 723     // Fix up axis
 724     FixAxis(leftright, 1/yd/1.7, "Right/Left", 4);
 725
 726     leftright->SetMaximum(1+TMath::Max(.12,1.05*max));
 727     leftright->SetMinimum(1-TMath::Max(.15,1.05*max));
 728     p3->Clear();
 729     leftright->DrawClone("nostack e1");
 730
 731
 732     // Make a legend
 733     TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
 734     l2->SetNColumns(2);
 735     l2->SetFillColor(0);
 736     l2->SetFillStyle(0);
 737     l2->SetBorderSize(0);
 738     l2->SetTextFont(132);
 739 #ifndef __CINT__
 740     if (dummy) {
 741       TList* prims = l2->GetListOfPrimitives();
 742       TIter next(prims);
 743       TLegendEntry* o = 0;
 744       while ((o = static_cast<TLegendEntry*>(next()))) {
 745         TString lbl(o->GetLabel());
 746         if (lbl != "dummy") continue;
 747         prims->Remove(o);
 748         break;
 749       }
 750     }
 751 #endif
 752     // Make a nice band from 0.9 to 1.1
 753     TGraphErrors* band = new TGraphErrors(2);
 754     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
 755     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
 756     band->SetPointError(0, 0, .05);
 757     band->SetPointError(1, 0, .05);
 758     band->SetFillColor(kYellow+2);
 759     band->SetFillStyle(3002);
 760     band->SetLineStyle(2);
 761     band->SetLineWidth(1);
 762     band->Draw("3 same");
 763     band->DrawClone("X L same");
 764
 765     leftright->DrawClone("nostack e1 same");
 766     if (isBottom) {
 767       fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName());
 768       p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
 769                                 fRangeParam));
 770     }
 771     else {
 772       fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName());
 773       fRangeParam->fSlave2Pad  = p3;
 774     }
 775   }
 776   /** @} */
 777   //==================================================================
 778   /**
 779    * @{
 780    * @name Data utility functions
 781    */
 782   /**
 783    * Get a result from the passed list
 784    *
 785    * @param list List to search
 786    * @param name Object name to search for
 787    *
 788    * @return
 789    */
 790   TH1* FetchResult(const TList* list, const char* name) const
 791   {
 792     if (!list) return 0;
 793
 794     TList* all = static_cast<TList*>(list->FindObject("all"));
 795     if (!all) {
 796       Warning("GetResult", "No 'all' list find in %s", list->GetName());
 797       // list->ls();
 798       return 0;
 799     }
 800
 801     TH1* ret = static_cast<TH1*>(all->FindObject(name));
 802     if (!ret) {
 803       // all->ls();
 804       Warning("GetResult", "Histogram %s not found", name);
 805     }
 806
 807     return ret;
 808   }
 809   //__________________________________________________________________
 810   /**
 811    * Add a histogram to the stack after possibly rebinning it
 812    *
 813    * @param stack   Stack to add to
 814    * @param hist    histogram
 815    * @param option  Draw options
 816    *
 817    * @return Maximum of histogram
 818    */
 819   Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const
 820   {
 821     // Check if we have input
 822     if (!hist) return 0;
 823
 824     // Rebin if needed
 825     Rebin(hist);
 826
 827     stack->Add(hist, option);
 828     return hist->GetMaximum();
 829   }
 830   //__________________________________________________________________
 831   /**
 832    * Add a histogram to the stack after possibly rebinning it
 833    *
 834    * @param stack   Stack to add to
 835    * @param hist    histogram
 836    * @param option  Draw options
 837    * @param sym     On return, the data symmetriced (added to stack)
 838    *
 839    * @return Maximum of histogram
 840    */
 841   Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option,
 842                         TH1*& sym) const
 843   {
 844     // Check if we have input
 845     if (!hist) return 0;
 846
 847     // Rebin if needed
 848     Rebin(hist);
 849     stack->Add(hist, option);
 850
 851     // Now symmetrice the histogram
 852     sym = Symmetrice(hist);
 853     stack->Add(sym, option);
 854
 855     return hist->GetMaximum();
 856   }
 857
 858   //__________________________________________________________________
 859   /**
 860    * Rebin a histogram
 861    *
 862    * @param h     Histogram to rebin
 863    * @param rebin Rebinning factor
 864    *
 865    * @return
 866    */
 867   virtual void Rebin(TH1* h) const
 868   {
 869     if (fRebin <= 1) return;
 870
 871     Int_t nBins = h->GetNbinsX();
 872     if(nBins % fRebin != 0) {
 873       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
 874               "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
 875       return;
 876     }
 877
 878     // Make a copy
 879     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
 880     tmp->Rebin(fRebin);
 881     tmp->SetDirectory(0);
 882     tmp->Reset();
 883     // The new number of bins
 884     Int_t nBinsNew = nBins / fRebin;
 885     for(Int_t i = 1;i<= nBinsNew; i++) {
 886       Double_t content = 0;
 887       Double_t sumw    = 0;
 888       Double_t wsum    = 0;
 889       Int_t    nbins   = 0;
 890       for(Int_t j = 1; j<=fRebin;j++) {
 891         Int_t    bin = (i-1)*fRebin + j;
 892         Double_t c   =  h->GetBinContent(bin);
 893
 894         if (c <= 0) continue;
 895
 896         if (fCutEdges) {
 897           if (h->GetBinContent(bin+1)<=0 ||
 898               h->GetBinContent(bin-1)) {
 899             Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
 900                     bin, c, h->GetName(),
 901                     bin+1, h->GetBinContent(bin+1),
 902                     bin-1, h->GetBinContent(bin-1));
 903             continue;
 904           }
 905         }
 906         Double_t e =  h->GetBinError(bin);
 907         Double_t w =  1 / (e*e); // 1/c/c
 908         content    += c;
 909         sumw       += w;
 910         wsum       += w * c;
 911         nbins++;
 912       }
 913
 914       if(content > 0 && nbins > 1 ) {
 915         tmp->SetBinContent(i, wsum / sumw);
 916         tmp->SetBinError(i,1./TMath::Sqrt(sumw));
 917       }
 918     }
 919
 920     // Finally, rebin the histogram, and set new content
 921     h->Rebin(fRebin);
 922     h->Reset();
 923     for(Int_t i = 1; i<= nBinsNew; i++) {
 924       h->SetBinContent(i,tmp->GetBinContent(i));
 925       h->SetBinError(i,  tmp->GetBinError(i));
 926     }
 927
 928     delete tmp;
 929   }
 930   //__________________________________________________________________
 931   /**
 932    * Make an extension of @a h to make it symmetric about 0
 933    *
 934    * @param h Histogram to symmertrice
 935    *
 936    * @return Symmetric extension of @a h
 937    */
 938   TH1* Symmetrice(const TH1* h) const
 939   {
 940     Int_t nBins = h->GetNbinsX();
 941     TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
 942     s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
 943     s->Reset();
 944     s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
 945     s->SetMarkerColor(h->GetMarkerColor());
 946     s->SetMarkerSize(h->GetMarkerSize());
 947     s->SetMarkerStyle(h->GetMarkerStyle()+4);
 948     s->SetFillColor(h->GetFillColor());
 949     s->SetFillStyle(h->GetFillStyle());
 950     s->SetDirectory(0);
 951
 952     // Find the first and last bin with data
 953     Int_t first = nBins+1;
 954     Int_t last  = 0;
 955     for (Int_t i = 1; i <= nBins; i++) {
 956       if (h->GetBinContent(i) <= 0) continue;
 957       first = TMath::Min(first, i);
 958       last  = TMath::Max(last,  i);
 959     }
 960
 961     Double_t xfirst = h->GetBinCenter(first-1);
 962     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 963     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 964     for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 965       s->SetBinContent(j, h->GetBinContent(i));
 966       s->SetBinError(j, h->GetBinError(i));
 967     }
 968     // Fill in overlap bin
 969     s->SetBinContent(l2+1, h->GetBinContent(first));
 970     s->SetBinError(l2+1, h->GetBinError(first));
 971     return s;
 972   }
 973   //__________________________________________________________________
 974   /**
 975    * Calculate the left-right asymmetry of input histogram
 976    *
 977    * @param h   Input histogram
 978    * @param max On return, the maximum distance from 1 of the histogram
 979    *
 980    * @return Asymmetry
 981    */
 982   TH1* Asymmetry(TH1* h, Double_t& max)
 983   {
 984     if (!h) return 0;
 985
 986     TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
 987     // Int_t    oBins = h->GetNbinsX();
 988     // Double_t high  = h->GetXaxis()->GetXmax();
 989     // Double_t low   = h->GetXaxis()->GetXmin();
 990     // Double_t dBin  = (high - low) / oBins;
 991     // Int_t    tBins = Int_t(2*high/dBin+.5);
 992     // ret->SetBins(tBins, -high, high);
 993     ret->Reset();
 994     ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
 995     ret->SetYTitle("Right/Left");
 996     Int_t nBins = h->GetNbinsX();
 997     for (Int_t i = 1; i <= nBins; i++) {
 998       Double_t x = h->GetBinCenter(i);
 999       if (x > 0) break;
1000
1001       Double_t c1 = h->GetBinContent(i);
1002       Double_t e1 = h->GetBinError(i);
1003       if (c1 <= 0) continue;
1004
1005       Int_t    j  = h->FindBin(-x);
1006       if (j <= 0 || j > nBins) continue;
1007
1008       Double_t c2 = h->GetBinContent(j);
1009       Double_t e2 = h->GetBinError(j);
1010
1011       Double_t c12 = c1*c1;
1012       Double_t e   = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1013
1014       Int_t    k   = ret->FindBin(x);
1015       ret->SetBinContent(k, c2/c1);
1016       ret->SetBinError(k, e);
1017     }
1018     max = TMath::Max(max, RatioMax(ret));
1019
1020     return ret;
1021   }
1022   //__________________________________________________________________
1023   /**
1024    * Transform a graph into a histogram
1025    *
1026    * @param g
1027    *
1028    * @return
1029    */
1030   TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1031   {
1032     Int_t    nBins = g->GetN();
1033     TArrayF  bins(nBins+1);
1034     Double_t dx = 0;
1035     for (Int_t i = 0; i < nBins; i++) {
1036       Double_t x   = g->GetX()[i];
1037       Double_t exl = g->GetEXlow()[i];
1038       Double_t exh = g->GetEXhigh()[i];
1039       bins.fArray[i]   = x-exl;
1040       bins.fArray[i+1] = x+exh;
1041       Double_t dxi = exh+exl;
1042       if (i == 0) dx  = dxi;
1043       else if (dxi != dx) dx = 0;
1044     }
1045     TString name(g->GetName());
1046     TString title(g->GetTitle());
1047     TH1D* h = 0;
1048     if (dx != 0) {
1049       h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1050     }
1051     else {
1052       h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1053     }
1054     h->SetMarkerStyle(g->GetMarkerStyle());
1055     h->SetMarkerColor(g->GetMarkerColor());
1056     h->SetMarkerSize(g->GetMarkerSize());
1057
1058     return h;
1059   }
1060   /* @} */
1061   //==================================================================
1062   /**
1063    * @{
1064    * @name Ratio utility functions
1065    */
1066   /**
1067    * Get the maximum diviation from 1 in the passed ratio
1068    *
1069    * @param h Ratio histogram
1070    *
1071    * @return Max diviation from 1
1072    */
1073   Double_t RatioMax(TH1* h) const
1074   {
1075     Int_t    nBins = h->GetNbinsX();
1076     Double_t ret   = 0;
1077     for (Int_t i = 1; i <= nBins; i++) {
1078       Double_t c = h->GetBinContent(i);
1079       if (c == 0) continue;
1080       Double_t e = h->GetBinError(i);
1081       Double_t d = TMath::Abs(1-c-e);
1082       ret        = TMath::Max(d, ret);
1083     }
1084     return ret;
1085   }
1086   //__________________________________________________________________
1087   /**
1088    * Make the ratio of h1 to h2
1089    *
1090    * @param h1 First histogram (numerator)
1091    * @param h2 Second histogram (denominator)
1092    *
1093    * @return h1 / h2
1094    */
1095   TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1096   {
1097     const TH1* h1 = dynamic_cast<const TH1*>(o1);
1098     if (h1) {
1099       const TH1* h2 = dynamic_cast<const TH1*>(o2);
1100       if (h2) return Ratio(h1,h2,max);
1101
1102       const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1103       if (g2) return Ratio(h1,g2,max);
1104     }
1105
1106     const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1107     if (g1) {
1108       const TGraphAsymmErrors* g2 = dynamic_cast<const TGraphAsymmErrors*>(o2);
1109       if (g2) return Ratio(g1, g2, max);
1110     }
1111
1112     Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)",
1113             o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
1114     return 0;
1115   }
1116   //__________________________________________________________________
1117   /**
1118    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
1119    * centers of @a h
1120    *
1121    * @param h  Numerator
1122    * @param g  Divisor
1123    *
1124    * @return h/g
1125    */
1126   TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const
1127   {
1128     if (!h || !g) return 0;
1129
1130     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
1131     ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
1132     ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
1133     ret->Reset();
1134     ret->SetMarkerStyle(g->GetMarkerStyle());
1135     ret->SetMarkerColor(h->GetMarkerColor());
1136     ret->SetMarkerSize(0.9*g->GetMarkerSize());
1137     Double_t xlow  = g->GetX()[0];
1138     Double_t xhigh = g->GetX()[g->GetN()-1];
1139     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1140
1141     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1142       Double_t c = h->GetBinContent(i);
1143       if (c <= 0) continue;
1144
1145       Double_t x = h->GetBinCenter(i);
1146       if (x < xlow || x > xhigh) continue;
1147
1148       Double_t f = g->Eval(x);
1149       if (f <= 0) continue;
1150
1151       ret->SetBinContent(i, c / f);
1152       ret->SetBinError(i, h->GetBinError(i) / f);
1153     }
1154     if (ret->GetEntries() <= 0) {
1155       delete ret;
1156       ret = 0;
1157     }
1158     else
1159       max = TMath::Max(RatioMax(ret), max);
1160     return ret;
1161   }
1162   //__________________________________________________________________
1163   /**
1164    * Make the ratio of h1 to h2
1165    *
1166    * @param h1 First histogram (numerator)
1167    * @param h2 Second histogram (denominator)
1168    *
1169    * @return h1 / h2
1170    */
1171   TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const
1172   {
1173     if (!h1 || !h2) return 0;
1174     TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s",
1175                                                h1->GetName(),
1176                                                h2->GetName())));
1177     t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
1178     t1->Divide(h2);
1179     t1->SetMarkerColor(h1->GetMarkerColor());
1180     t1->SetMarkerStyle(h2->GetMarkerStyle());
1181     max = TMath::Max(RatioMax(t1), max);
1182     return t1;
1183   }
1184   //__________________________________________________________________
1185   /**
1186    * Calculate the ratio of two graphs - g1 / g2
1187    *
1188    * @param g1 Numerator
1189    * @param g2 Denominator
1190    *
1191    * @return g1 / g2 in a histogram
1192    */
1193   TH1* Ratio(const TGraphAsymmErrors* g1,
1194              const TGraphAsymmErrors* g2, Double_t& max) const
1195   {
1196     Int_t    nBins = g1->GetN();
1197     TArrayF  bins(nBins+1);
1198     Double_t dx   = 0;
1199     for (Int_t i = 0; i < nBins; i++) {
1200       Double_t x   = g1->GetX()[i];
1201       Double_t exl = g1->GetEXlow()[i];
1202       Double_t exh = g1->GetEXhigh()[i];
1203       bins.fArray[i]   = x-exl;
1204       bins.fArray[i+1] = x+exh;
1205       Double_t dxi = exh+exl;
1206       if (i == 0) dx  = dxi;
1207       else if (dxi != dx) dx = 0;
1208     }
1209     TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
1210     TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
1211     TH1* h = 0;
1212     if (dx != 0) {
1213       h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1214     }
1215     else {
1216       h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
1217     }
1218     h->SetMarkerStyle(g2->GetMarkerStyle());
1219     h->SetMarkerColor(g1->GetMarkerColor());
1220     h->SetMarkerSize(0.9*g2->GetMarkerSize());
1221
1222     Double_t low  = g2->GetX()[0];
1223     Double_t high = g2->GetX()[g2->GetN()-1];
1224     if (low > high) { Double_t t = low; low = high; high = t; }
1225     for (Int_t i = 0; i < nBins; i++) {
1226       Double_t x  = g1->GetX()[i];
1227       if (x < low-dx || x > high+dx) continue;
1228       Double_t c1 = g1->GetY()[i];
1229       Double_t e1 = g1->GetErrorY(i);
1230       Double_t c2 = g2->Eval(x);
1231
1232       h->SetBinContent(i+1, c1 / c2);
1233       h->SetBinError(i+1, e1 / c2);
1234     }
1235     max = TMath::Max(RatioMax(h), max);
1236     return h;
1237   }
1238   /* @} */
1239   //==================================================================
1240   /**
1241    * @{
1242    * @name Graphics utility functions
1243    */
1244   /**
1245    * Find an X axis in a pad
1246    *
1247    * @param p     Pad
1248    * @param name  Histogram to find axis for
1249    *
1250    * @return Found axis or null
1251    */
1252   TAxis* FindXAxis(TVirtualPad* p, const char* name)
1253   {
1254     TObject* o = p->GetListOfPrimitives()->FindObject(name);
1255     if (!o) {
1256       Warning("FindXAxis", "%s not found in pad", name);
1257       return 0;
1258     }
1259     THStack* stack = dynamic_cast<THStack*>(o);
1260     if (!stack) {
1261       Warning("FindXAxis", "%s is not a THStack", name);
1262       return 0;
1263     }
1264     if (!stack->GetHistogram()) {
1265       Warning("FindXAxis", "%s has no histogram", name);
1266       return 0;
1267     }
1268     TAxis* ret = stack->GetHistogram()->GetXaxis();
1269     return ret;
1270   }
1271
1272   //__________________________________________________________________
1273   /**
1274    * Fix the apperance of the axis in a stack
1275    *
1276    * @param stack  stack of histogram
1277    * @param s      Scaling factor
1278    * @param ytitle Y axis title
1279    * @param force  Whether to draw the stack first or not
1280    * @param ynDiv  Divisions on Y axis
1281    */
1282   void FixAxis(THStack* stack, Double_t s, const char* ytitle,
1283                Int_t ynDiv=210, Bool_t force=true)
1284   {
1285     if (force) stack->Draw("nostack e1");
1286
1287     TH1* h = stack->GetHistogram();
1288     if (!h) return;
1289
1290     h->SetXTitle("#eta");
1291     h->SetYTitle(ytitle);
1292     TAxis* xa = h->GetXaxis();
1293     TAxis* ya = h->GetYaxis();
1294     if (xa) {
1295       xa->SetTitle("#eta");
1296       // xa->SetTicks("+-");
1297       xa->SetTitleSize(s*xa->GetTitleSize());
1298       xa->SetLabelSize(s*xa->GetLabelSize());
1299       xa->SetTickLength(s*xa->GetTickLength());
1300     }
1301     if (ya) {
1302       ya->SetTitle(ytitle);
1303       ya->SetDecimals();
1304       // ya->SetTicks("+-");
1305       ya->SetNdivisions(ynDiv);
1306       ya->SetTitleSize(s*ya->GetTitleSize());
1307       ya->SetTitleOffset(ya->GetTitleOffset()/s);
1308       ya->SetLabelSize(s*ya->GetLabelSize());
1309     }
1310   }
1311   /* @} */
1312
1313
1314
1315   //__________________________________________________________________
1316   Bool_t      fShowOthers;   // Show other data
1317   Bool_t      fShowAlice;    // Show ALICE published data
1318   Bool_t      fShowRatios;   // Show ratios
1319   Bool_t      fShowLeftRight;// Show asymmetry
1320   UShort_t    fRebin;        // Rebinning factor
1321   Bool_t      fCutEdges;     // Whether to cut edges
1322   TString     fTitle;        // Title on plot
1323   TNamed*     fTrigString;   // Trigger string (read, or set)
1324   TNamed*     fSNNString;    // Energy string (read, or set)
1325   TNamed*     fSysString;    // Collision system string (read or set)
1326   TAxis*      fVtxAxis;      // Vertex cuts (read or set)
1327   TList       fNumerators;   // Numerators in ratios
1328   TList       fDenominators; // Denominators in ratios
1329   THStack*    fResults;      // Stack of results
1330   THStack*    fRatios;       // Stack of ratios
1331   // TH1*        fForward;      // Results
1332   // TH1*        fForwardMC;    // MC results
1333   // TH1*        fForwardHHD;   // Old results
1334   // TH1*        fTruth;        // MC truth
1335   // TH1*        fCentral;      // Central region data
1336   // TH1*        fForwardSym;   // Symmetric extension
1337   // TH1*        fForwardMCSym; // Symmetric extension
1338   // TH1*        fForwardHHDSym;// Symmetric extension
1339   TH1*        fTriggers;     // Number of triggers
1340   RangeParam* fRangeParam;   // Parameter object for range zoom
1341
1342 };
1343
1344 //=== Stuff for auto zooming =========================================
1345 void UpdateRange(dNdetaDrawer::RangeParam* p)
1346 {
1347   if (!p) {
1348     Warning("UpdateRange", "No parameters %p", p);
1349     return;
1350   }
1351   if (!p->fMasterAxis) {
1352     Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
1353     return;
1354   }
1355   Int_t    first = p->fMasterAxis->GetFirst();
1356   Int_t    last  = p->fMasterAxis->GetLast();
1357   Double_t x1    = p->fMasterAxis->GetBinCenter(first);
1358   Double_t x2    = p->fMasterAxis->GetBinCenter(last);
1359   //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
1360
1361   if (p->fSlave1Axis) {
1362     Int_t i1 = p->fSlave1Axis->FindBin(x1);
1363     Int_t i2 = p->fSlave1Axis->FindBin(x2);
1364     p->fSlave1Axis->SetRange(i1, i2);
1365     p->fSlave1Pad->Modified();
1366     p->fSlave1Pad->Update();
1367   }
1368   if (p->fSlave2Axis) {
1369     Int_t i1 = p->fSlave2Axis->FindBin(x1);
1370     Int_t i2 = p->fSlave2Axis->FindBin(x2);
1371     p->fSlave2Axis->SetRange(i1, i2);
1372     p->fSlave2Pad->Modified();
1373     p->fSlave2Pad->Update();
1374   }
1375   TCanvas*  c = gPad->GetCanvas();
1376   c->cd();
1377 }
1378
1379 //____________________________________________________________________
1380 void RangeExec(dNdetaDrawer::RangeParam* p)
1381 {
1382   // Event types:
1383   //  51:   Mouse move
1384   //  53:
1385   //  1:    Button down
1386   //  21:   Mouse drag
1387   //  11:   Mouse release
1388   // dNdetaDrawer::RangeParam* p =
1389   //   reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
1390   Int_t event     = gPad->GetEvent();
1391   TObject *select = gPad->GetSelected();
1392   if (event == 53) {
1393     UpdateRange(p);
1394     return;
1395   }
1396   if (event != 11 || !select || select != p->fMasterAxis) return;
1397   UpdateRange(p);
1398 }
1399
1400 //=== Steering function ==============================================
1401 void
1402 DrawdNdeta(const char* filename="forward_dndeta.root",
1403            Int_t       flags=0xf,
1404            const char* title="",
1405            UShort_t    rebin=5,
1406            Bool_t      cutEdges=false,
1407            UShort_t    sNN=0,
1408            UShort_t    sys=0,
1409            UShort_t    trg=1,
1410            Float_t     vzMin=999,
1411            Float_t     vzMax=-999)
1412 {
1413   dNdetaDrawer d;
1414   d.SetRebin(rebin);
1415   d.SetCutEdges(cutEdges);
1416   d.SetTitle(title);
1417   d.SetShowOthers(flags & 0x1);
1418   d.SetShowAlice(flags & 0x2);
1419   d.SetShowRatios(flags & 0x4);
1420   d.SetShowLeftRight(flags & 0x8);
1421   // Do the below if your input data does not contain these settings
1422   if (sNN > 0) d.SetSNN(sNN);     // Collision energy per nucleon pair (GeV)
1423   if (sys > 0) d.SetSys(sys);     // Collision system (1:pp, 2:PbPB)
1424   if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
1425   if (vzMin < 999 && vzMax > -999)
1426     d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
1427   d.Run(filename);
1428 }
1429 //____________________________________________________________________
1430 //
1431 // EOF
1432 //
1433