PWGLF/FORWARD/analysis2/DrawdNdeta.C

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