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