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