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