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