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