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