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