PWGLF/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 <i>very</i> long - sigh - the joy of drawing
  12  * things nicely in ROOT
  13  *
  14  * @ingroup pwglf_forward_dndeta
  15  */
  16 #include <TH1.h>
  17 #include <TColor.h>
  18 #include <THStack.h>
  19 #include <TGraphErrors.h>
  20 #include <TGraphAsymmErrors.h>
  21 #include <TMultiGraph.h>
  22 #include <TFile.h>
  23 #include <TList.h>
  24 #include <TString.h>
  25 #include <TError.h>
  26 #include <TSystem.h>
  27 #include <TROOT.h>
  28 #include <TMath.h>
  29 #include <TCanvas.h>
  30 #include <TPad.h>
  31 #include <TStyle.h>
  32 #include <TLegend.h>
  33 #include <TLegendEntry.h>
  34 #include <TLatex.h>
  35 #include <TImage.h>
  36 #include <TRandom.h>
  37 #include <TParameter.h>
  38 #include <fstream>
  39 #include <iostream>
  40 /** Systematic error color */
  41 #define SYSERR_COLOR kBlue-10
  42 /** Systematic error style */
  43 #define SYSERR_STYLE 1001
  44
  45 Double_t myFunc(Double_t* xp, Double_t* pp);
  46
  47 /**
  48  * Class to draw dN/deta results
  49  *
  50  * @ingroup pwglf_forward_tasks_dndeta
  51  * @ingroup pwglf_forward_dndeta
  52  */
  53 struct dNdetaDrawer
  54 {
  55   /**
  56    * POD of data for range zooming
  57    */
  58   struct RangeParam
  59   {
  60     TAxis*       fMasterAxis; // Master axis
  61     TAxis*       fSlave1Axis; // First slave axis
  62     TVirtualPad* fSlave1Pad;  // First slave pad
  63     TAxis*       fSlave2Axis; // Second slave axis
  64     TVirtualPad* fSlave2Pad;  // Second slave pad
  65   };
  66   //__________________________________________________________________
  67   /**
  68    * Constructor
  69    *
  70    */
  71   dNdetaDrawer()
  72     : // Options
  73     fShowRatios(false),    // Show ratios
  74     fShowLeftRight(false), // Show asymmetry
  75     fShowRings(false),     // Show rings too
  76     fExport(false),        // Export data to script
  77     fCutEdges(false),      // Whether to cut edges
  78     fRemoveOuters(false),  // Whether to remove outers
  79     fShowOthers(0),        // Show other data
  80     fMirror(false),
  81     fForceMB(false),
  82     fAddExec(false),
  83     fSaveAsPDF(false),
  84     fSaveAsPNG(false),
  85     fSaveAsROOT(false),
  86     fSaveAsScript(false),
  87     // Settings
  88     fRebin(0),             // Rebinning factor
  89     fFwdSysErr(0.076),     // Systematic error in forward range
  90     fCenSysErr(0),         // Systematic error in central range
  91     fTitle(""),            // Title on plot
  92     fBase(""),             // Optional base name of output files
  93     fClusterScale(""),     // Scaling of clusters to tracklets
  94     // Read (or set) information
  95     fTrigString(0),        // Trigger string (read, or set)
  96     fNormString(0),        // Normalisation string (read, or set)
  97     fSNNString(0),         // Energy string (read, or set)
  98     fSysString(0),         // Collision system string (read or set)
  99     fVtxAxis(0),           // Vertex cuts (read or set)
 100     fCentAxis(0),          // Centrality axis
 101     fTriggerEff(1),        // Trigger efficency
 102     fExtTriggerEff(false), // True if fTriggerEff was read
 103     // Resulting plots
 104     fResults(0),           // Stack of results
 105     fRatios(0),            // Stack of ratios
 106     fLeftRight(0),         // Left-right asymmetry
 107     fOthers(0),            // Older data
 108     fTriggers(0),          // Number of triggers
 109     fTruth(0),             // Pointer to truth
 110     fRangeParam(0)         // Parameter object for range zoom
 111   {
 112     fRangeParam = new RangeParam;
 113     fRangeParam->fMasterAxis = 0;
 114     fRangeParam->fSlave1Axis = 0;
 115     fRangeParam->fSlave1Pad  = 0;
 116     fRangeParam->fSlave2Axis = 0;
 117     fRangeParam->fSlave2Pad  = 0;
 118   }
 119   /**
 120    * Cpoy constructor
 121    */
 122   dNdetaDrawer(const dNdetaDrawer&) {}
 123   /**
 124    * Assignment operator
 125    *
 126    *
 127    * @return Reference to this object
 128    */
 129   dNdetaDrawer& operator=(const dNdetaDrawer&) { return *this; }
 130
 131   //__________________________________________________________________
 132   /**
 133    * Destructor
 134    */
 135   virtual ~dNdetaDrawer()
 136   {
 137     if (fRatios  && fRatios->GetHists())  fRatios->GetHists()->Delete();
 138     if (fResults && fResults->GetHists()) fResults->GetHists()->Delete();
 139
 140     if (fTrigString) { delete fTrigString; fTrigString = 0; }
 141     if (fSNNString)  { delete fSNNString;  fSNNString  = 0; }
 142     if (fSysString)  { delete fSysString;  fSysString  = 0; }
 143     if (fVtxAxis)    { delete fVtxAxis;    fVtxAxis    = 0; }
 144     if (fCentAxis)   { delete fCentAxis;   fCentAxis   = 0; }
 145     if (fResults)    { delete fResults;    fResults    = 0; }
 146     if (fRatios)     { delete fRatios;     fRatios     = 0; }
 147     if (fOthers)     { delete fOthers;     fOthers     = 0; }
 148     if (fTriggers)   { delete fTriggers;   fTriggers   = 0; }
 149     fRangeParam = 0;
 150   }
 151
 152   //==================================================================
 153   /**
 154    * @{
 155    * @name Set parameters
 156    */
 157   void SetOld(Bool_t use=true) { fOld = use; }
 158   /**
 159    * Show other (UA5, CMS, ...) data
 160    *
 161    * @param x Whether to show or not
 162    */
 163   void SetShowOthers(UShort_t x)    { fShowOthers = x; }
 164   //__________________________________________________________________
 165   /**
 166    * Whether to show ratios or not.  If there's nothing to compare to,
 167    * the ratio panel will be implicitly disabled
 168    *
 169    * @param x Whether to show or not
 170    */
 171   void SetShowRatios(Bool_t x)    { fShowRatios = x; }
 172   //__________________________________________________________________
 173   /**
 174    *
 175    * Whether to show the left/right asymmetry
 176    *
 177    * @param x To show or not
 178    */
 179   void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
 180   //__________________________________________________________________
 181   /**
 182    * Whether to show rings
 183    *
 184    * @param x To show or not
 185    */
 186   void SetShowRings(Bool_t x) { fShowRings = x; }
 187   //__________________________________________________________________
 188   /**
 189    * Whether to export results to a script
 190    *
 191    * @param x Wheter to export results to a script
 192    */
 193   void SetExport(Bool_t x)     { fExport = x; }
 194   //__________________________________________________________________
 195   /**
 196    * Set the rebinning factor
 197    *
 198    * @param x Rebinning factor (must be a divisor in the number of bins)
 199    */
 200   void SetRebin(UShort_t x)       { fRebin = x; }
 201   //__________________________________________________________________
 202   /**
 203    * Wheter to cut away the edges
 204    *
 205    * @param x Whether or not to cut away edges
 206    */
 207   void SetCutEdges(Bool_t x)      { fCutEdges = x; }
 208   //__________________________________________________________________
 209   /**
 210    * Set the title of the plot
 211    *
 212    * @param x Title
 213    */
 214   void SetTitle(TString x)        { fTitle = x; }
 215   //__________________________________________________________________
 216   /**
 217    * Set the base name of the output files
 218    *
 219    * @param x Base name
 220    */
 221   void SetBase(TString x) { fBase = x; }
 222   void SetSaveAsPNG(Bool_t x) { fSaveAsPNG = x; }
 223   void SetSaveAsPDF(Bool_t x) { fSaveAsPDF = x; }
 224   void SetSaveAsROOT(Bool_t x) { fSaveAsROOT = x; }
 225   void SetSaveAsScript(Bool_t x) { fSaveAsScript = x; }
 226   //__________________________________________________________________
 227   /**
 228    * Set the systematic error in the forward region
 229    *
 230    * @param e Systematic error in the forward region
 231    */
 232   void SetForwardSysError(Double_t e=0) { fFwdSysErr = e; }
 233   //__________________________________________________________________
 234   /**
 235    * Set the systematic error in the forward region
 236    *
 237    * @param e Systematic error in the forward region
 238    */
 239   void SetCentralSysError(Double_t e=0) { fCenSysErr = e; }
 240   /**
 241    * Force the plot of minimum bias, even if centrality dependent data
 242    * is present
 243    *
 244    * @param force if true, force minimum bias
 245    */
 246   void SetForceMB(Bool_t force=true) { fForceMB = force; }
 247   /**
 248    * Force the plot of minimum bias, even if centrality dependent data
 249    * is present
 250    *
 251    * @param add if true, force minimum bias
 252    */
 253   void SetAddExec(Bool_t add=true) { fAddExec = add; }
 254   /**
 255    * Mirror data to regions with no coverage (@f$-5.0<\eta<-3.5@f$)
 256    *
 257    * @param mirror If true, mirror data
 258    */
 259   void SetMirror(Bool_t mirror=true) { fMirror = mirror; }
 260   /**
 261    * Set the 'Final MC' correction file.  This is needed if the
 262    * secondary maps where produced using the old code
 263    *
 264    * @param file Filename
 265    */
 266   void SetFinalMC(const TString& file) { fFinalMC = file; }
 267   /**
 268    * Set the file that contains the empirical correction.  This is
 269    * needed when the secondary maps was generated with an in-accurate
 270    * geometry, and when we're analysing data at nominal interaction
 271    * points
 272    *
 273    * @param file Filename
 274    */
 275   void SetEmpirical(const TString& file) { fEmpirical = file; }
 276   /* @} */
 277   //==================================================================
 278   /**
 279    * @{
 280    * @name Override settings from input
 281    */
 282   /**
 283    * Override setting from file
 284    *
 285    * @param sNN Center of mass energy per nucleon pair (GeV)
 286    */
 287   void SetSNN(UShort_t sNN)
 288   {
 289     fSNNString = new TNamed("sNN", Form("%04dGeV", sNN));
 290     fSNNString->SetUniqueID(sNN);
 291   }
 292   //__________________________________________________________________
 293   /**
 294    * Set the collision system
 295    * - 1: pp
 296    * - 2: PbPb
 297    *
 298    * @param sys collision system
 299    */
 300   void SetSys(UShort_t sys)
 301   {
 302     fSysString = new TNamed("sys", (sys == 1 ? "pp" :
 303                                     sys == 2 ? "PbPb" :
 304                                     sys == 3 ? "pPb" :
 305                                     "unknown"));
 306     fSysString->SetUniqueID(sys);
 307   }
 308   //__________________________________________________________________
 309   /**
 310    * Set the vertex range in centimeters
 311    *
 312    * @param vzMin Min @f$ v_z@f$
 313    * @param vzMax Max @f$ v_z@f$
 314    */
 315   void SetVertexRange(Double_t vzMin, Double_t vzMax)
 316   {
 317     fVtxAxis = new TAxis(10, vzMin, vzMax);
 318     fVtxAxis->SetName("vtxAxis");
 319     fVtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", vzMin, vzMax));
 320   }
 321   //__________________________________________________________________
 322   /**
 323    * Set the trigger mask (overrides what's in the file)
 324    *
 325    * @param trig Trigger mask (0x1: INEL, 0x2: INEL>0, 0x4: NSD)
 326    */
 327   void SetTrigger(UShort_t trig)
 328   {
 329     fTrigString = new TNamed("trigString", (trig & 0x1 ? "INEL" :
 330                                             trig & 0x2 ? "INEL>0" :
 331                                             trig & 0x4 ? "NSD" :
 332                                             "unknown"));
 333     fTrigString->SetUniqueID(trig);
 334   }
 335   //__________________________________________________________________
 336   /**
 337    * Set the trigger efficiency - if set, then scale result histograms
 338    * by this factor
 339    *
 340    * @param eff @f$\varepsilon_{T}@f$
 341    */
 342   void SetTriggerEfficiency(Float_t eff)
 343   {
 344     fTriggerEff = eff;
 345     fExtTriggerEff = false;
 346   }
 347
 348   //==================================================================
 349   /**
 350    * @{
 351    * @name Main steering functions
 352    */
 353   /**
 354    * Run the code to produce the final result.
 355    *
 356    * @param filename  File containing the data
 357    */
 358   void Run(const char* filename="forward_dndeta.root")
 359   {
 360     Double_t max = 0, rmax=0, amax=0;
 361
 362     gStyle->SetPalette(1);
 363
 364     // --- Open input file -------------------------------------------
 365     TFile* file = TFile::Open(filename, "READ");
 366     if (!file) {
 367       Error("Run", "Cannot open %s", filename);
 368       return;
 369     }
 370     Info("Run", "Drawing results from %s", file->GetName());
 371
 372     // --- Get forward list ------------------------------------------
 373     TList* forward = static_cast<TList*>(file->Get("ForwardResults"));
 374     if (!forward) {
 375       Error("Run", "Couldn't find list ForwardResults");
 376       return;
 377     }
 378     TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
 379     if (!sums) {
 380       Error("Run", "Couldn't find list ForwardSums");
 381       return;
 382     }
 383     TParameter<bool>* p =
 384       static_cast<TParameter<bool>*>(sums->FindObject("empirical"));
 385     if (p && p->GetVal() && !fEmpirical.IsNull()) {
 386       Warning("Run", "Empirical correction already applied");
 387       fEmpirical = "__task__";
 388     }
 389     // --- Get information on the run --------------------------------
 390     FetchInformation(forward);
 391
 392     // --- Print settings --------------------------------------------
 393     Info("Run", "Settings for the drawer:\n"
 394          "   Show ratios:                      %5s\n"
 395          "   Show Left/right:                  %5s\n"
 396          "   Show rings:                       %5s\n"
 397          "   Export to file:                   %5s\n"
 398          "   Cut edges when rebinning:         %5s\n"
 399          "   Remove outer rings:               %5s\n"
 400          "   Mirror to un-covered regions:     %5s\n"
 401          "   Force minimum bias:               %5s\n"
 402          "   Show other results:               0x%03x\n"
 403          "   Rebinning factor:                 %5d\n"
 404          "   Forward systematic error:         %5.1f%%\n"
 405          "   Central systematic error:         %5.1f%%\n"
 406          "   Title on plot:                    %s\n"
 407          "   Scaling of clusters to tracklets: %s\n"
 408          "   Final MC correction file:         %s\n"
 409          "   Empirical correction file:        %s",
 410          (fShowRatios    ? "yes" : "no"),
 411          (fShowLeftRight ? "yes" : "no"),
 412          (fShowRings     ? "yes" : "no"),
 413          (fExport        ? "yes" : "no"),
 414          (fCutEdges      ? "yes" : "no"),
 415          (fRemoveOuters  ? "yes" : "no"),
 416          (fMirror        ? "yes" : "no"),
 417          (fForceMB       ? "yes" : "no"),
 418          fShowOthers, fRebin, (100*fFwdSysErr), (100*fCenSysErr),
 419          fTitle.Data(), fClusterScale.Data(), fFinalMC.Data(),
 420          fEmpirical.Data());
 421
 422     // --- Set the macro pathand load other data script --------------
 423     TString savPath(gROOT->GetMacroPath());
 424     gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2",
 425                              gROOT->GetMacroPath()));
 426     // Always recompile
 427     if (!gROOT->GetClass("RefData"))
 428       gROOT->LoadMacro("OtherData.C++");
 429     gROOT->SetMacroPath(savPath);
 430
 431     // --- Get the central results -----------------------------------
 432     TList* clusters = static_cast<TList*>(file->Get("CentralResults"));
 433     if (!clusters) Warning("Run", "Couldn't find list CentralResults");
 434
 435     // --- Get the central results -----------------------------------
 436     TList* mcTruth = static_cast<TList*>(file->Get("MCTruthResults"));
 437     if (!mcTruth) Warning("Run", "Couldn't find list MCTruthResults");
 438
 439     // --- Make our containtes ---------------------------------------
 440     fResults   = new THStack("results", "Results");
 441     fRatios    = new THStack("ratios",  "Ratios");
 442     fLeftRight = new THStack("asymmetry", "Left-right asymmetry");
 443     fOthers    = new TMultiGraph();
 444
 445     // --- Try to open the final MC file, and find relevant lists ----
 446     TList* forwardMC = 0;
 447     // TList* centralMC = 0;
 448     if (!fFinalMC.IsNull()) {
 449       TFile* finalMC = TFile::Open(fFinalMC, "READ");
 450       if (!finalMC) {
 451         Warning("Run", "Failed to open file %s for final MC corrections",
 452                 fFinalMC.Data());
 453       }
 454       else {
 455         forwardMC = static_cast<TList*>(finalMC->Get("ForwardResults"));
 456         if (!forwardMC)
 457           Warning("Run", "Couldn't find list ForwardResults for final MC");
 458 #if 0
 459         centralMC = static_cast<TList*>(finalMC->Get("CentralResults"));
 460         if (!centralMC)
 461           Warning("Run", "Couldn't find list CentralResults for final MC");
 462 #endif
 463       }
 464     }
 465     if (!forwardMC) fFinalMC = "";
 466
 467     // --- Try to get the emperical correction -----------------------
 468     TGraphErrors* empCorr = 0;
 469     if (!fEmpirical.IsNull() && !fEmpirical.EqualTo("__task__")) {
 470       if (gSystem->AccessPathName(fEmpirical.Data())) { // Not found here
 471         fEmpirical =
 472           gSystem->ExpandPathName(Form("$ALICE_ROOT/PWGLF/FORWARD/"
 473                                        "corrections/Empirical/%s",
 474                                        fEmpirical.Data()));
 475         if (gSystem->AccessPathName(fEmpirical.Data())) { // Not found here
 476           Warning("Run", "Couldn't get empirical correction file");
 477           fEmpirical = "";
 478         }
 479       }
 480       if (!fEmpirical.IsNull()) {
 481         TFile* empirical = TFile::Open(fEmpirical, "READ");
 482         if (!empirical) {
 483           Warning("Run", "couldn't open empirical correction file: %s",
 484                   fEmpirical.Data());
 485           fEmpirical = "";
 486         }
 487         const char* empPath = "fmdfull/average";
 488         empCorr = static_cast<TGraphErrors*>(empirical->Get(empPath));
 489         if (!empCorr) {
 490           Warning("Run", "Didn't find the graph %s in %s",
 491                   empPath, fEmpirical.Data());
 492           fEmpirical = "";
 493         }
 494       }
 495     }
 496     if (!empCorr && !fEmpirical.EqualTo("__task__")) fEmpirical = "";
 497
 498     // --- Loop over input data --------------------------------------
 499     TObjArray truths;
 500     FetchResults(mcTruth,  0, 0, "MCTruth", max, rmax, amax,truths);
 501     TObjArray* fwdA = FetchResults(forward,  forwardMC, empCorr, "Forward",
 502                                    max, rmax, amax,truths);
 503     TObjArray* cenA = FetchResults(clusters, 0, 0, "Central",
 504                                    max, rmax, amax,truths);
 505
 506     // --- Get trigger information -----------------------------------
 507     // TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
 508     if (sums) {
 509       TList* all = 0;
 510       if (fOld) all = sums;
 511       else      all = static_cast<TList*>(sums->FindObject("all"));
 512       if (all) {
 513         fTriggers = FetchResult(all, "triggers");
 514         if (!fTriggers) all->ls();
 515       }
 516       else  {
 517         Warning("Run", "List all not found in ForwardSums");
 518         sums->ls();
 519       }
 520     }
 521     else {
 522       Warning("Run", "No ForwardSums directory found in %s", file->GetName());
 523       file->ls();
 524     }
 525
 526     // --- Check our stacks ------------------------------------------
 527     if (!fResults->GetHists() ||
 528         fResults->GetHists()->GetEntries() <= 0) {
 529       Error("Run", "No histograms in result stack!");
 530       return;
 531     }
 532     if (!fOthers->GetListOfGraphs() ||
 533         fOthers->GetListOfGraphs()->GetEntries() <= 0) {
 534       Warning("Run", "No other data found - disabling that");
 535       fShowOthers = 0;
 536     }
 537     if (!fRatios->GetHists() ||
 538         fRatios->GetHists()->GetEntries() <= 0) {
 539       Warning("Run", "No ratio data found - disabling that");
 540       // fRatios->ls();
 541       fShowRatios = false;
 542     }
 543     if (!fLeftRight->GetHists() ||
 544         fLeftRight->GetHists()->GetEntries() <= 0) {
 545       Warning("Run", "No left/right data found - disabling that");
 546       // fLeftRight->ls();
 547       fShowLeftRight = false;
 548     }
 549     if (fFwdSysErr > 0) {
 550       if (fCenSysErr <= 0) fCenSysErr = fFwdSysErr;
 551       for (Int_t i = 0; i < fwdA->GetEntriesFast(); i++) {
 552         TH1* fwd = static_cast<TH1*>(fwdA->At(i));
 553         TH1* cen = static_cast<TH1*>(cenA->At(i));
 554         CorrectForward(fwd);
 555         CorrectCentral(cen);
 556         Double_t low, high;
 557         TH1* tmp = Merge(cen, fwd, low, high);
 558         TF1* f   = FitMerged(tmp, low, high);
 559         MakeSysError(tmp, cen, fwd, f);
 560         delete f;
 561         Info("", "Adding systematic error histogram %s",
 562              tmp->GetName());
 563         fResults->GetHists()->AddFirst(tmp, "e5");
 564
 565         if (!fMirror) continue;
 566
 567         TH1* tmp2 = Symmetrice(tmp);
 568         tmp2->SetFillColor(tmp->GetFillColor());
 569         tmp2->SetFillStyle(tmp->GetFillStyle());
 570         tmp2->SetMarkerStyle(tmp->GetMarkerStyle());
 571         tmp2->SetLineWidth(tmp->GetLineWidth());
 572         fResults->GetHists()->AddFirst(tmp2, "e5");
 573         fResults->Modified();
 574       }
 575     }
 576     delete fwdA;
 577     delete cenA;
 578
 579     // --- Close the input file --------------------------------------
 580     file->Close();
 581
 582
 583
 584     // --- Plot results ----------------------------------------------
 585     Plot(max, rmax, amax);
 586   }
 587
 588   //__________________________________________________________________
 589   /**
 590    * Fetch the information on the run from the results list
 591    *
 592    * @param results  Results list
 593    */
 594   void FetchInformation(const TList* results)
 595   {
 596     if (!fTrigString)
 597       fTrigString = static_cast<TNamed*>(results->FindObject("trigger"));
 598     if (!fNormString)
 599       fNormString = static_cast<TNamed*>(results->FindObject("scheme"));
 600     if (!fSNNString)
 601       fSNNString  = static_cast<TNamed*>(results->FindObject("sNN"));
 602     if (!fSysString)
 603       fSysString  = static_cast<TNamed*>(results->FindObject("sys"));
 604     if (!fVtxAxis)
 605       fVtxAxis    = static_cast<TAxis*>(results->FindObject("vtxAxis"));
 606     if (!fCentAxis)
 607       fCentAxis   = static_cast<TAxis*>(results->FindObject("centAxis"));
 608     if (fTriggerEff <= 0 || (1-fTriggerEff)<1e-6) {
 609       TParameter<double>* eff =
 610         static_cast<TParameter<double>*>(results->FindObject("triggerEff"));
 611       if (eff) {
 612         fTriggerEff = eff->GetVal();
 613         fExtTriggerEff = true;
 614       }
 615       if (fTriggerEff <= 0) SetTriggerEfficiency(1);
 616     }
 617
 618     TNamed* options = static_cast<TAxis*>(results->FindObject("options"));
 619     if (!fTrigString) fTrigString = new TNamed("trigger", "unknown");
 620     if (!fNormString) fNormString = new TNamed("scheme", "unknown");
 621     if (!fSNNString)  fSNNString  = new TNamed("sNN", "unknown");
 622     if (!fSysString)  fSysString  = new TNamed("sys", "unknown");
 623     if (!fVtxAxis) {
 624       fVtxAxis    = new TAxis(1,0,0);
 625       fVtxAxis->SetName("vtxAxis");
 626       fVtxAxis->SetTitle("v_{z} range unspecified");
 627     }
 628
 629     TString centTxt("none");
 630     if (fCentAxis) {
 631       Int_t nCent = fCentAxis->GetNbins();
 632       centTxt = Form("%d bins", nCent);
 633       for (Int_t i = 0; i <= nCent; i++)
 634         centTxt.Append(Form("%c%d", i == 0 ? ' ' : '-',
 635                             int(fCentAxis->GetXbins()->At(i))));
 636     }
 637     Info("FetchInformation",
 638          "Initialized for\n"
 639          "   Trigger:       %-30s  (0x%x)\n"
 640          "   Efficiency:    %-6.4f\n"
 641          "   sqrt(sNN):     %-30s  (%dGeV)\n"
 642          "   System:        %-30s  (%d)\n"
 643          "   Vz range:      %-30s  (%f,%f)\n"
 644          "   Normalization: %-30s  (%d)\n"
 645          "   Centrality:    %s\n"
 646          "   Options:       %s",
 647          fTrigString->GetTitle(), fTrigString->GetUniqueID(),
 648          fTriggerEff,
 649          fSNNString->GetTitle(),  fSNNString->GetUniqueID(),
 650          fSysString->GetTitle(),  fSysString->GetUniqueID(),
 651          fVtxAxis->GetTitle(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax(),
 652          fNormString->GetTitle(), fNormString->GetUniqueID(),
 653          centTxt.Data(), (options ? options->GetTitle() : "none"));
 654     if (fSysString->GetUniqueID() == 3) {
 655       Info("FetchResults", "Left/Right assymmetry, mirror, and systematic "
 656            "errors explicitly disabled for pPb");
 657       fShowLeftRight = false;
 658       fMirror        = false;
 659       fFwdSysErr     = 0;
 660       fCenSysErr     = 0;
 661     }
 662   }
 663   //__________________________________________________________________
 664   TMultiGraph* FetchOthers(UShort_t centLow, UShort_t centHigh)
 665   {
 666     TMultiGraph* thisOther = 0;
 667     if (fShowOthers == 0) return 0;
 668
 669     UShort_t sys   = (fSysString  ? fSysString->GetUniqueID() : 0);
 670     UShort_t trg   = (fTrigString ? fTrigString->GetUniqueID() : 0);
 671     UShort_t snn   = (fSNNString  ? fSNNString->GetUniqueID() : 0);
 672     Long_t   ret   =
 673       gROOT->ProcessLine(Form("RefData::GetData(%d,%d,%d,%d,%d,%d);",
 674                               sys,snn,trg,centLow,centHigh,fShowOthers));
 675     if (!ret) return 0;
 676
 677     thisOther = reinterpret_cast<TMultiGraph*>(ret);
 678     return thisOther;
 679   }
 680   //__________________________________________________________________
 681   /**
 682    * Get the results from the top-level list
 683    *
 684    * @param list    List
 685    * @param mcList  List of histograms from MC
 686    * @param empCorr Emperical correction if any
 687    * @param name    name
 688    * @param max     On return, maximum of data
 689    * @param rmax    On return, maximum of ratios
 690    * @param amax    On return, maximum of left-right comparisons
 691    * @param truths  List of MC truths to compare to.
 692    *
 693    * @return Array of results
 694    */
 695   TObjArray*
 696   FetchResults(const TList*  list,
 697                const TList*  mcList,
 698                TGraphErrors* empCorr,
 699                const char*   name,
 700                Double_t&     max,
 701                Double_t&     rmax,
 702                Double_t&     amax,
 703                TObjArray&    truths)
 704   {
 705     if (!list) return 0;
 706     UShort_t   n = HasCent() ? fCentAxis->GetNbins() : 0;
 707     // Info("FetchResults","got %d centrality bins", n);
 708     if (n == 0) {
 709       TH1*  h = FetchOne(list, mcList, empCorr, name, "all",
 710                          FetchOthers(0,0), -1000, 0,
 711                          max, rmax, amax, fTruth);
 712       if (!h) return 0;
 713       TObjArray* a = new TObjArray;
 714       // Info("FetchResults", "Adding %s to result stack", h->GetName());
 715       a->AddAt(h, 0);
 716       return a;
 717     }
 718
 719     TObjArray* a = new TObjArray;
 720     truths.Expand(n);
 721     for (UShort_t i = 0; i < n; i++) {
 722       UShort_t centLow  = fCentAxis->GetBinLowEdge(i+1);
 723       UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1);
 724       TString  lname    = Form("cent%03d_%03d", centLow, centHigh);
 725       Int_t    col      = GetCentralityColor(i+1);
 726       TString  centTxt  = Form("%3d%%-%3d%% central", centLow, centHigh);
 727
 728       TH1* tt = static_cast<TH1*>(truths.At(i));
 729       TH1* ot = tt;
 730       TH1* h  = FetchOne(list, mcList, empCorr, name, lname,
 731                          FetchOthers(centLow,centHigh), col,
 732                          centTxt.Data(), max, rmax, amax, fTruth);
 733       if (!h) continue;
 734       if (ot != tt) {
 735         //Info("FetchResults", "old truth=%p new truth=%p (%s)", ot, tt, name);
 736         truths.AddAt(tt, i);
 737       }
 738       // Info("FetchResults", "Adding %p to result stack", h);
 739       a->AddAt(h, i);
 740     }
 741     return a;
 742   }
 743   //__________________________________________________________________
 744   /**
 745    * Get the color for a centrality bin
 746    *
 747    * @param bin Centrality bin
 748    *
 749    * @return Color
 750    */
 751   Int_t GetCentralityColor(Int_t bin) const
 752   {
 753     if (fCentAxis->GetNbins() < 6) {
 754       switch (bin) {
 755       case 1: return kRed+2;
 756       case 2: return kGreen+2;
 757       case 3: return kBlue+1;
 758       case 4: return kCyan+1;
 759       case 5: return kMagenta+1;
 760       case 6: return kYellow+2;
 761       }
 762     }
 763     UShort_t centLow  = fCentAxis->GetBinLowEdge(bin);
 764     UShort_t centHigh = fCentAxis->GetBinUpEdge(bin);
 765     Float_t  fc       = (centLow+double(centHigh-centLow)/2) / 100;
 766     Int_t    nCol     = gStyle->GetNumberOfColors();
 767     Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
 768     Int_t    col      = gStyle->GetColorPalette(icol);
 769     //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col);
 770     return col;
 771   }
 772   //__________________________________________________________________
 773   /**
 774    * Set attributed on a histogram
 775    *
 776    * @param h     Histogram
 777    * @param color Color
 778    */
 779   void SetAttributes(TH1* h, Int_t color)
 780   {
 781     if (!h) return;
 782     if (color < 0) return;
 783     // h->SetLineColor(color);
 784     h->SetMarkerColor(color);
 785     // h->SetFillColor(color);
 786   }
 787   //__________________________________________________________________
 788   /**
 789    * Set attributed on a graph
 790    *
 791    * @param g     Graph
 792    * @param color Color
 793    */
 794   void SetAttributes(TGraph* g, Int_t color)
 795   {
 796     if (!g) return;
 797     if (color < 0) return;
 798     // g->SetLineColor(color);
 799     g->SetMarkerColor(color);
 800     // g->SetFillColor(color);
 801   }
 802   //__________________________________________________________________
 803   /**
 804    * Modify the title
 805    *
 806    */
 807   void ModifyTitle(TNamed* h, const char* /*centTxt*/)
 808   {
 809     if (!h) return;
 810
 811     TString title(h->GetTitle());
 812     title.ReplaceAll("ALICE ","");
 813     if (title.Contains("Central"))
 814       title.ReplaceAll("Central", "SPD clusters");
 815     if (title.Contains("Forward"))
 816       title.ReplaceAll("Forward", "FMD");
 817     h->SetTitle(title);
 818
 819
 820     return;
 821     // if (!centTxt || !h) return;
 822     // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt));
 823   }
 824   //__________________________________________________________________
 825   TH1* FetchOne(const TList*  list,
 826                 const TList*  mcList,
 827                 TGraphErrors* empCorr,
 828                 const char*   name,
 829                 const char*   folderName,
 830                 TMultiGraph*  others,
 831                 Int_t         col,
 832                 const char*   txt,
 833                 Double_t&     max,
 834                 Double_t&     rmax,
 835                 Double_t&     amax,
 836                 TH1*&         truth)
 837   {
 838     TList* folder = 0;
 839     if (fOld) folder = const_cast<TList*>(list);
 840     else      folder = static_cast<TList*>(list->FindObject(folderName));
 841     if (!folder) {
 842       Error("FetchResults", "Couldn't find list '%s' in %s",
 843             folderName, list->GetName());
 844       return 0;
 845     }
 846     TList* mcFolder = 0;
 847     if (mcList) {
 848       mcFolder = static_cast<TList*>(mcList->FindObject(folderName));
 849       if (!mcFolder)
 850         Warning("FetchResults",
 851                 "Didn't find the list '%s' in %s for final MC correction",
 852                 folderName, mcList->GetName());
 853     }
 854     TObject* normCalc = folder->FindObject("normCalc");
 855     if (normCalc) Info("FetchOne", "%s:\n%s", folderName, normCalc->GetTitle());
 856     TH1* h = FetchResults(folder, mcFolder, empCorr, name,
 857                           others, col, txt, max, rmax, amax, truth);
 858     return h;
 859   }
 860   //__________________________________________________________________
 861   /**
 862    * Fetch results for a particular centrality bin
 863    *
 864    * @param list       List
 865    * @param mcList     List of MC results
 866    * @param empCorr    Emperical correction if any
 867    * @param name       Name
 868    * @param thisOther  Other graphs
 869    * @param color      Color
 870    * @param centTxt    Centrality text
 871    * @param max        On return, data maximum
 872    * @param rmax       On return, ratio maximum
 873    * @param amax       On return, left-right maximum
 874    * @param truth      MC truth to compare to or possibly update
 875    *
 876    * @return Histogram of results
 877    */
 878   TH1* FetchResults(const TList*  list,
 879                     const TList*  mcList,
 880                     TGraphErrors* empCorr,
 881                     const char*   name,
 882                     TMultiGraph*  thisOther,
 883                     Int_t         color,
 884                     const char*   centTxt,
 885                     Double_t&     max,
 886                     Double_t&     rmax,
 887                     Double_t&     amax,
 888                     TH1*&         truth)
 889   {
 890
 891     TH1* dndeta      = FetchResult(list, Form("dndeta%s", name));
 892     TH1* dndetaMC    = FetchResult(list, Form("dndeta%sMC", name));
 893     TH1* dndetaTruth = FetchResult(list, "dndetaTruth");
 894
 895     if (mcList && FetchResult(mcList, "finalMCCorr"))
 896       Warning("FetchResults", "dNdeta already corrected for final MC");
 897     else
 898       CorrectFinalMC(dndeta, mcList);
 899
 900     CorrectEmpirical(dndeta, empCorr);
 901     CorrectTriggerEff(dndeta);
 902     CorrectTriggerEff(dndetaMC);
 903
 904     TH1* dndetaSym   = 0;
 905     TH1* dndetaMCSym = 0;
 906     SetAttributes(dndeta,     color);
 907     SetAttributes(dndetaMC,   HasCent() ? color : color+2);
 908     SetAttributes(dndetaTruth,color);
 909     SetAttributes(dndetaSym,  color);
 910     SetAttributes(dndetaMCSym,HasCent() ? color : color+2);
 911     if (dndetaMC && HasCent())
 912       dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2);
 913     if (dndetaMCSym && HasCent())
 914       dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2);
 915     if (dndetaTruth && HasCent()) {
 916       dndetaTruth->SetMarkerStyle(34);
 917       dndetaTruth->SetMarkerColor(kYellow-1);
 918     }
 919     if (dndetaTruth) {
 920       dndetaTruth->SetLineColor(kBlack);
 921       dndetaTruth->SetFillColor(kBlack);
 922       dndetaTruth->SetFillStyle(3002);
 923       // dndetaTruth->SetLineColor(kBlack);
 924     }
 925     ModifyTitle(dndeta,     centTxt);
 926     ModifyTitle(dndetaMC,   centTxt);
 927     ModifyTitle(dndetaTruth,centTxt);
 928     ModifyTitle(dndetaSym,  centTxt);
 929     ModifyTitle(dndetaMCSym,centTxt);
 930
 931
 932     max = TMath::Max(max, AddHistogram(fResults, dndetaTruth, "e5"));
 933     max = TMath::Max(max, AddHistogram(fResults, dndetaMC,    dndetaMCSym));
 934     max = TMath::Max(max, AddHistogram(fResults, dndeta,      dndetaSym));
 935
 936     if (dndetaTruth) {
 937       truth = dndetaTruth;
 938     }
 939     else {
 940       if (fShowRings) {
 941         THStack* rings = static_cast<THStack*>(list->FindObject("dndetaRings"));
 942         if (rings) {
 943           TIter next(rings->GetHists());
 944           TH1*  hist = 0;
 945           while ((hist = static_cast<TH1*>(next())))
 946             max = TMath::Max(max, AddHistogram(fResults, hist));
 947         }
 948       }
 949       // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f",
 950       //      dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max);
 951
 952       if (fShowLeftRight) {
 953         fLeftRight->Add(Asymmetry(dndeta,    amax));
 954         fLeftRight->Add(Asymmetry(dndetaMC,  amax));
 955       }
 956
 957       if (thisOther) {
 958         TIter next(thisOther->GetListOfGraphs());
 959         TGraph* g = 0;
 960         while ((g = static_cast<TGraph*>(next()))) {
 961           fRatios->Add(Ratio(dndeta,    g, rmax));
 962           fRatios->Add(Ratio(dndetaSym, g, rmax));
 963           SetAttributes(g, color);
 964           ModifyTitle(g, centTxt);
 965           if (!fOthers->GetListOfGraphs() ||
 966               !fOthers->GetListOfGraphs()->FindObject(g->GetName())) {
 967             max = TMath::Max(max,TMath::MaxElement(g->GetN(), g->GetY()));
 968             fOthers->Add(g);
 969           }
 970         }
 971         // fOthers->Add(thisOther);
 972       }
 973     }
 974     if (dndetaMC) {
 975       fRatios->Add(Ratio(dndeta,    dndetaMC,    rmax));
 976       fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax));
 977     }
 978     if (truth) {
 979       fRatios->Add(Ratio(dndeta,      truth, rmax));
 980       fRatios->Add(Ratio(dndetaSym,   truth, rmax));
 981     }
 982     return dndeta;
 983   }
 984   //__________________________________________________________________
 985   void CorrectFinalMC(TH1* dndeta, const TList* mcList)
 986   {
 987     if (!dndeta) return;
 988     if (!mcList) return;
 989
 990     TH1* dndetaMC    = FetchResult(mcList, dndeta->GetName());
 991     TH1* dndetaTruth = FetchResult(mcList, "dndetaTruth");
 992     if (!dndetaMC || !dndetaTruth) return;
 993
 994     TH1* corr = static_cast<TH1*>(dndetaMC->Clone("finalMCCorr"));
 995     corr->Divide(dndetaTruth);
 996
 997     Info("CorrectFinalMC", "Correcting dN/deta with final MC correction");
 998     dndeta->Divide(corr);
 999   }
1000   //__________________________________________________________________
1001   void CorrectEmpirical(TH1* dndeta, const TGraphErrors* empCorr)
1002   {
1003     if (!dndeta) return;
1004     if (!empCorr) return;
1005
1006     Info("CorrectEmpirical", "Doing empirical correction of dN/deta");
1007     TAxis* xAxis = dndeta->GetXaxis();
1008     for (Int_t i = 1; i <= xAxis->GetNbins(); i++) {
1009       Double_t x = xAxis->GetBinCenter(i);
1010       Double_t y = dndeta->GetBinContent(i);
1011       Double_t c = empCorr->Eval(x);
1012       dndeta->SetBinContent(i, y / c);
1013     }
1014   }
1015   //__________________________________________________________________
1016   void CorrectTriggerEff(TH1* dndeta)
1017   {
1018     if (fExtTriggerEff) return;
1019     if (!dndeta) return;
1020     if (fTriggerEff <= 0 || fTriggerEff >= 1) return;
1021     dndeta->Scale(fTriggerEff);
1022   }
1023   //__________________________________________________________________
1024   /**
1025    * Plot the results
1026    * @param max        Max value
1027    * @param rmax       Maximum diviation from 1 of ratios
1028    * @param amax       Maximum diviation from 1 of asymmetries
1029    */
1030   void Plot(Double_t     max,
1031             Double_t     rmax,
1032             Double_t     amax)
1033   {
1034     gStyle->SetOptTitle(0);
1035     gStyle->SetTitleFont(kFont, "xyz");
1036     gStyle->SetLabelFont(kFont, "xyz");
1037
1038     Int_t    h  = 800;
1039     Int_t    w  = 800; // h / TMath::Sqrt(2);
1040     Double_t y1 = 0;
1041     Double_t y2 = 0;
1042     Double_t y3 = 0;
1043     if (!fShowRatios)    w  *= 1.3;
1044     else                 y1 =  0.3;
1045     if (!fShowLeftRight) w  *= 1.3;
1046     else {
1047       Double_t y11 = y1;
1048       y1 = (y11 > 0.0001 ? 0.4 : 0.2);
1049       y2 = (y11 > 0.0001 ? 0.2 : 0.3);
1050     }
1051     TCanvas* c = new TCanvas("Results", "Results", w, h);
1052     c->SetFillColor(0);
1053     c->SetBorderSize(0);
1054     c->SetBorderMode(0);
1055
1056     PlotResults(max, y1);
1057     c->cd();
1058
1059     PlotRatios(rmax, y2, y1);
1060     c->cd( );
1061
1062     PlotLeftRight(amax, y3, y2);
1063     c->cd();
1064
1065
1066     Int_t   vMin = fVtxAxis->GetXmin();
1067     Int_t   vMax = fVtxAxis->GetXmax();
1068     TString trg(fTrigString->GetTitle());
1069     Int_t   nev  = 0;
1070     if (fTriggers) nev = fTriggers->GetBinContent(1);
1071     trg          = trg.Strip(TString::kBoth);
1072     trg.ReplaceAll(" ", "_");
1073     trg.ReplaceAll(">", "Gt");
1074     trg.ReplaceAll("&", "AND");
1075     trg.ReplaceAll("|", "OR");
1076     if (fBase.IsNull())
1077       fBase = "dndeta_<sys>_<snn>_<trig>_<ipmin><ipmax>cm_<nev>ev";
1078     fBase.ReplaceAll("<sys>",   fSysString->GetTitle());
1079     fBase.ReplaceAll("<snn>",   fSNNString->GetTitle());
1080     fBase.ReplaceAll("<trig>",  trg.Data());
1081     fBase.ReplaceAll("<ipmin>", Form("%c%02d",vMin<0?'m':'p',TMath::Abs(vMin)));
1082     fBase.ReplaceAll("<ipmax>", Form("%c%02d",vMax<0?'m':'p',TMath::Abs(vMax)));
1083     fBase.ReplaceAll("<nev>",   Form("%09d",  nev));
1084     if (fSaveAsPNG)   c->SaveAs(Form("%s.png",  fBase.Data()));
1085     if (fSaveAsROOT)  c->SaveAs(Form("%s.root", fBase.Data()));
1086     if (fSaveAsScript)c->SaveAs(Form("%s.C",    fBase.Data()));
1087     if (fSaveAsPDF)   c->SaveAs(Form("%s.pdf",  fBase.Data()));
1088     TString exp(fBase);
1089     exp.ReplaceAll("dndeta", "export");
1090     Export(exp);
1091   }
1092   //__________________________________________________________________
1093   /**
1094    * Build main legend
1095    *
1096    * @param stack    Stack to include
1097    * @param mg       (optional) Multi graph to include
1098    * @param x1       Lower X coordinate in the range [0,1]
1099    * @param y1       Lower Y coordinate in the range [0,1]
1100    * @param x2       Upper X coordinate in the range [0,1]
1101    * @param y2       Upper Y coordinate in the range [0,1]
1102    * @param forceCol If non-zero, force this many columns
1103    */
1104   void BuildLegend(THStack* stack, TMultiGraph* mg,
1105                    Double_t x1, Double_t y1, Double_t x2, Double_t y2,
1106                    Int_t forceCol=0)
1107   {
1108     TLegend* l = new TLegend(x1,y1,x2,y2);
1109     Int_t nCol = forceCol;
1110     if (nCol <= 0) nCol = HasCent() ? 1 : 2;
1111     l->SetNColumns(nCol);
1112     l->SetFillColor(0);
1113     l->SetFillStyle(0);
1114     l->SetBorderSize(0);
1115     l->SetTextFont(kFont);
1116     l->SetTextColor(kAliceBlue);
1117
1118     // Loop over items in stack and get unique items, while ignoring
1119     // mirrored data and systematic error bands
1120     TIter    next(stack->GetHists());
1121     TH1*     hist = 0;
1122     TObjArray unique;
1123     unique.SetOwner();
1124     Bool_t   sysErrSeen = false;
1125     while ((hist = static_cast<TH1*>(next()))) {
1126       TString t(hist->GetTitle());
1127       TString n(hist->GetName());
1128       n.ToLower();
1129       if (t.Contains("mirrored")) continue;
1130       if (n.Contains("syserror")) { sysErrSeen = true; continue; }
1131       if (unique.FindObject(t.Data())) continue;
1132       TObjString* s1 = new TObjString(hist->GetTitle());
1133       s1->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) |
1134                       ((hist->GetMarkerColor() & 0xFFFF) <<  0));
1135       unique.Add(s1);
1136       // l->AddEntry(hist, hist->GetTitle(), "pl");
1137     }
1138     if (mg) {
1139       // If we have other stuff, scan for unique names
1140       TIter nexto(mg->GetListOfGraphs());
1141       TGraph* g = 0;
1142       while ((g = static_cast<TGraph*>(nexto()))) {
1143         TString n(g->GetTitle());
1144         if (n.Contains("mirrored")) continue;
1145         if (unique.FindObject(n.Data())) continue;
1146         TObjString* s2 = new TObjString(n);
1147         s2->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) |
1148                          ((g->GetMarkerColor() & 0xFFFF) <<  0));
1149         unique.Add(s2);
1150         // l->AddEntry(hist, hist->GetTitle(), "pl");
1151       }
1152     }
1153
1154     // Add legend entries for unique items only
1155     TIter nextu(&unique);
1156     TObject* s = 0;
1157     Int_t i = 0;
1158     while ((s = nextu())) {
1159       TLegendEntry* dd = l->AddEntry(Form("data%2d", i++),
1160                                      s->GetName(), "lp");
1161       Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF;
1162       Int_t color = (s->GetUniqueID() >>  0) & 0xFFFF;
1163       dd->SetLineColor(kBlack);
1164       if (HasCent()) dd->SetMarkerColor(kBlack);
1165       else           dd->SetMarkerColor(color);
1166       dd->SetMarkerStyle(style);
1167     }
1168     if (sysErrSeen) {
1169       // Add entry for systematic errors
1170       TLegendEntry* d0 = l->AddEntry("d0", Form("%4.1f%% Systematic error",
1171                                                 100*fFwdSysErr), "f");
1172       d0->SetLineColor(SYSERR_COLOR);
1173       d0->SetMarkerColor(SYSERR_COLOR);
1174       d0->SetFillColor(SYSERR_COLOR);
1175       d0->SetFillStyle(SYSERR_STYLE);
1176       d0->SetMarkerStyle(0);
1177       d0->SetLineWidth(0);
1178       i++;
1179     }
1180     if (nCol == 2 && i % 2 == 1)  {
1181       // To make sure the 'data' and 'mirrored' entries are on a line
1182       // by themselves
1183       TLegendEntry* dd = l->AddEntry("dd", "   ", "");
1184       dd->SetTextSize(0);
1185       dd->SetFillStyle(0);
1186       dd->SetFillColor(0);
1187       dd->SetLineWidth(0);
1188       dd->SetLineColor(0);
1189       dd->SetMarkerSize(0);
1190     }
1191     if (fMirror) {
1192       // Add entry for 'data'
1193       TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp");
1194       d1->SetLineColor(kBlack);
1195       d1->SetMarkerColor(kBlack);
1196       d1->SetMarkerStyle(20);
1197
1198       // Add entry for 'mirrored data'
1199       TLegendEntry* d2 = l->AddEntry("d2", "Mirrored data", "lp");
1200       d2->SetLineColor(kBlack);
1201       d2->SetMarkerColor(kBlack);
1202       d2->SetMarkerStyle(24);
1203     }
1204
1205     l->Draw();
1206   }
1207   //__________________________________________________________________
1208   /**
1209    * Build centrality legend
1210    *
1211    * @param x1      Lower X coordinate in the range [0,1]
1212    * @param y1      Lower Y coordinate in the range [0,1]
1213    * @param x2      Upper X coordinate in the range [0,1]
1214    * @param y2      Upper Y coordinate in the range [0,1]
1215    */
1216   void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
1217   {
1218     if (!HasCent()) return;
1219
1220     TLegend* l = new TLegend(x1,y1,x2,y2);
1221     l->SetNColumns(1);
1222     l->SetFillColor(0);
1223     l->SetFillStyle(0);
1224     l->SetBorderSize(0);
1225     l->SetTextFont(kFont);
1226     l->SetTextColor(kAliceBlue);
1227
1228     Int_t n = fCentAxis->GetNbins();
1229     for (Int_t i = 1; i <= n; i++) {
1230       Double_t low = fCentAxis->GetBinLowEdge(i);
1231       Double_t upp = fCentAxis->GetBinUpEdge(i);
1232       TLegendEntry* e = l->AddEntry(Form("dummy%02d", i),
1233                                     Form("%3d%% - %3d%%",
1234                                          int(low), int(upp)), "pl");
1235       e->SetMarkerColor(GetCentralityColor(i));
1236     }
1237     l->Draw();
1238   }
1239   //__________________________________________________________________
1240   /**
1241    * Plot the results
1242    *
1243    * @param max       Maximum
1244    * @param yd        Bottom position of pad
1245    */
1246   void PlotResults(Double_t max, Double_t yd)
1247   {
1248     // Make a sub-pad for the result itself
1249     TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
1250     p1->SetTopMargin(0.05);
1251     p1->SetBorderSize(0);
1252     p1->SetBorderMode(0);
1253     p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
1254     p1->SetRightMargin(kRightMargin);
1255     if (fShowLeftRight || fShowRatios) p1->SetGridx();
1256     p1->SetTicks(1,1);
1257     p1->SetNumber(1);
1258     p1->Draw();
1259     p1->cd();
1260
1261     Info("PlotResults", "Plotting results with max=%f", max);
1262     fResults->SetMaximum(1.15*max);
1263     fResults->SetMinimum(yd > 0.00001 ? -0.02*max : 0);
1264     // fResults->SetMinimum(yd > 0.00001 ? -0.02*max : 0);
1265
1266     FixAxis(fResults, (1-yd)*(yd > .001 ? 1 : .9 / 1.2),
1267             "#frac{1}{#it{N}}#kern[.1]{#frac{d#it{N}_{ch}}{d#it{#eta}}}");
1268
1269     p1->Clear();
1270     fResults->DrawClone("nostack e1");
1271
1272     fRangeParam->fSlave1Axis = fResults->GetXaxis();
1273     fRangeParam->fSlave1Pad  = p1;
1274
1275     // Draw other data
1276     if (fShowOthers != 0) {
1277       TGraphAsymmErrors* o      = 0;
1278       TIter              nextG(fOthers->GetListOfGraphs());
1279       while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
1280         o->DrawClone("same p");
1281     }
1282
1283     // Make a legend in the main result pad
1284     BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5,
1285                     .35, 1-p1->GetTopMargin()-.01-.1);
1286     Double_t x1 = (HasCent() ? .7 : .15);
1287     Double_t x2 = (HasCent() ? 1-p1->GetRightMargin()-.01: .90);
1288     Double_t y1 = (HasCent() ? .5: p1->GetBottomMargin()+.01);
1289     Double_t y2 = (HasCent() ? 1-p1->GetTopMargin()-.01-.15 : .35);
1290
1291     BuildLegend(fResults, fOthers, x1, y1, x2, y2);
1292
1293     // Put a title on top
1294     fTitle.ReplaceAll("@", " ");
1295     TLatex* tit = new TLatex(0.10, .99, fTitle.Data());
1296     tit->SetNDC();
1297     tit->SetTextFont(kFont);
1298     tit->SetTextAlign(13);
1299     tit->SetTextSize(0.045);
1300     tit->SetTextColor(kAlicePurple);
1301     tit->Draw();
1302
1303     // TColor::GetColor(41,73,156);
1304     Double_t x         = .95;
1305     Double_t y         = .93;
1306     // Put a nice label in the plot
1307     TString     eS;
1308     UShort_t    snn = fSNNString->GetUniqueID();
1309     const char* sys = fSysString->GetTitle();
1310     if (snn == 2750) snn = 2760;
1311     if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
1312     else            eS = Form("%3dGeV", snn);
1313     TLatex* tt = new TLatex(x, y, Form("%s #sqrt{s%s}=%s, %s",
1314                                        sys,
1315                                        (HasCent() ? "_{NN}" : ""),
1316                                        eS.Data(),
1317                                        HasCent() ? "by centrality" :
1318                                        fTrigString->GetTitle()));
1319     tt->SetTextColor(kAliceBlue);
1320     tt->SetNDC();
1321     tt->SetTextFont(kFont);
1322     tt->SetTextAlign(33);
1323     tt->Draw();
1324     y -= tt->GetTextSize() + .01;
1325
1326     // Put number of accepted events on the plot
1327     Int_t nev = 0;
1328     if (fTriggers) nev = fTriggers->GetBinContent(1);
1329     TLatex* et = new TLatex(x, y, Form("%d events", nev));
1330     et->SetTextColor(kAliceBlue);
1331     et->SetNDC();
1332     et->SetTextFont(kFont);
1333     et->SetTextAlign(33);
1334     et->Draw();
1335     y -= et->GetTextSize() + .01;
1336
1337     // Put number of accepted events on the plot
1338     if (fVtxAxis) {
1339       TLatex* vt = new TLatex(x, y, fVtxAxis->GetTitle());
1340       vt->SetNDC();
1341       vt->SetTextFont(kFont);
1342       vt->SetTextAlign(33);
1343       vt->SetTextColor(kAliceBlue);
1344       vt->Draw();
1345       y -= vt->GetTextSize() + .01;
1346     }
1347     // results->Draw("nostack e1 same");
1348
1349     TString corrs;
1350     if (!fEmpirical.IsNull()) corrs.Append("Emperical");
1351     if (!fFinalMC.IsNull())   {
1352       if (!corrs.IsNull()) corrs.Append("+");
1353       corrs.Append("Final MC");
1354     }
1355
1356     if (!corrs.IsNull()) {
1357       corrs.Append(" correction");
1358       if (corrs.Index("+") != kNPOS) corrs.Append("s");
1359       TLatex* em = new TLatex(x, y, corrs);
1360       em->SetNDC();
1361       em->SetTextFont(kFont);
1362       em->SetTextAlign(33);
1363       em->SetTextColor(kAliceBlue);
1364       em->Draw();
1365       y -= em->GetTextSize() + .01;
1366     }
1367
1368     if (fTriggerEff > 0 && fTriggerEff <= 1 && !HasCent()) {
1369       TLatex* ef = new TLatex(x, y, Form("#varepsilon_{%s} = %5.3f",
1370                                          fTrigString->GetTitle(),
1371                                          fTriggerEff));
1372       ef->SetNDC();
1373       ef->SetTextFont(kFont);
1374       ef->SetTextAlign(33);
1375       ef->SetTextColor(kAliceBlue);
1376       ef->Draw();
1377       y -= ef->GetTextSize() + .01;
1378     }
1379
1380     fRangeParam->fSlave1Axis = FindXAxis(p1, fResults->GetName());
1381     fRangeParam->fSlave1Pad  = p1;
1382
1383
1384     // Mark the plot as preliminary
1385     TLatex* pt = new TLatex(.12, .93, "Work in progress");
1386     pt->SetNDC();
1387     pt->SetTextFont(62);
1388     // pt->SetTextSize();
1389     pt->SetTextColor(kAliceRed);
1390     pt->SetTextAlign(13);
1391     pt->Draw();
1392     TDatime now;
1393     TLatex* dt = new TLatex(.12, .88, now.AsSQLString());
1394     dt->SetNDC();
1395     dt->SetTextFont(42);
1396     dt->SetTextSize(0.04);
1397     dt->SetTextColor(kAliceBlue); // kAliceRed);
1398     dt->SetTextAlign(13);
1399     dt->Draw();
1400
1401     TString savPath(gROOT->GetMacroPath());
1402     gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2/scripts",
1403                              gROOT->GetMacroPath()));
1404     // Always recompile
1405     if (!gROOT->GetClass("AliceLogo"))
1406       gROOT->LoadMacro("AliceLogo.C++");
1407     gROOT->SetMacroPath(savPath);
1408
1409     if (gROOT->GetClass("AliceLogo")) {
1410       p1->cd();
1411       p1->Range(0,0,1,1);
1412       gROOT->ProcessLine("AliceLogo* al = new AliceLogo();");
1413       gROOT->ProcessLine("al->Draw(0,.88,.5,.2, 0, 0);");
1414     }
1415     p1->cd();
1416   }
1417   //__________________________________________________________________
1418   /**
1419    * Plot the ratios
1420    *
1421    * @param max     Maximum diviation from 1
1422    * @param y1      Lower y coordinate of pad
1423    * @param y2      Upper y coordinate of pad
1424    */
1425   void PlotRatios(Double_t max, Double_t y1, Double_t y2)
1426   {
1427     if (fShowRatios == 0) return;
1428
1429     bool isBottom = (y1 < 0.0001);
1430     Double_t yd = y2 - y1;
1431     // Make a sub-pad for the result itself
1432     TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
1433     p2->SetTopMargin(0.001);
1434     p2->SetRightMargin(kRightMargin);
1435     p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1436     p2->SetGridx();
1437     p2->SetTicks(1,1);
1438     p2->SetNumber(2);
1439     p2->Draw();
1440     p2->cd();
1441
1442     // Fix up axis
1443     FixAxis(fRatios, yd, "Ratios", 7);
1444
1445     fRatios->SetMaximum(1+TMath::Max(.22,1.05*max));
1446     fRatios->SetMinimum(1-TMath::Max(.32,1.05*max));
1447     p2->Clear();
1448     fRatios->DrawClone("nostack e1");
1449
1450
1451     // Make a legend
1452     BuildLegend(fRatios, 0, .15,p2->GetBottomMargin()+.01,.9,
1453                 isBottom ? .6 : .4, 2);
1454 #if 0
1455     TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
1456                                   isBottom ? .6 : .4);
1457     l2->SetNColumns(2);
1458     l2->SetFillColor(0);
1459     l2->SetFillStyle(0);
1460     l2->SetBorderSize(0);
1461     l2->SetTextFont(kFont);
1462 #endif
1463     // Make a nice band from 0.9 to 1.1
1464     TGraphErrors* band = new TGraphErrors(2);
1465     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1466     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
1467     band->SetPointError(0, 0, .1);
1468     band->SetPointError(1, 0, .1);
1469     band->SetFillColor(kYellow+2);
1470     band->SetFillStyle(3002);
1471     band->SetLineStyle(2);
1472     band->SetLineWidth(1);
1473     band->Draw("3 same");
1474     band->DrawClone("X L same");
1475
1476     // Replot the ratios on top
1477     fRatios->DrawClone("nostack e1 same");
1478
1479     if (fAddExec) {
1480       if (isBottom) {
1481         fRangeParam->fMasterAxis = FindXAxis(p2, fRatios->GetName());
1482         p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1483                                   fRangeParam));
1484       }
1485       else {
1486         fRangeParam->fSlave2Axis = FindXAxis(p2, fRatios->GetName());
1487         fRangeParam->fSlave2Pad  = p2;
1488       }
1489     }
1490   }
1491   //__________________________________________________________________
1492   /**
1493    * Plot the asymmetries
1494    *
1495    * @param max     Maximum diviation from 1
1496    * @param y1      Lower y coordinate of pad
1497    * @param y2      Upper y coordinate of pad
1498    */
1499   void PlotLeftRight(Double_t max, Double_t y1, Double_t y2)
1500   {
1501     if (!fShowLeftRight) return;
1502
1503     bool isBottom = (y1 < 0.0001);
1504     Double_t yd = y2 - y1;
1505     // Make a sub-pad for the result itself
1506     TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
1507     p3->SetTopMargin(0.001);
1508     p3->SetRightMargin(kRightMargin);
1509     p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
1510     p3->SetGridx();
1511     p3->SetTicks(1,1);
1512     p3->SetNumber(2);
1513     p3->Draw();
1514     p3->cd();
1515
1516     // Fix up axis
1517     FixAxis(fLeftRight, yd, "Right/Left", 4);
1518
1519     fLeftRight->SetMaximum(1+TMath::Max(.12,1.05*max));
1520     fLeftRight->SetMinimum(1-TMath::Max(.15,1.05*max));
1521     p3->Clear();
1522     fLeftRight->DrawClone("nostack e1");
1523
1524
1525     // Make a legend
1526     Double_t xx1 = (HasCent() ? .7                           : .15);
1527     Double_t xx2 = (HasCent() ? 1-p3->GetRightMargin()-.01   : .90);
1528     Double_t yy1 = p3->GetBottomMargin()+.01;
1529     Double_t yy2 = (HasCent() ? 1-p3->GetTopMargin()-.01-.15 : .5);
1530     BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2);
1531     // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
1532     // l2->SetNColumns(2);
1533     // l2->SetFillColor(0);
1534     // l2->SetFillStyle(0);
1535     // l2->SetBorderSize(0);
1536     // l2->SetTextFont(132);
1537
1538     // Make a nice band from 0.9 to 1.1
1539     TGraphErrors* band = new TGraphErrors(2);
1540     band->SetPoint(0, fResults->GetXaxis()->GetXmin(), 1);
1541     band->SetPoint(1, fResults->GetXaxis()->GetXmax(), 1);
1542     band->SetPointError(0, 0, .05);
1543     band->SetPointError(1, 0, .05);
1544     band->SetFillColor(kYellow+2);
1545     band->SetFillStyle(3002);
1546     band->SetLineStyle(2);
1547     band->SetLineWidth(1);
1548     band->Draw("3 same");
1549     band->DrawClone("X L same");
1550
1551     fLeftRight->DrawClone("nostack e1 same");
1552     if (isBottom) {
1553       fRangeParam->fMasterAxis = FindXAxis(p3, fLeftRight->GetName());
1554       p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
1555                                 fRangeParam));
1556     }
1557     else {
1558       fRangeParam->fSlave2Axis = FindXAxis(p3, fLeftRight->GetName());
1559       fRangeParam->fSlave2Pad  = p3;
1560     }
1561   }
1562   /** @} */
1563   //==================================================================
1564   /**
1565    * @{
1566    * @name Data utility functions
1567    */
1568   /**
1569    * Get a result from the passed list
1570    *
1571    * @param list List to search
1572    * @param name Object name to search for
1573    *
1574    * @return Histogram
1575    */
1576   TH1* FetchResult(const TList* list, const char* name) const
1577   {
1578     if (!list) return 0;
1579
1580     TH1* ret = static_cast<TH1*>(list->FindObject(name));
1581 #if 0
1582     if (!ret) {
1583       // all->ls();
1584       Warning("GetResult", "Histogram %s not found", name);
1585     }
1586 #endif
1587
1588     return ret;
1589   }
1590   //__________________________________________________________________
1591   /**
1592    * Add a histogram to the stack after possibly rebinning it
1593    *
1594    * @param stack   Stack to add to
1595    * @param hist    histogram
1596    * @param option  Draw options
1597    *
1598    * @return Maximum of histogram
1599    */
1600   Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option="") const
1601   {
1602     // Check if we have input
1603     if (!hist) return 0;
1604
1605     // Rebin if needed
1606     Rebin(hist);
1607
1608     stack->Add(hist, option);
1609     return hist->GetMaximum();
1610   }
1611   //__________________________________________________________________
1612   /**
1613    * Add a histogram to the stack after possibly rebinning it
1614    *
1615    * @param stack   Stack to add to
1616    * @param hist    histogram
1617    * @param option  Draw options
1618    * @param sym     On return, the data symmetriced (added to stack)
1619    *
1620    * @return Maximum of histogram
1621    */
1622   Double_t AddHistogram(THStack* stack, TH1* hist, TH1*& sym,
1623                         Option_t* option="") const
1624   {
1625     // Check if we have input
1626     if (!hist) return 0;
1627
1628     // Rebin if needed
1629     Rebin(hist);
1630     stack->Add(hist, option);
1631
1632     // Now symmetrice the histogram
1633     if (fMirror) {
1634       sym = Symmetrice(hist);
1635       stack->Add(sym, option);
1636     }
1637
1638     return hist->GetMaximum();
1639   }
1640
1641   //__________________________________________________________________
1642   /**
1643    * Rebin a histogram
1644    *
1645    * @param h     Histogram to rebin
1646    */
1647   virtual void Rebin(TH1* h) const
1648   {
1649     if (fRebin <= 1) return;
1650
1651     Int_t nBins = h->GetNbinsX();
1652     if(nBins % fRebin != 0) {
1653       Warning("Rebin", "Rebin factor %d is not a devisor of current number "
1654               "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
1655       return;
1656     }
1657
1658     // Make a copy
1659     TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
1660     tmp->Rebin(fRebin);
1661     tmp->SetDirectory(0);
1662     tmp->Reset();
1663     // The new number of bins
1664     Int_t nBinsNew = nBins / fRebin;
1665     for(Int_t i = 1;i<= nBinsNew; i++) {
1666       Double_t content = 0;
1667       Double_t sumw    = 0;
1668       Double_t wsum    = 0;
1669       Int_t    nbins   = 0;
1670       for(Int_t j = 1; j<=fRebin;j++) {
1671         Int_t    bin = (i-1)*fRebin + j;
1672         Double_t c   =  h->GetBinContent(bin);
1673
1674         if (c <= 0) continue;
1675
1676         if (fCutEdges) {
1677           if (h->GetBinContent(bin+1)<=0 ||
1678               h->GetBinContent(bin-1)<=0) {
1679             Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
1680                     bin, c, h->GetName(),
1681                     bin+1, h->GetBinContent(bin+1),
1682                     bin-1, h->GetBinContent(bin-1));
1683             continue;
1684           }
1685         }
1686         Double_t e =  h->GetBinError(bin);
1687         Double_t w =  1 / (e*e); // 1/c/c
1688         content    += c;
1689         sumw       += w;
1690         wsum       += w * c;
1691         nbins++;
1692       }
1693
1694       if(content > 0 && nbins > 1 ) {
1695         tmp->SetBinContent(i, wsum / sumw);
1696         tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1697       }
1698     }
1699
1700     // Finally, rebin the histogram, and set new content
1701     h->Rebin(fRebin);
1702     h->Reset();
1703     for(Int_t i = 1; i<= nBinsNew; i++) {
1704       h->SetBinContent(i,tmp->GetBinContent(i));
1705       h->SetBinError(i,  tmp->GetBinError(i));
1706     }
1707
1708     delete tmp;
1709   }
1710   //__________________________________________________________________
1711   /**
1712    * Make an extension of @a h to make it symmetric about 0
1713    *
1714    * @param h Histogram to symmertrice
1715    *
1716    * @return Symmetric extension of @a h
1717    */
1718   TH1* Symmetrice(const TH1* h) const
1719   {
1720     Int_t nBins = h->GetNbinsX();
1721     TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1722     s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1723     s->Reset();
1724     s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1725     s->SetMarkerColor(h->GetMarkerColor());
1726     s->SetMarkerSize(h->GetMarkerSize());
1727     s->SetMarkerStyle(h->GetMarkerStyle()+4);
1728     s->SetFillColor(h->GetFillColor());
1729     s->SetFillStyle(h->GetFillStyle());
1730     s->SetDirectory(0);
1731
1732     // Find the first and last bin with data
1733     Int_t first = nBins+1;
1734     Int_t last  = 0;
1735     for (Int_t i = 1; i <= nBins; i++) {
1736       if (h->GetBinContent(i) <= 0) continue;
1737       first = TMath::Min(first, i);
1738       last  = TMath::Max(last,  i);
1739     }
1740
1741     Double_t xfirst = h->GetBinCenter(first-1);
1742     Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
1743     Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
1744     for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1745       s->SetBinContent(j, h->GetBinContent(i));
1746       s->SetBinError(j, h->GetBinError(i));
1747     }
1748     // Fill in overlap bin
1749     s->SetBinContent(l2+1, h->GetBinContent(first));
1750     s->SetBinError(l2+1, h->GetBinError(first));
1751     return s;
1752   }
1753   //__________________________________________________________________
1754   /**
1755    * Calculate the left-right asymmetry of input histogram
1756    *
1757    * @param h   Input histogram
1758    * @param max On return, the maximum distance from 1 of the histogram
1759    *
1760    * @return Asymmetry
1761    */
1762   TH1* Asymmetry(TH1* h, Double_t& max)
1763   {
1764     if (!h) return 0;
1765
1766     TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
1767     // Int_t    oBins = h->GetNbinsX();
1768     // Double_t high  = h->GetXaxis()->GetXmax();
1769     // Double_t low   = h->GetXaxis()->GetXmin();
1770     // Double_t dBin  = (high - low) / oBins;
1771     // Int_t    tBins = Int_t(2*high/dBin+.5);
1772     // ret->SetBins(tBins, -high, high);
1773     ret->SetDirectory(0);
1774     ret->Reset();
1775     ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
1776     ret->SetYTitle("Right/Left");
1777     Int_t nBins = h->GetNbinsX();
1778     for (Int_t i = 1; i <= nBins; i++) {
1779       Double_t x = h->GetBinCenter(i);
1780       if (x > 0) break;
1781
1782       Double_t c1 = h->GetBinContent(i);
1783       Double_t e1 = h->GetBinError(i);
1784       if (c1 <= 0) continue;
1785
1786       Int_t    j  = h->FindBin(-x);
1787       if (j <= 0 || j > nBins) continue;
1788
1789       Double_t c2 = h->GetBinContent(j);
1790       Double_t e2 = h->GetBinError(j);
1791
1792       Double_t c12 = c1*c1;
1793       Double_t e   = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
1794
1795       Int_t    k   = ret->FindBin(x);
1796       ret->SetBinContent(k, c2/c1);
1797       ret->SetBinError(k, e);
1798     }
1799     max = TMath::Max(max, RatioMax(ret));
1800
1801     return ret;
1802   }
1803   //__________________________________________________________________
1804   /**
1805    * Transform a graph into a histogram
1806    *
1807    * @param g Graph
1808    *
1809    * @return Histogram
1810    */
1811   TH1* Graph2Hist(const TGraphAsymmErrors* g) const
1812   {
1813     Int_t    nBins = g->GetN();
1814     TArrayF  bins(nBins+1);
1815     Double_t dx = 0;
1816     for (Int_t i = 0; i < nBins; i++) {
1817       Double_t x   = g->GetX()[i];
1818       Double_t exl = g->GetEXlow()[i];
1819       Double_t exh = g->GetEXhigh()[i];
1820       bins.fArray[i]   = x-exl;
1821       bins.fArray[i+1] = x+exh;
1822       Double_t dxi = exh+exl;
1823       if (i == 0) dx  = dxi;
1824       else if (dxi != dx) dx = 0;
1825     }
1826     TString name(g->GetName());
1827     TString title(g->GetTitle());
1828     TH1D* h = 0;
1829     if (dx != 0) {
1830       h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
1831     }
1832     else {
1833       h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
1834     }
1835     h->SetMarkerStyle(g->GetMarkerStyle());
1836     h->SetMarkerColor(g->GetMarkerColor());
1837     h->SetMarkerSize(g->GetMarkerSize());
1838
1839     return h;
1840   }
1841   /* @} */
1842   //==================================================================
1843   /**
1844    * @{
1845    * @name Ratio utility functions
1846    */
1847   /**
1848    * Get the maximum diviation from 1 in the passed ratio
1849    *
1850    * @param h Ratio histogram
1851    *
1852    * @return Max diviation from 1
1853    */
1854   Double_t RatioMax(TH1* h) const
1855   {
1856     Int_t    nBins = h->GetNbinsX();
1857     Double_t ret   = 0;
1858     for (Int_t i = 1; i <= nBins; i++) {
1859       Double_t c = h->GetBinContent(i);
1860       if (c == 0) continue;
1861       Double_t e = h->GetBinError(i);
1862       Double_t d = TMath::Abs(1-c-e);
1863       ret        = TMath::Max(d, ret);
1864     }
1865     return ret;
1866   }
1867   //__________________________________________________________________
1868   /**
1869    * Make the ratio of h1 to h2
1870    *
1871    * @param o1  First object (numerator)
1872    * @param o2  Second object (denominator)
1873    * @param max Maximum diviation from 1
1874    *
1875    * @return o1 / o2
1876    */
1877   TH1* Ratio(const TObject* o1, const TObject* o2, Double_t& max) const
1878   {
1879     if (!o1 || !o2) return 0;
1880
1881     TH1*        r  = 0;
1882     const TAttMarker* m1 = 0;
1883     const TAttMarker* m2 = 0;
1884     const TH1* h1 = dynamic_cast<const TH1*>(o1);
1885     if (h1) {
1886       m1 = h1;
1887       const TH1* h2 = dynamic_cast<const TH1*>(o2);
1888       if (h2) {
1889         m2 = h2;
1890         r = RatioHH(h1,h2);
1891       }
1892       else {
1893         const TGraph* g2 = dynamic_cast<const TGraph*>(o2);
1894         if (g2) {
1895           m2 = g2;
1896           r = RatioHG(h1,g2);
1897         }
1898       }
1899     }
1900     else {
1901       const TGraphAsymmErrors* g1 = dynamic_cast<const TGraphAsymmErrors*>(o1);
1902       if (g1) {
1903         m1 = g1;
1904         const TGraphAsymmErrors* g2 =
1905           dynamic_cast<const TGraphAsymmErrors*>(o2);
1906         if (g2) {
1907           m2 = g2;
1908           r = RatioGG(g1, g2);
1909         }
1910       }
1911     }
1912     if (!r) {
1913       Warning("Ratio", "Don't know how to divide a %s (%s) with a %s (%s)",
1914               o1->ClassName(), o1->GetName(), o2->ClassName(), o2->GetName());
1915       return 0;
1916     }
1917     // Check that the histogram isn't empty
1918     if (r->GetEntries() <= 0) {
1919       delete r;
1920       r = 0;
1921     }
1922     if (r) {
1923       r->SetMarkerStyle(m2->GetMarkerStyle());
1924       r->SetMarkerColor(m1->GetMarkerColor());
1925       if (TString(o2->GetName()).Contains("truth", TString::kIgnoreCase))
1926         r->SetMarkerStyle(m1->GetMarkerStyle());
1927       r->SetMarkerSize(0.9*m1->GetMarkerSize());
1928       r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName()));
1929       r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle()));
1930       r->SetDirectory(0);
1931       max = TMath::Max(RatioMax(r), max);
1932     }
1933
1934     return r;
1935   }
1936   //__________________________________________________________________
1937   /**
1938    * Compute the ratio of @a h to @a g.  @a g is evaluated at the bin
1939    * centers of @a h
1940    *
1941    * @param h  Numerator
1942    * @param g  Divisor
1943    *
1944    * @return h/g
1945    */
1946   TH1* RatioHG(const TH1* h, const TGraph* g) const
1947   {
1948     if (!h || !g) return 0;
1949
1950     TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
1951     ret->Reset();
1952     Double_t xlow  = g->GetX()[0];
1953     Double_t xhigh = g->GetX()[g->GetN()-1];
1954     if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
1955
1956     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
1957       Double_t c = h->GetBinContent(i);
1958       if (c <= 0) continue;
1959
1960       Double_t x = h->GetBinCenter(i);
1961       if (x < xlow || x > xhigh) continue;
1962
1963       Double_t f = g->Eval(x);
1964       if (f <= 0) continue;
1965
1966       ret->SetBinContent(i, c / f);
1967       ret->SetBinError(i, h->GetBinError(i) / f);
1968     }
1969     return ret;
1970   }
1971   //__________________________________________________________________
1972   /**
1973    * Make the ratio of h1 to h2
1974    *
1975    * @param h1 First histogram (numerator)
1976    * @param h2 Second histogram (denominator)
1977    *
1978    * @return h1 / h2
1979    */
1980   TH1* RatioHH(const TH1* h1, const TH1* h2) const
1981   {
1982     if (!h1 || !h2) return 0;
1983     TH1* t1 = static_cast<TH1*>(h1->Clone("tmp"));
1984     t1->Divide(h2);
1985     return t1;
1986   }
1987   //__________________________________________________________________
1988   /**
1989    * Calculate the ratio of two graphs - g1 / g2
1990    *
1991    * @param g1 Numerator
1992    * @param g2 Denominator
1993    *
1994    * @return g1 / g2 in a histogram
1995    */
1996   TH1* RatioGG(const TGraphAsymmErrors* g1,
1997                const TGraphAsymmErrors* g2) const
1998   {
1999     Int_t    nBins = g1->GetN();
2000     TArrayF  bins(nBins+1);
2001     Double_t dx   = 0;
2002     for (Int_t i = 0; i < nBins; i++) {
2003       Double_t x   = g1->GetX()[i];
2004       Double_t exl = g1->GetEXlow()[i];
2005       Double_t exh = g1->GetEXhigh()[i];
2006       bins.fArray[i]   = x-exl;
2007       bins.fArray[i+1] = x+exh;
2008       Double_t dxi = exh+exl;
2009       if (i == 0) dx  = dxi;
2010       else if (dxi != dx) dx = 0;
2011     }
2012     TH1* h = 0;
2013     if (dx != 0) {
2014       h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]);
2015     }
2016     else {
2017       h = new TH1F("tmp", "tmp", nBins, bins.fArray);
2018     }
2019
2020     Double_t low  = g2->GetX()[0];
2021     Double_t high = g2->GetX()[g2->GetN()-1];
2022     if (low > high) { Double_t t = low; low = high; high = t; }
2023     for (Int_t i = 0; i < nBins; i++) {
2024       Double_t x  = g1->GetX()[i];
2025       if (x < low-dx || x > high+dx) continue;
2026       Double_t c1 = g1->GetY()[i];
2027       Double_t e1 = g1->GetErrorY(i);
2028       Double_t c2 = g2->Eval(x);
2029
2030       h->SetBinContent(i+1, c1 / c2);
2031       h->SetBinError(i+1, e1 / c2);
2032     }
2033     return h;
2034   }
2035   /* @} */
2036   //==================================================================
2037   /**
2038    * @{
2039    * @name Graphics utility functions
2040    */
2041   /**
2042    * Find an X axis in a pad
2043    *
2044    * @param p     Pad
2045    * @param name  Histogram to find axis for
2046    *
2047    * @return Found axis or null
2048    */
2049   TAxis* FindXAxis(TVirtualPad* p, const char* name)
2050   {
2051     TObject* o = p->GetListOfPrimitives()->FindObject(name);
2052     if (!o) {
2053       Warning("FindXAxis", "%s not found in pad", name);
2054       return 0;
2055     }
2056     THStack* stack = dynamic_cast<THStack*>(o);
2057     if (!stack) {
2058       Warning("FindXAxis", "%s is not a THStack", name);
2059       return 0;
2060     }
2061     if (!stack->GetHistogram()) {
2062       Warning("FindXAxis", "%s has no histogram", name);
2063       return 0;
2064     }
2065     TAxis* ret = stack->GetHistogram()->GetXaxis();
2066     return ret;
2067   }
2068
2069   //__________________________________________________________________
2070   /**
2071    * Fix the apperance of the axis in a stack
2072    *
2073    * @param stack  stack of histogram
2074    * @param yd     How the canvas is cut
2075    * @param ytitle Y axis title
2076    * @param ynDiv  Divisions on Y axis
2077    * @param force  Whether to draw the stack first or not
2078    */
2079   void FixAxis(THStack* stack, Double_t yd, const char* ytitle,
2080                Int_t ynDiv=210, Bool_t force=true)
2081   {
2082     if (!stack) {
2083       Warning("FixAxis", "No stack passed for %s!", ytitle);
2084       return;
2085     }
2086     if (force) stack->Draw("nostack e1");
2087
2088     TH1* h = stack->GetHistogram();
2089     if (!h) {
2090       Warning("FixAxis", "Stack %s has no histogram", stack->GetName());
2091       return;
2092     }
2093     Double_t s = 1/yd/1.2;
2094     // Info("FixAxis", "for %s, s=1/%f=%f", stack->GetName(), yd, s);
2095
2096     h->SetXTitle("#it{#eta}");
2097     h->SetYTitle(ytitle);
2098     TAxis* xa = h->GetXaxis();
2099     TAxis* ya = h->GetYaxis();
2100
2101     // Int_t   npixels = 20
2102     // Float_t dy = gPad->AbsPixeltoY(0) - gPad->AbsPixeltoY(npixels);
2103     // Float_t ts = dy/(gPad->GetY2() - gPad->GetY1());
2104
2105     if (xa) {
2106       // xa->SetTitle(h->GetXTitle());
2107       // xa->SetTicks("+-");
2108       xa->SetTitleSize(s*xa->GetTitleSize());
2109       xa->SetLabelSize(s*xa->GetLabelSize());
2110       xa->SetTickLength(s*xa->GetTickLength());
2111       // xa->SetTitleOffset(xa->GetTitleOffset()/s);
2112
2113       if (stack != fResults) {
2114         TAxis* rxa = fResults->GetXaxis();
2115         xa->Set(rxa->GetNbins(), rxa->GetXmin(), rxa->GetXmax());
2116       }
2117     }
2118     if (ya) {
2119       // ya->SetTitle(h->GetYTitle());
2120       ya->SetDecimals();
2121       // ya->SetTicks("+-");
2122       ya->SetNdivisions(ynDiv);
2123       ya->SetTitleSize(s*ya->GetTitleSize());
2124       ya->SetTitleOffset(1.15*ya->GetTitleOffset()/s);
2125       ya->SetLabelSize(s*ya->GetLabelSize());
2126     }
2127   }
2128   //__________________________________________________________________
2129   /**
2130    * Merge two histograms into one
2131    *
2132    * @param cen    Central part
2133    * @param fwd    Forward part
2134    * @param xlow   On return, lower eta bound
2135    * @param xhigh  On return, upper eta bound
2136    *
2137    * @return Newly allocated histogram or null
2138    */
2139   TH1*
2140   Merge(const TH1* cen, const TH1* fwd, Double_t& xlow, Double_t& xhigh)
2141   {
2142     TH1* tmp = static_cast<TH1*>(fwd->Clone("tmp"));
2143     TString name(fwd->GetName());
2144     name.ReplaceAll("Forward", "Merged");
2145     tmp->SetName(name);
2146
2147     // tmp->SetMarkerStyle(28);
2148     // tmp->SetMarkerColor(kBlack);
2149     tmp->SetDirectory(0);
2150     xlow  = 100;
2151     xhigh = -100;
2152     for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
2153       Double_t cc = cen->GetBinContent(i);
2154       Double_t cf = fwd->GetBinContent(i);
2155       Double_t ec = cen->GetBinError(i);
2156       Double_t ef = fwd->GetBinError(i);
2157       Double_t nc = cf;
2158       Double_t ne = ef;
2159       if (cc < 0.001 && cf < 0.01) continue;
2160       xlow  = TMath::Min(tmp->GetXaxis()->GetBinLowEdge(i),xlow);
2161       xhigh = TMath::Max(tmp->GetXaxis()->GetBinUpEdge(i),xhigh);
2162       if (cc > 0.001) {
2163         nc = cc;
2164         ne = ec;
2165         if (cf > 0.001) {
2166           nc  = (cf + cc) / 2;
2167           ne  = TMath::Sqrt(ec*ec + ef*ef);
2168         }
2169       }
2170       tmp->SetBinContent(i, nc);
2171       tmp->SetBinError(i, ne);
2172     }
2173     return tmp;
2174   }
2175   //____________________________________________________________________
2176   /**
2177    * Fit  @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$ to histogram data
2178    *
2179    * @param tmp    Histogram
2180    * @param xlow   Lower x bound
2181    * @param xhigh  Upper x bound
2182    *
2183    * @return Fitted function
2184    */
2185   TF1*
2186   FitMerged(TH1* tmp, Double_t xlow, Double_t xhigh)
2187   {
2188     TF1* tmpf  = new TF1("tmpf", "gaus", xlow, xhigh);
2189     tmp->Fit(tmpf, "NQ", "");
2190     tmp->SetDirectory(0);
2191
2192     TF1* fit = new TF1("f", myFunc, xlow, xhigh, 4);
2193     fit->SetParNames("a_{1}", "a_{2}", "#sigma_{1}", "#sigma_{2}");
2194     fit->SetParameters(tmpf->GetParameter(0),
2195                        .2,
2196                        tmpf->GetParameter(2),
2197                        tmpf->GetParameter(2)/4);
2198     fit->SetParLimits(3, 0, 100);
2199     fit->SetParLimits(4, 0, 100);
2200     tmp->Fit(fit,"0W","");
2201
2202     delete tmpf;
2203     return fit;
2204   }
2205   //____________________________________________________________________
2206   /**
2207    * Make band of systematic errors
2208    *
2209    * @param tmp Histogram
2210    * @param cen Central
2211    * @param fwd Forward
2212    * @param fit Fit
2213    */
2214   void
2215   MakeSysError(TH1* tmp, TH1* cen, TH1* fwd, TF1* fit)
2216   {
2217     for (Int_t i = 1; i <= tmp->GetNbinsX(); i++) {
2218       Double_t tc = tmp->GetBinContent(i);
2219       if (tc < 0.01) continue;
2220       Double_t fc = fwd->GetBinContent(i);
2221       Double_t cc = cen->GetBinContent(i);
2222       Double_t sysErr = fFwdSysErr;
2223       if (cc > .01 && fc > 0.01)
2224         sysErr = (fFwdSysErr+fCenSysErr) / 2;
2225       else if (cc > .01)
2226         sysErr = fCenSysErr;
2227       Double_t x = tmp->GetXaxis()->GetBinCenter(i);
2228       Double_t y = fit->Eval(x);
2229       tmp->SetBinContent(i, y);
2230       tmp->SetBinError(i,sysErr*y);
2231     }
2232     TString name(tmp->GetName());
2233     name.ReplaceAll("Merged", "SysError");
2234     tmp->SetName(name);
2235     tmp->SetMarkerColor(SYSERR_COLOR);
2236     tmp->SetLineColor(SYSERR_COLOR);
2237     tmp->SetFillColor(SYSERR_COLOR);
2238     tmp->SetFillStyle(SYSERR_STYLE);
2239     tmp->SetMarkerStyle(0);
2240     tmp->SetLineWidth(0);
2241   }
2242   void CorrectForward(TH1* h) const
2243   {
2244     if (!fRemoveOuters) return;
2245
2246     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
2247       Double_t eta = h->GetBinCenter(i);
2248       if (TMath::Abs(eta) < 2.3) {
2249         h->SetBinContent(i, 0);
2250         h->SetBinError(i, 0);
2251       }
2252     }
2253   }
2254   void CorrectCentral(TH1* h) const
2255   {
2256     if (fClusterScale.IsNull()) return;
2257     TString t(h->GetTitle());
2258     Info("CorrectCentral", "Replacing Central with Tracklets in %s", t.Data());
2259     t.ReplaceAll("Central", "Tracklets");
2260     h->SetTitle(t);
2261
2262     TF1* cf = new TF1("clusterScale", fClusterScale);
2263     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
2264       Double_t eta = h->GetBinCenter(i);
2265       Double_t f   = cf->Eval(eta);
2266       Double_t c   = h->GetBinContent(i);
2267       if (f < .1) f = 1;
2268       h->SetBinContent(i, c / f);
2269     }
2270     delete cf;
2271   }
2272   //____________________________________________________________________
2273   void Export(const char* basename)
2274   {
2275     TString bname(basename);
2276     bname.ReplaceAll(" ", "_");
2277     bname.ReplaceAll("-", "_");
2278     TString fname(Form("%s.C", bname.Data()));
2279
2280     std::ofstream outf(fname.Data());
2281     if (!outf) {
2282       Error("Export", "Failed to open output file %s", fname.Data());
2283       return;
2284     }
2285     outf << "// Create by dNdetaDrawer\n"
2286          << "void " << bname << "(THStack* stack, TLegend* l, Int_t m)\n"
2287          << "{"
2288          << "   Int_t ma[] = { 24, 25, 26, 32,\n"
2289          << "                  20, 21, 22, 33,\n"
2290          << "                  34, 30, 29, 0, \n"
2291          << "                  23, 27 };\n"
2292          << "   Int_t mm = ((m < 20 || m > 34) ? 0 : ma[m-20]);\n\n";
2293     TList* hists = fResults->GetHists();
2294     TIter  next(hists);
2295     TH1*   hist = 0;
2296     while ((hist = static_cast<TH1*>(next()))) {
2297       TString hname = hist->GetName();
2298       hname.Append(Form("_%04x", (gRandom->Integer(0xffff) & 0xffff)));
2299       hist->SetName(hname);
2300       hist->GetListOfFunctions()->Clear();
2301       hist->SavePrimitive(outf, "nodraw");
2302       bool mirror = hname.Contains("mirror");
2303       bool syserr = hname.Contains("SysError");
2304       if (!syserr)
2305         outf << "   " << hname << "->SetMarkerStyle("
2306              << (mirror ? "mm" : "m") << ");\n";
2307       else
2308         outf << "   " << hname << "->SetMarkerStyle(1);\n";
2309       outf << "   stack->Add(" << hname
2310            << (syserr ? ",\"e5\"" : "") << ");\n\n";
2311     }
2312     UShort_t    snn = fSNNString->GetUniqueID();
2313     // const char* sys = fSysString->GetTitle();
2314     TString eS;
2315     if      (snn == 2750)     snn = 2760;
2316     if      (snn < 1000)      eS = Form("%3dGeV", snn);
2317     else if (snn % 1000 == 0) eS = Form("%dTeV", snn/1000);
2318     else                      eS = Form("%4.2fTeV", float(snn)/1000);
2319     outf << "  if (l) {\n"
2320          << "    TLegendEntry* e = l->AddEntry(\"\",\"" << eS << "\",\"pl\");\n"
2321          << "    e->SetMarkerStyle(m);\n"
2322          << "    e->SetMarkerColor(kBlack);\n"
2323          << "  }\n"
2324          << "}\n" << std::endl;
2325   }
2326   /* @} */
2327   /**
2328    * Check if we have centrality dependent information, and we're not
2329    * forcing to use minimum bias
2330    *
2331    *
2332    * @return True if we should do centrality dependent ploting
2333    */
2334   Bool_t HasCent() const { return fCentAxis && !fForceMB; }
2335
2336
2337
2338   //__________________________________________________________________
2339   /**
2340    * @{
2341    * @name Options
2342    */
2343   Bool_t       fShowRatios;   // Show ratios
2344   Bool_t       fShowLeftRight;// Show asymmetry
2345   Bool_t       fShowRings;    // Show rings too
2346   Bool_t       fExport;       // Export results to file
2347   Bool_t       fCutEdges;     // Whether to cut edges
2348   Bool_t       fRemoveOuters; // Whether to remove outers
2349   UShort_t     fShowOthers;   // Show other data
2350   Bool_t       fMirror;       // Whether to mirror
2351   Bool_t       fForceMB;      // Force min-bias
2352   Bool_t       fAddExec;      // Add code to do combined zooms
2353   Bool_t       fSaveAsPDF;    // Flag whether to save as PDF
2354   Bool_t       fSaveAsPNG;    // Flag whether to save as PNG
2355   Bool_t       fSaveAsROOT;   // Flag whether to save as ROOT file
2356   Bool_t       fSaveAsScript; // Flag whether to save as ROOT script
2357   /* @} */
2358   /**
2359    * @{
2360    * @name Settings
2361    */
2362   UShort_t     fRebin;        // Rebinning factor
2363   Double_t     fFwdSysErr;    // Systematic error in forward range
2364   Double_t     fCenSysErr;    // Systematic error in central range
2365   TString      fTitle;        // Title on plot
2366   TString      fBase;         // Base name of output
2367   TString      fClusterScale; // Scaling of clusters to tracklets
2368   TString      fFinalMC;      // Final MC correction file name
2369   TString      fEmpirical;    // Empirical correction file name
2370   /* @} */
2371   /**
2372    * @{
2373    * @name Read (or set) information
2374    */
2375   TNamed*      fTrigString;   // Trigger string (read, or set)
2376   TNamed*      fNormString;   // Normalisation string (read, or set)
2377   TNamed*      fSNNString;    // Energy string (read, or set)
2378   TNamed*      fSysString;    // Collision system string (read or set)
2379   TAxis*       fVtxAxis;      // Vertex cuts (read or set)
2380   TAxis*       fCentAxis;     // Centrality axis
2381   Float_t      fTriggerEff;   // Trigger efficiency
2382   Bool_t       fExtTriggerEff;// True if read externally
2383   /* @} */
2384   /**
2385    * @{
2386    * @name Resulting plots
2387    */
2388   THStack*     fResults;      // Stack of results
2389   THStack*     fRatios;       // Stack of ratios
2390   THStack*     fLeftRight;    // Left-right asymmetry
2391   TMultiGraph* fOthers;       // Older data
2392   TH1*         fTriggers;     // Number of triggers
2393   TH1*         fTruth;        // Pointer to truth
2394   /* @} */
2395   RangeParam*  fRangeParam;   // Parameter object for range zoom
2396   Bool_t       fOld;
2397
2398   static const Float_t kRightMargin;
2399   static const Int_t   kFont;
2400   static const Int_t   kAliceBlue;
2401   static const Int_t   kAliceRed;
2402   static const Int_t   kAlicePurple;
2403   static const Int_t   kAliceYellow;
2404 };
2405
2406 const Float_t dNdetaDrawer::kRightMargin = 0.02;
2407 const Int_t   dNdetaDrawer::kFont        = 42; // 132 for serif
2408 const Int_t   dNdetaDrawer::kAliceBlue   = TColor::GetColor(40,   58, 68);
2409 const Int_t   dNdetaDrawer::kAliceRed    = TColor::GetColor(226,   0, 26);
2410 const Int_t   dNdetaDrawer::kAlicePurple = TColor::GetColor(202,  71, 67);
2411 const Int_t   dNdetaDrawer::kAliceYellow = TColor::GetColor(238, 125, 17);
2412
2413
2414 //____________________________________________________________________
2415 /**
2416  * Function to calculate
2417  * @f[
2418  *  g(x;A_1,A_2,\sigma_1,\sigma_2) =
2419  *       A_1\left(\frac{1}{2\pi\sigma_1}e^{(x/\sigma_1)^2} -
2420  *           A_2\frac{1}{2\pi\sigma_2}e^{(x/\sigma_2)^2}\right)
2421  * @f]
2422  *
2423  * @param xp Pointer to x array
2424  * @param pp Pointer to parameter array
2425  *
2426  * @return @f$g(x;A_1,A_2,\sigma_1,\sigma_2)@f$
2427  */
2428 Double_t myFunc(Double_t* xp, Double_t* pp)
2429 {
2430   Double_t x  = xp[0];
2431   Double_t a1 = pp[0];
2432   Double_t a2 = pp[1];
2433   Double_t s1 = pp[2];
2434   Double_t s2 = pp[3];
2435   return a1*(TMath::Gaus(x, 0, s1) - a2 * TMath::Gaus(x, 0, s2));
2436 }
2437
2438 //=== Stuff for auto zooming =========================================
2439 /**
2440  * Update canvas range
2441  *
2442  * @param p Parameter
2443  */
2444 void UpdateRange(dNdetaDrawer::RangeParam* p)
2445 {
2446   if (!p) {
2447     Warning("UpdateRange", "No parameters %p", p);
2448     return;
2449   }
2450   if (!p->fMasterAxis) {
2451     Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
2452     return;
2453   }
2454   Int_t    first = p->fMasterAxis->GetFirst();
2455   Int_t    last  = p->fMasterAxis->GetLast();
2456   Double_t x1    = p->fMasterAxis->GetBinCenter(first);
2457   Double_t x2    = p->fMasterAxis->GetBinCenter(last);
2458
2459   if (p->fSlave1Axis) {
2460     Int_t i1 = p->fSlave1Axis->FindBin(x1);
2461     Int_t i2 = p->fSlave1Axis->FindBin(x2);
2462     p->fSlave1Axis->SetRange(i1, i2);
2463     p->fSlave1Pad->Modified();
2464     p->fSlave1Pad->Update();
2465   }
2466   if (p->fSlave2Axis) {
2467     Int_t i1 = p->fSlave2Axis->FindBin(x1);
2468     Int_t i2 = p->fSlave2Axis->FindBin(x2);
2469     p->fSlave2Axis->SetRange(i1, i2);
2470     p->fSlave2Pad->Modified();
2471     p->fSlave2Pad->Update();
2472   }
2473   TCanvas*  c = gPad->GetCanvas();
2474   c->cd();
2475 }
2476
2477 //____________________________________________________________________
2478 /**
2479  * Called when user changes X range
2480  *
2481  * @param p Parameter
2482  */
2483 void RangeExec(dNdetaDrawer::RangeParam* p)
2484 {
2485   // Event types:
2486   //  51:   Mouse move
2487   //  53:
2488   //  1:    Button down
2489   //  21:   Mouse drag
2490   //  11:   Mouse release
2491   // dNdetaDrawer::RangeParam* p =
2492   //   reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
2493   Int_t event     = gPad->GetEvent();
2494   TObject *select = gPad->GetSelected();
2495   if (event == 53) {
2496     UpdateRange(p);
2497     return;
2498   }
2499   if (event != 11 || !select || select != p->fMasterAxis) return;
2500   UpdateRange(p);
2501 }
2502
2503 //=== Steering functions
2504 //==============================================
2505 /**
2506  * Display usage information
2507  *
2508  */
2509 void
2510 Usage()
2511 {
2512   printf("Usage: DrawdNdeta(FILE,TITLE,REBIN,OTHERS,FLAGS,"
2513          "SNN,SYS,TRIG,IPZMIN,IPZMAX,BASE,FMT)\n"
2514          "  const char* FILE   File name to open (\"forward_dndeta.root\")\n"
2515          "  const char* TITLE  Title to put on plot (\"\")\n"
2516          "  UShort_t    REBIN  Rebinning factor (1)\n"
2517          "  UShort_t    OTHERS Other data to draw - more below (0x7)\n"
2518          "  UShort_t    FLAGS  Visualisation flags - more below (0x7)\n"
2519          "  UShort_t    SYS    (optional) 1:pp, 2:PbPb, 3:pPb\n"
2520          "  UShort_t    SNN    (optional) sqrt(s_NN) in GeV\n"
2521          "  UShort_t    TRIG   (optional) 1: INEL, 2: INEL>0, 4: NSD, ...\n"
2522          "  Float_t     EFF    (optional) Trigger efficiency\n"
2523          "  Float_t     IPZMIN (optional) Least z coordinate of IP\n"
2524          "  Float_t     IPZMAX (optional) Largest z coordinate of IP\n"
2525          "  const char* BASE   (optional) base name of output files\n"
2526          "  UShort_t    FMT    (optional) Output formats\n"
2527          "\n"
2528          " OTHERS is a bit mask of\n"
2529          "  0x1     Show UA5 data (INEL,NSD, ppbar, 900GeV)\n"
2530          "  0x2     Show CMS data (NSD, pp)\n"
2531          "  0x4     Show published ALICE data (INEL,INEL>0,NSD, pp)\n"
2532          "  0x8     Show event genertor data\n"
2533          "\n"
2534          " FLAGS is a bit mask of\n"
2535          "  0x1     Show ratios of data to other data and possibly MC\n"
2536          "  0x2     Show left-right asymmetry\n"
2537          "  0x4     Show systematic error band\n"
2538          "  0x8     Show individual ring results (INEL only)\n"
2539          "  0x10    Cut edges when rebinning\n"
2540          "  0x20    Remove FMDxO points\n"
2541          "  0x40    Do not make our own canvas\n"
2542          "  0x80    Force use of MB\n"
2543          "  0x100   Mirror data\n"
2544          "  0x200   Apply `final MC' correction\n"
2545          "  0x400   Apply `Emperical' correction\n"
2546          "  0x800   Export results to script\n"
2547          "  0x1000  Add code to do combined zooms on eta axis\n"
2548          "  0x2000  Assume old-style input\n"
2549          "\n"
2550          " FMD is a bit mask of\n"
2551          "  0x1     Make PNG output\n"
2552          "  0x2     Make PDF output\n"
2553          "  0x4     Make ROOT file output\n"
2554          "  0x8     Make ROOT script output\n"
2555          "\n"
2556          "0x200 requires the file forward_dndetamc.root\n"
2557          "0x400 requires the file EmpiricalCorrection.root\n"
2558          "To specify that you want ratios, force MB, apply empirical "
2559          "correction, and export to script, set flags to\n\n"
2560          "  0x1|0x80|0x400|0x800=0xC81\n\n"
2561          );
2562 }
2563
2564 //____________________________________________________________________
2565 /**
2566  * Draw @f$ dN/d\eta@f$
2567  *
2568  * @param filename  File name
2569  * @param title     Title
2570  * @param rebin     Rebinning factor
2571  * @param others    What other data to show
2572  * @param flags     Flags
2573  * @param sNN       (optional) Collision energy [GeV]
2574  * @param sys       (optional) Collision system (1: pp, 2: PbPb)
2575  * @param trg       (optional) Trigger (1: INEL, 2: INEL>0, 4: NSD)
2576  * @param vzMin     Least @f$ v_z@f$
2577  * @param vzMax     Largest @f$ v_z@f$
2578  *
2579  * @ingroup pwglf_forward_dndeta
2580  */
2581 void
2582 DrawdNdeta(const char* filename="forward_dndeta.root",
2583            const char* title="",
2584            UShort_t    rebin=5,
2585            UShort_t    others=0x7,
2586            UShort_t    flags=0x187,
2587            UShort_t    sNN=0,
2588            UShort_t    sys=0,
2589            UShort_t    trg=0,
2590            Float_t     eff=0,
2591            Float_t     vzMin=999,
2592            Float_t     vzMax=-999,
2593            const char* base="",
2594            UShort_t    outflg=0xf)
2595 {
2596   TString fname(filename);
2597   fname.ToLower();
2598   if (fname.CompareTo("help") == 0 ||
2599       fname.CompareTo("--help") == 0) {
2600     Usage();
2601     return;
2602   }
2603   dNdetaDrawer* pd = new dNdetaDrawer;
2604   dNdetaDrawer& d = *pd;
2605   d.SetRebin(rebin);
2606   d.SetTitle(title);
2607   d.SetShowOthers(others);
2608   d.SetShowRatios(flags & 0x1);
2609   d.SetShowLeftRight(flags & 0x2);
2610   d.SetForwardSysError(flags & 0x4 ? 0.076 : 0);
2611   d.SetShowRings(flags & 0x8);
2612   d.SetCutEdges(flags & 0x10);
2613   d.fRemoveOuters = (flags & 0x20);
2614   d.SetExport(flags & 0x40);
2615   d.SetForceMB(flags & 0x80);
2616   d.SetMirror(flags & 0x100);
2617   d.SetFinalMC(flags & 0x200 ? "forward_dndetamc.root" : "");
2618   d.SetEmpirical(flags & 0x400 ? "EmpiricalCorrection.root" : "");
2619   d.SetExport(flags & 0x800);
2620   d.SetAddExec(flags & 0x1000);
2621   d.SetOld(flags & 0x2000);
2622   d.SetBase(base);
2623   d.SetSaveAsPNG(outflg & 0x1);
2624   d.SetSaveAsPDF(outflg & 0x2);
2625   d.SetSaveAsROOT(outflg & 0x4);
2626   d.SetSaveAsScript(outflg & 0x8);
2627   // d.fClusterScale = "1.06 -0.003*x +0.0119*x*x";
2628   // Do the below if your input data does not contain these settings
2629   if (sNN > 0) d.SetSNN(sNN);     // Collision energy per nucleon pair (GeV)
2630   if (sys > 0) d.SetSys(sys);     // Collision system (1:pp, 2:PbPB)
2631   if (trg > 0) d.SetTrigger(trg); // Collision trigger (1:INEL, 2:INEL>0, 4:NSD)
2632   if (eff > 0) d.SetTriggerEfficiency(eff); // Trigger efficiency
2633   if (vzMin < 999 && vzMax > -999)
2634     d.SetVertexRange(vzMin,vzMax); // Collision vertex range (cm)
2635   d.Run(filename);
2636 }
2637 //____________________________________________________________________
2638 //
2639 // EOF
2640 //
2641