PWGLF/FORWARD/analysis2/DrawdNdeta.C

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