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