PWG2/FORWARD/analysis2/DrawdNdeta.C

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