Fix some documentation issues

[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / UnfoldMultDists.C

/**
 * @file   UnfoldMultDists.C
 * @author Christian Holm Christensen <cholm@nbi.dk>
 * @date   Tue Nov 12 09:25:52 2013
 * 
 * @brief  A class to do unfolding 
 * 
 * 
 * @ingroup pwglf_forward_multdist
 */
#include <TFile.h>
#include <TList.h>
#include <TH1.h>
#include <TH2.h>
#include <THStack.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TClass.h>
#include <TRegexp.h>
#include <TMath.h>
#include <TParameter.h>
#include <TMultiGraph.h>
#include <TGraphAsymmErrors.h>
#include "RooUnfold.h"
#include "RooUnfoldResponse.h"
#include <fstream>

/**
 * Class to do unfolding of raw histograms produced by AliForwardMultDists 
 * 
 * @ingroup pwglf_forward_multdist
 */
struct Unfolder
{
  /**
   * Colours used 
   * 
   */
  enum {
    kColorMeasured = kOrange-2, 
    kColorTruth    = kBlue-3, 
    kColorAccepted = kMagenta-3,
    kColorTrgVtx   = kBlack, 
    kColorUnfolded = kOrange+2,
    kColorCorrected= kRed+2, 
    kColorError    = kBlue-10,
    kColorALICE    = kPink+1, 
    kColorCMS      = kGreen+2
  };

  /** 
   * Constructor 
   */
  Unfolder() {}
  /** 
   * Get a top collection from a file
   * 
   * @param fileName Name of file 
   * @param results  Wheter it's the results collection or not 
   * 
   * @return Collection or null
   */
  static TCollection* GetTop(const TString& fileName, Bool_t results=false)
  {
    TFile* file = TFile::Open(fileName, "READ");
    if (!file) { 
      Error("GetTop", "Failed to open %s", fileName.Data());
      return 0;
    }
    TCollection* ret = 0;
    TString cName(Form("ForwardMult%s", results ? "Results" : "Sums"));
    file->GetObject(cName, ret);
    if (!ret) 
      Error("GetTop", "Failed to get collection %s from %s", 
	    cName.Data(), fileName.Data());
    file->Close();
    return ret;
  }
  /** 
   * Get an object from a collection 
   * 
   * @param c    Collection
   * @param name Name of object
   * @param cl   Possible class to check against
   * @param verbose  Be verbose
   * 
   * @return Pointer to object or null
   */
  static TObject* GetObject(TCollection* c, const TString& name, 
			    TClass* cl, Bool_t verbose=true)
  {
    TObject* o = c->FindObject(name);
    if (!o) { 
      if (verbose)
	Warning("GetObject", "%s not found in %s", name.Data(), c->GetName());
      return 0;
    }
    if (cl && !o->IsA()->InheritsFrom(cl)) {
      if (verbose) 
	Warning("GetCollection", "%s is not a %s but a %s", 
		name.Data(), cl->GetName(), o->ClassName());
      return 0;
    }
    return o;
  }
  /** 
   * Get a collection 
   * 
   * @param c        Parent collection
   * @param name     Name of object to findf
   * @param verbose  Be verbose
   * 
   * @return 
   */
  static TCollection* GetCollection(TCollection*   c, 
				    const TString& name, 
				    Bool_t         verbose=-true)
  {
    return static_cast<TCollection*>(GetObject(c, name, 
					       TCollection::Class(),
					       verbose));
  }
  /** 
   * Get a 1D histogram from a collection
   * 
   * @param c    Collection
   * @param name Nanme of histogram
   * @param verbose  Be verbose
   * 
   * @return Pointer to object or null
   */
  static TH1* GetH1(TCollection* c, const TString& name, Bool_t verbose=true) 
  {
    return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verbose));
  }
  /** 
   * Get a 2D histogram from a collection
   * 
   * @param c    Collection
   * @param name Nanme of histogram
   * @param verbose  Be verbose
   * 
   * @return Pointer to object or null
   */
  static TH2* GetH2(TCollection* c, const TString& name, Bool_t verbose=true) 
  {
    return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verbose));
  }
  /** 
   * Get an unsigned short parameter from the collection 
   * 
   * @param c    Collection
   * @param name Parameter name 
   * @param v    Value
   * 
   * @return Value 
   */
  static void GetParameter(TCollection* c, const TString& name, UShort_t& v)
  {
    TObject* o = GetObject(c, name, TParameter<int>::Class(), true);
    v = (!o ? 0 : o->GetUniqueID());
  }
  /** 
   * Get an unsigned short parameter from the collection 
   * 
   * @param c    Collection
   * @param name Parameter name 
   * @param v    Value
   * 
   * @return Value 
   */
  static void GetParameter(TCollection* c, const TString& name, ULong_t& v)
  {
    TObject* o = GetObject(c, name, TParameter<long>::Class(), true);
    v = (!o ? 0 : o->GetUniqueID());
  }
  /** 
   * Get an unsigned short parameter from the collection 
   * 
   * @param c    Collection
   * @param name Parameter name 
   * @param v    Value
   * 
   * @return Value 
   */
  static void GetParameter(TCollection* c, const TString& name, Double_t& v)
  {
    TObject* o = GetObject(c, name, TParameter<double>::Class(), true);
    v = (!o ? 0 : static_cast<TParameter<double>*>(o)->GetVal());
  }
  /** 
   * Get the method identifier 
   * 
   * @param method Method 
   * 
   * @return Method identifier 
   */    
  static UInt_t MethodId(TString& method) 
  {
    struct Method { 
      UInt_t  id;
      TString name;
    };
    const Method methods[] = { {RooUnfold::kNone,    "None"},
			       {RooUnfold::kBayes,   "Bayes"},
			       {RooUnfold::kSVD,     "SVD"},
			       {RooUnfold::kBinByBin,"BinByBin"},
			       {RooUnfold::kTUnfold, "TUnfold"},
			       {RooUnfold::kInvert,  "Invert"},
			       {RooUnfold::kDagostini,"Dagostini"}, 
			       {0xDeadBeef,           "unknown"} };
    const Method* pMethod = methods;
    while (pMethod->id != 0xDeadBeef) {
      if (method.BeginsWith(pMethod->name, TString::kIgnoreCase)) {
	method = pMethod->name;
	break;
      }
      pMethod++;
    }
    if (pMethod->id == 0xDeadBeef) 
      Error("MethodId", "Unknown unfolding method: %s", method.Data());

    return pMethod->id;
  }
  
  /** 
   * Run the unfolding and correction task. 
   *
   * The @a measuredFile is assumed to have the structure 
   *
   * @verbatim
   *     /- ForwardMultSums         (TCollection)
   *          |- [type]             (TCollection)
   *          |    |- [bin]         (TCollection) 
   *          |    |    `- rawDist  (TH1)
   *          |    |- [bin]
   *          |    ...
   *          |- [type]
   *          ...
   * @endverbatim 
   * 
   * and @a corrFile is assumed to have the structure 
   *
   * @verbatim
   *     /- ForwardMultResults            (TCollection)
   *          |- [type]                   (TCollection)
   *          |    |- [bin]               (TCollection) 
   *          |    |    |- truth          (TH1)
   *          |    |    |- truthAccepted  (TH1)
   *          |    |    |- triggerVertex  (TH1)
   *          |    |    `- response       (TH2)
   *          |    |- [bin]
   *          |    ...
   *          |- [type]
   *          ...
   * @endverbatim 
   *
   * where @c [type] is one of <i>symmetric</i>, <i>positive</i>,
   * <i>negative</i>, or <i>other</i>, and [bin] is the @f$ \eta@f$
   * bin named like
   *
   * @verbatim
   *   [bin]          := [eta_spec] _ [eta_spec]
   *   [eta_spec]     := [sign_char] [integer_part] d [decimal_part]
   *   [sign_part]    := p          positive eta 
   *                  |  m          negative eta 
   *   [integer_part] := [number]
   *   [decimal_part] := [number] [number]
   *   [number]       := 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 
   * @endverbatim
   *
   * That is, the bin @f$ -3\le\eta\ge3@f$ is labeled
   * <b>m3d00_p3d00</b>, @f$ 0\le\eta\ge2.5@f$ is <b>p0d00_p2d50</b> 
   *
   * @a measuredFile and @a corrFile can point to the same file.  If
   * @a corrFile is not specified, it is assumed that @a measuredFile
   * has the expected @a corrFile @e in @e addition to the
   * expected content of that file.
   * 
   * @param measuredFile Name of file containing measured data
   * @param corrFile     Name of file containing correction data
   * @param method       Unfolding method to use 
   * @param regParam     Regularization parameter 
   */  
  void Run(const TString& measuredFile, const TString& corrFile,
	   const TString& method="Bayes", Double_t regParam=4)
  {
    // Get the input collections 
    if (measuredFile.IsNull()) {
      Error("Run", "No measurements given");
      return;
    }
    TCollection* mTop = GetTop(measuredFile, false);
    TCollection* cTop = GetTop(corrFile.IsNull() ? measuredFile : corrFile, 
			       true);
    if (!mTop || !cTop) return;

    // Get some info from the input collection 
    UShort_t sys;
    UShort_t sNN; 
    ULong_t  trig; 
    Double_t minZ; 
    Double_t maxZ;
    GetParameter(mTop, "sys",     sys);	  
    GetParameter(mTop, "sNN",     sNN);	  
    GetParameter(mTop, "trigger", trig);
    GetParameter(mTop, "minIpZ",  minZ); 
    GetParameter(mTop, "maxIpZ",  maxZ); 
    if (sys == 0 || sNN == 0) 
      Warning("Run", "System (%d) and/or collision energy (%d) unknown", 
	      sys, sNN);
    
    // Open the output file 
    TFile* out = TFile::Open("forward_unfolded.root", "RECREATE");
    if (!out) { 
      Error("Run", "Failed to open output file");
      return;
    }

    // Decode method option and store in file 
    TString meth(method);
    UInt_t  mId = MethodId(meth);
    if (mId == 0xDeadBeef) return;

    // Store information 
    SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
		    corrFile.IsNull());

    // Load other data 
    TString savPath(gROOT->GetMacroPath());
    gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2/scripts",
                             gROOT->GetMacroPath()));
    // Always recompile 
    if (!gROOT->GetClass("OtherPNch"))
      gROOT->LoadMacro("OtherPNchData.C++");
    gROOT->SetMacroPath(savPath);

    // Loop over the input 
    const char*  inputs[] = { "symmetric", "positive", "negative", 0 };
    const char** pinput   = inputs;
    while (*pinput) { 
      TCollection* mInput = GetCollection(mTop, *pinput, false);
      TCollection* cInput = GetCollection(cTop, *pinput, false);
      if (mInput && cInput)
	ProcessType(mInput, cInput, mId, regParam, out,
		    sys, sNN);
      pinput++;
    }      

    out->Write();
    // out->ls();
    out->Close();

    SaveSummarize();
  }
  /** 
   * Append an & to a string and the next term.
   * 
   * @param trg  Output string
   * @param what Term
   */
  static void AppendAnd(TString& trg, const TString& what)
  {
    if (!trg.IsNull()) trg.Append(" & ");
    trg.Append(what);
  }
  /** 
   * Store some information on the output
   * 
   * @param dir      Where to store
   * @param method   Method used
   * @param mId      Method identifier 
   * @param regParam Regularization parameter 
   * @param sys      Collision system
   * @param sNN      Center of mass energy 
   * @param trigger  Trigger mask 
   * @param minIpZ   Least z coordinate of interaction point
   * @param maxIpZ   Largest z coordinate of interaction point
   * @param self     Self-consistency check
   */
  void SaveInformation(TDirectory* dir, 
		       const TString& method,
		       Int_t          mId,
		       Double_t       regParam, 
		       UShort_t       sys, 
		       UShort_t       sNN,
		       UInt_t         trigger,
		       Double_t       minIpZ, 
		       Double_t       maxIpZ, 
		       Bool_t         self) const
  {
    dir->cd();

    TParameter<bool>* pM = new TParameter<bool>("self", self);
    pM->SetBit(BIT(19));
    pM->Write();

    TNamed* outMeth = new TNamed("method", method.Data());
    outMeth->SetUniqueID(mId);
    outMeth->Write();

    TParameter<double>* pR = new TParameter<double>("regParam", regParam);
    pR->SetBit(BIT(19));
    pR->Write();
    
    TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
    TNamed* pS = new TNamed("sys", tS.Data()); pS->SetUniqueID(sys);
    pS->Write();
    
    TString tE;
    if      (sNN <  1000) tE = Form("%dGeV", sNN);
    else if (sNN <  3000) tE = Form("%4.2fTeV", float(sNN)/1000);
    else                  tE = Form("%dTeV", sNN/1000);
    TNamed* pE = new TNamed("sNN", tE.Data()); pE->SetUniqueID(sNN);
    pE->Write();
    
    TString tT;
    /** 
     * Bits of the trigger pattern
     */
    enum { 
      /** In-elastic collision */
      kInel        = 0x0001, 
      /** In-elastic collision with at least one SPD tracklet */
      kInelGt0     = 0x0002, 
      /** Non-single diffractive collision */
      kNSD         = 0x0004, 
      /** Empty bunch crossing */
      kEmpty       = 0x0008, 
      /** A-side trigger */
      kA           = 0x0010, 
      /** B(arrel) trigger */
      kB           = 0x0020, 
      /** C-side trigger */
      kC           = 0x0080,  
      /** Empty trigger */
      kE           = 0x0100,
      /** pileup from SPD */
      kPileUp      = 0x0200,    
      /** true NSD from MC */
      kMCNSD       = 0x0400,    
      /** Offline MB triggered */
      kOffline     = 0x0800,
      /** At least one SPD cluster */ 
      kNClusterGt0 = 0x1000,
      /** V0-AND trigger */
      kV0AND       = 0x2000, 
      /** Satellite event */
      kSatellite   = 0x4000
    };
    if ((trigger & kInel)        != 0x0) AppendAnd(tT, "INEL");
    if ((trigger & kInelGt0)     != 0x0) AppendAnd(tT, "INEL>0");
    if ((trigger & kNSD)         != 0x0) AppendAnd(tT, "NSD");
    if ((trigger & kV0AND)       != 0x0) AppendAnd(tT, "V0AND");
    if ((trigger & kA)           != 0x0) AppendAnd(tT, "A");
    if ((trigger & kB)           != 0x0) AppendAnd(tT, "B");
    if ((trigger & kC)           != 0x0) AppendAnd(tT, "C");
    if ((trigger & kE)           != 0x0) AppendAnd(tT, "E");
    if ((trigger & kMCNSD)       != 0x0) AppendAnd(tT, "MCNSD");
    if ((trigger & kNClusterGt0) != 0x0) AppendAnd(tT, "NCluster>0");
    if ((trigger & kSatellite)   != 0x0) AppendAnd(tT, "Satellite");
    TNamed* pT = new TNamed("trigger", tT.Data()); pT->SetUniqueID(trigger);
    pT->Write();
    
    TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
    pY->SetBit(BIT(19));
    pY->Write();

    TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
    pZ->SetBit(BIT(19));
    pZ->Write();
  }
  /** 
   * Save a script to do a summary of this step 
   * 
   */		       
  void SaveSummarize()
  {
    std::ofstream f("SummarizeUnfold.C");
    f << "void SummarizeUnfold(const char* file=\"forward_unfolded.root\")\n"
      << "{\n"
      << "  const char* fwd=\"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n"
      << "  gROOT->LoadMacro(Form(\"%s/DrawUnfoldedSummary.C\",fwd));\n"
      << "  DrawUnfoldedSummary(file);\n"
      << "}\n"
      << "// EOF" << std::endl;
    f.close();
  }
  /** 
   * Process a single type - i.e., one of <i>symmetric</i>,
   * <i>positive</i>, <i>negative</i>, or <i>other</i> - by looping
   * over all contained objects and call ProcessBin for each found
   * bin.
   * 
   * @param measured     Input collection of measured data
   * @param corrections  Input collection of correction data
   * @param method       Unfolding method to use 
   * @param regParam     Regularisation parameter
   * @param out          Output directory. 
   * @param sys          Collision system
   * @param sNN          Collision energy 
   */
  void ProcessType(TCollection* measured, 
		   TCollection* corrections,
		   UInt_t       method,
		   Double_t     regParam,
		   TDirectory*  out,
		   UShort_t     sys, 
		   UShort_t     sNN)
  {
    Printf("  Processing %s ...", measured->GetName());
    TDirectory* dir = out->mkdir(measured->GetName());
    
    // Make some summary stacks 
    THStack*  allMeasured  = new THStack("measured",      
					 "Measured P(#it{N}_{ch})");
    THStack*  allTruth     = new THStack("truth",        
					 "MC 'truth' P(#it{N}_{ch})");
    THStack*  allTruthA    = new THStack("truthAccepted",
					 "MC 'truth' accepted P(#it{N}_{ch})");
    THStack*  allUnfolded  = new THStack("unfolded",
					 "Unfolded P(#it{N}_{ch})");
    THStack*  allCorrected = new THStack("corrected",
					 "Corrected P(#it{N}_{ch})");
    THStack*  allRatio     = (sys != 1 ? 0 : 
			      new THStack("ratios", "Ratios to other"));
    TMultiGraph* allALICE  = (sys != 1 ? 0 : 
			      new TMultiGraph("alice", "ALICE Published"));
    TMultiGraph* allCMS    = (sys != 1 ? 0 : 
			      new TMultiGraph("cms", "CMS Published"));

    // Loop over the list of objects. 
    static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
    TIter          next(measured);
    TObject*       o = 0;
    Int_t          i = 0;
    Double_t       r = regParam;
    while ((o = next())) {
      // Go back to where we where 
      dir->cd();
      
      // if not a collection, don't bother 
      if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
    
      // If it doesn't match our regular expression, don't bother 
      TString n(o->GetName());
      if (n.Index(regex) == kNPOS) { 
	// Warning("ScanType", "%s in %s doesn't match eta range regexp", 
	//         n.Data(), real->GetName());
	continue;
      }
      TCollection* mBin = static_cast<TCollection*>(o);
      TCollection* cBin = GetCollection(corrections, n.Data());
      if (!cBin) continue;

      THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
      if (!binS) continue;

      TH1* result = 0;
      Bin2Stack(binS, i, allMeasured, allTruth, allTruthA, 
		allUnfolded, allCorrected, result);

      TGraph* alice = 0;
      TGraph* cms   = 0;
      Other2Stack(o->GetName(), i, sNN, allALICE, allCMS, alice, cms);
      Ratio2Stack(i, result, alice, cms, allRatio);
      i++;
    }
    dir->Add(allMeasured);
    dir->Add(allTruth);
    dir->Add(allTruthA);
    dir->Add(allUnfolded);
    dir->Add(allCorrected);
    if (allALICE && allALICE->GetListOfGraphs()) {
      if (allALICE->GetListOfGraphs()->GetEntries() > 0)
	dir->Add(allALICE);
      else 
	delete allALICE;
    }
    if (allCMS && allCMS->GetListOfGraphs()) {
      if (allCMS->GetListOfGraphs()->GetEntries() > 0) 
	dir->Add(allCMS);
      else 
	delete allCMS;
    }
    if (allRatio && allRatio->GetHists()) { 
      if (allRatio->GetHists()->GetEntries() > 0) 
	dir->Add(allRatio);
      else 
	delete allRatio;
    }
  }
  /** 
   * Process a single eta bin 
   * 
   * @param measured     Input collection of measured data
   * @param corrections  Input collection of correction data
   * @param method       Unfolding method to use 
   * @param regParam     Regularisation parameter
   * @param out          Output directory. 
   *
   * @return Stack of histograms or null 
   */
  THStack* ProcessBin(TCollection* measured, 
		      TCollection* corrections, 
		      UInt_t       method,
		      Double_t     regParam, 
		      TDirectory*  out)
  {
    Printf("   Processing %s ...", measured->GetName());
    // Try to get the data 
    TH1* inRaw    = GetH1(measured, "rawDist");
    TH1* inTruth  = GetH1(corrections, "truth");
    TH1* inTruthA = GetH1(corrections, "truthAccepted");
    TH1* inTrgVtx = GetH1(corrections, "triggerVertex");
    TH2* inResp   = GetH2(corrections, "response");
    if (!inRaw || !inTruth || !inTruthA || !inTrgVtx || !inResp) 
      return 0;
    
    // Make output directory
    TDirectory* dir = out->mkdir(measured->GetName());
    dir->cd();

    // Copy the input to the output 
    TH1* outRaw    = static_cast<TH1*>(inRaw    ->Clone("measured"));
    TH1* outTruth  = static_cast<TH1*>(inTruth  ->Clone("truth"));
    TH1* outTruthA = static_cast<TH1*>(inTruthA ->Clone("truthAccepted"));
    TH1* outTrgVtx = static_cast<TH1*>(inTrgVtx ->Clone("triggerVertex"));
    TH2* outResp   = static_cast<TH2*>(inResp   ->Clone("response"));

    // Make our response matrix 
    RooUnfoldResponse matrix(0, 0, inResp);
    
    // Store regularization parameter 
    Double_t             r        = regParam;
    RooUnfold::Algorithm algo     = (RooUnfold::Algorithm)method;
    RooUnfold*           unfolder = RooUnfold::New(algo, &matrix, inRaw, r);
    unfolder->SetVerbose(1);

    // Do the unfolding and get the result
    TH1* res = unfolder->Hreco();
    res->SetDirectory(0);

    // Make a copy to store on the output 
    TH1* outUnfold = static_cast<TH1*>(res->Clone("unfolded"));
    TString tit(outUnfold->GetTitle());
    tit.ReplaceAll("Unfold Reponse matrix", "Unfolded P(#it{N}_{ch})");
    outUnfold->SetTitle(tit);

    // Clone the unfolded results and divide by the trigger/vertex
    // bias correction
    TH1* outCorr   = static_cast<TH1*>(outUnfold->Clone("corrected"));
    outCorr->Divide(inTrgVtx);
    tit.ReplaceAll("Unfolded", "Corrected");
    outCorr->SetTitle(tit);

    // Now normalize the output to integral=1 
    TH1*  hists[] = { outRaw, outUnfold, outCorr, 0 };
    TH1** phist   = hists;
    while (*phist) { 
      TH1* h = *phist;
      if (h) { 
	Double_t intg = h->Integral(1, h->GetXaxis()->GetXmax());
	h->Scale(1. / intg, "width");
      }
      phist++;
    }
    
    // And make ratios
    TH1* ratioTrue = static_cast<TH1*>(outCorr->Clone("ratioCorrTruth"));
    tit = ratioTrue->GetTitle();
    tit.ReplaceAll("Corrected", "Corrected/MC 'truth'");
    ratioTrue->SetTitle(tit);
    ratioTrue->Divide(outTruth);
    ratioTrue->SetYTitle("P_{corrected}(#it{N}_{ch})/P_{truth}(#it{N}_{ch})");

    TH1* ratioAcc  = static_cast<TH1*>(outUnfold->Clone("ratioUnfAcc"));
    tit = ratioAcc->GetTitle();
    tit.ReplaceAll("Unfolded", "Unfolded/MC selected");
    ratioAcc->SetTitle(tit);
    ratioAcc->Divide(outTruthA);
    ratioAcc->SetYTitle("P_{unfolded}(#it{N}_{ch})/P_{MC}(#it{N}_{ch})");
    

    // Make a stack 
    tit = measured->GetName();
    tit.ReplaceAll("m", "-");
    tit.ReplaceAll("p", "+");
    tit.ReplaceAll("d", ".");
    tit.ReplaceAll("_", "<#it{#eta}<");
    THStack* stack = new THStack("all", tit);
    stack->Add(outTruth,  "E2");
    stack->Add(outTruthA, "E2");
    stack->Add(outRaw,    "E1");
    stack->Add(outUnfold, "E1");
    stack->Add(outCorr,   "E1");
    dir->Add(stack);

    // Rest of the function is devoted to making the output look nice 
    outRaw   ->SetDirectory(dir); 
    outTruth ->SetDirectory(dir);  
    outTruthA->SetDirectory(dir);  
    outTrgVtx->SetDirectory(dir);  
    outResp  ->SetDirectory(dir);  
    outUnfold->SetDirectory(dir);   
    outCorr  ->SetDirectory(dir); 

    outRaw   ->SetMarkerStyle(20);  // Measured is closed
    outTruth ->SetMarkerStyle(24);  // MC is open
    outTruthA->SetMarkerStyle(24);  // MC is open
    outTrgVtx->SetMarkerStyle(20);  // Derived is closed
    outUnfold->SetMarkerStyle(20);  // Derived is closed   
    outCorr  ->SetMarkerStyle(20);  // Derived is closed 

    outRaw   ->SetMarkerSize(0.9); 
    outTruth ->SetMarkerSize(1.6);  
    outTruthA->SetMarkerSize(1.4);  
    outTrgVtx->SetMarkerSize(1.0);  
    outUnfold->SetMarkerSize(0.9);   
    outCorr  ->SetMarkerSize(1.0);
 
    outRaw   ->SetMarkerColor(kColorMeasured); 
    outTruth ->SetMarkerColor(kColorTruth);  
    outTruthA->SetMarkerColor(kColorAccepted);  
    outTrgVtx->SetMarkerColor(kColorTrgVtx);  
    outUnfold->SetMarkerColor(kColorUnfolded);   
    outCorr  ->SetMarkerColor(kColorCorrected); 

    outRaw   ->SetFillColor(kColorError);     
    outTruth ->SetFillColor(kColorError);  
    outTruthA->SetFillColor(kColorError);  
    outTrgVtx->SetFillColor(kColorError);  
    outUnfold->SetFillColor(kColorError);   
    outCorr  ->SetFillColor(kColorError); 

    outRaw   ->SetFillStyle(0); 
    outTruth ->SetFillStyle(1001);  
    outTruthA->SetFillStyle(1001);  
    outTrgVtx->SetFillStyle(0);  
    outUnfold->SetFillStyle(0);   
    outCorr  ->SetFillStyle(0);

    outRaw   ->SetLineColor(kBlack); 
    outTruth ->SetLineColor(kBlack);  
    outTruthA->SetLineColor(kBlack);  
    outTrgVtx->SetLineColor(kBlack);  
    outUnfold->SetLineColor(kBlack);   
    outCorr  ->SetLineColor(kBlack); 

    // Legend 
    TLegend* l = StackLegend(stack);
    l->AddEntry(outRaw,     "Raw",                 "lp");
    l->AddEntry(outTruth,   "MC 'truth'",          "fp");
    l->AddEntry(outTruthA,  "MC 'truth' accepted", "fp");
    l->AddEntry(outUnfold,  "Unfolded",            "lp");
    l->AddEntry(outCorr,    "Corrected",           "lp");

    return stack;
  }
  static void BinAttributes(Int_t i,
			    Int_t&    open, 
			    Int_t&    closed, 
			    Float_t&  size,
			    Double_t& factor) 
  {
    // --- Setup for markers -----------------------------------------
    const Int_t   nMarkers = 7;
    const Int_t   aClosed[] = { 20,  21,  22,  23,  29,  33,  34  };
    const Int_t   aOpen[]   = { 24,  25,  26,  32,  30,  27,  28  };
    const Float_t aSize[]   = { 1.1, 1.0, 1.2, 1.2, 1.2, 1.2, 1.0 };
    Int_t         j         = i % nMarkers;

    open   = aOpen[j];
    closed = aClosed[j];
    size   = aSize[j];
    factor = TMath::Power(10, i);
  }
  /** 
   * Add the bin histograms to our summary stacks 
   * 
   * @param bin       Bin stack
   * @param i         Current off-set in the stacks 
   * @param measured  All measured @f$ P(N_{ch})@f$ 
   * @param truth     All MC truth @f$ P(N_{ch})@f$ 
   * @param accepted  All MC accepted @f$ P(N_{ch})@f$ 
   * @param unfolded  All unfolded @f$ P(N_{ch})@f$ 
   * @param corrected All corrected @f$ P(N_{ch})@f$ 
   * @param result    The result in this bin
   */
  void Bin2Stack(const THStack* bin, Int_t i, 
		 THStack* measured, 
		 THStack* truth, 
		 THStack* accepted, 
		 THStack* unfolded,
		 THStack* corrected,
		 TH1*&    result)
  {
    Int_t open, closed;
    Double_t factor; 
    Float_t  size;
    BinAttributes(i, open, closed, size, factor);

    TIter next(bin->GetHists());
    TH1*  h = 0;
    while ((h = static_cast<TH1*>(next()))) {
      THStack* tmp = 0;
      Int_t    col = h->GetMarkerColor();
      Int_t    sty = 0;
      switch (col) { 
      case kColorMeasured:  tmp = measured;   sty = closed;  break;
      case kColorTruth:     tmp = truth;      sty = open;    break;
      case kColorAccepted:  tmp = accepted;   sty = open;    break;
      case kColorUnfolded:  tmp = unfolded;   sty = closed;  break;
      case kColorCorrected: tmp = corrected;  sty = closed;  break;
      default: continue; 
      }
      // Now clone, and add to the appropriate stack 
      TH1* cln = static_cast<TH1*>(h->Clone(h->GetName()));
      cln->SetDirectory(0);
      cln->SetMarkerStyle(sty);
      cln->SetMarkerSize(size);
      cln->Scale(factor); // Scale by 10^i
      if (col == kColorCorrected) result = cln;

      // Make sure we do not get the old legend 
      TObject* tst = cln->FindObject("legend");
      if (tst) cln->GetListOfFunctions()->Remove(tst);

      tmp->Add(cln, next.GetOption());
    }
    
    // Add entries to our stacks 
    TString   txt      = bin->GetTitle();
    if      (i == 0) txt.Append(" (#times1)");
    else if (i == 1) txt.Append(" (#times10)");
    else             txt.Append(Form(" (#times10^{%d})", i));
    THStack*  stacks[] = { measured, truth, accepted, unfolded, corrected, 0 };
    THStack** pstack   = stacks;
    while (*pstack) { 
      TLegend* leg = StackLegend(*pstack);
      pstack++;
      if (!leg) continue;
      
      TObject* dummy = 0;
      TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
      e->SetMarkerStyle(closed);
      e->SetMarkerSize(1.2*size);
      e->SetMarkerColor(kBlack);
      e->SetFillColor(0);
      e->SetFillStyle(0);
      e->SetLineColor(kBlack);
    }
  }
  /** 
   * Add distributions from other experiments to stacks 
   * 
   * @param name     Name of current bin 
   * @param i        Index of current bin
   * @param sNN      Center of mass energy 
   * @param allALICE Stack of ALICE data 
   * @param allCMS   Stack of CMS data 
   * @param alice    Possible ALICE result on return
   * @param cms      Possible CMS result on return
   */
  void Other2Stack(const TString& name, Int_t i,
		   UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
		   TGraph*& alice, TGraph*& cms) 
  {
    if (!allALICE && !allCMS) return;

    TString tmp(name);
    tmp.ReplaceAll("p", "+");
    tmp.ReplaceAll("m", "-");
    tmp.ReplaceAll("_", " ");
    tmp.ReplaceAll("d", ".");
    TObjArray* tokens = tmp.Tokenize(" ");
    if (!tokens || tokens->GetEntriesFast() < 2) { 
      Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
      if (tokens) tokens->Delete();
      return;
    }
    Double_t eta1 = static_cast<TObjString*>(tokens->At(0))->String().Atof();
    Double_t eta2 = static_cast<TObjString*>(tokens->At(1))->String().Atof();
    tokens->Delete();
    
    if (TMath::Abs(eta2-eta1) > 1e3) 
      // Not symmetric bin 
      return;
    Double_t aEta = TMath::Abs(eta1);

    Int_t open, closed;
    Double_t factor; 
    Float_t  size;
    BinAttributes(i, open, closed, size, factor);

    if (allALICE) {
      TGraphAsymmErrors* g = GetOther(1, aEta, sNN, factor);
      if (g) {
	g->SetMarkerStyle(closed);
	g->SetMarkerColor(kColorALICE);
	g->SetMarkerSize(size);
	allALICE->Add(g, "p same");
	alice = g;
      }
    }
    if (allCMS) {
      TGraphAsymmErrors* g = GetOther(0, aEta, sNN, factor);
      if (g) {
	g->SetMarkerStyle(closed);
	g->SetMarkerColor(kColorCMS);
	g->SetMarkerSize(size);
	allCMS->Add(g, "p same");
	cms = g;
      }
    }
  }
  /** 
   * Create ratios to other data 
   * 
   * @param ib      Bin number  
   * @param res     Result
   * @param alice   ALICE result if any
   * @param cms     CMS result if any
   * @param all     Stack to add ratio to 
   */
  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
  {
    if (!all || !res || !(alice || cms)) return;

    Int_t        off  = 5*ib;
    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
    TGraph**     pg   = gs;
    while (*pg) { 
      TGraph*     g = *pg;
      const char* n = (g == alice ? "ALICE" : "CMS");

      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
      TString tit(r->GetTitle());
      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
      r->SetTitle(tit);
      r->SetMarkerColor(g->GetMarkerColor());
      r->SetLineColor(g->GetLineColor());

      TObject* tst = r->FindObject("legend");
      if (tst) r->GetListOfFunctions()->Remove(tst);

      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
	Double_t c = r->GetBinContent(i);
	Double_t e = r->GetBinError(i);
	Double_t o = g->Eval(r->GetBinCenter(i));
	if (o < 1e-12) { 
	  r->SetBinContent(i, 0);
	  r->SetBinError(i, 0);
	  continue;
	}
	r->SetBinContent(i, (c - o) / o + off);
	r->SetBinError(i, e / o);
      }
      all->Add(r);
      pg++;
    }
    TLegend* leg = StackLegend(all);
    if (!leg) return;
      
    TString   txt      = res->GetTitle();
    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
    if      (ib == 0) txt.Append(" "); // (#times1)");
    // else if (ib == 1) txt.Append(" (#times10)");
    else              txt.Append(Form(" (+%d)", off));

    TObject* dummy = 0;
    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
    e->SetMarkerStyle(res->GetMarkerStyle());
    e->SetMarkerSize(res->GetMarkerSize());
    e->SetMarkerColor(kBlack);
    e->SetFillColor(0);
    e->SetFillStyle(0);
    e->SetLineColor(kBlack);
  }

  /** 
   * Get or create a stack legend.  This is done by adding a TLegend
   * object to the list of functions for the first histogram in the
   * stack.
   * 
   * @param stack Stack to get the legend from/modify 
   * 
   * @return The legend object or null
   */
  TLegend* StackLegend(THStack* stack) 
  {
    TList* hists = stack->GetHists();
    if (!hists) return 0;
    
    TObject* first = hists->First();
    if (!first) return 0;

    TH1*    hist = static_cast<TH1*>(first);
    TList*  list = hist->GetListOfFunctions();
    TObject* o   = list->FindObject("legend");
    if (o) return static_cast<TLegend*>(o);
    
    TLegend* l = new TLegend(0.65, 0.65, 0.9, 0.9, "", "NDC");
    l->SetName("legend");
    l->SetBorderSize(0);
    l->SetFillColor(0);
    l->SetFillStyle(0);
    l->SetTextFont(42);
    list->Add(l);
    
    return l;
  }
  
  /* =================================================================
   *
   * Measurements from other sources, such as published ALICE, CMS, ...
   */
  TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
			      Int_t factor)
  {
    TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);", 
                      type, eta, sNN);
    TGraphAsymmErrors* g = 
      reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
    if (!g) { 
      // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
      //         type, eta, sNN);
      return 0;
    }
  

    for (Int_t j = 0; j < g->GetN(); j++) { 
      g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
      g->SetPointError(j, g->GetEXlow()[j], g->GetEXhigh()[j], 
		       g->GetEYlow()[j]*factor, g->GetEYhigh()[j]*factor);
    }
    return g;
  }    
};

void
UnfoldMultDists(const TString& method="Bayes", 
		Double_t       regParam=1e-30,
		const TString& measured="forward_multdists.root", 
		const TString& corrections="")
{
  Unfolder m;
  m.Run(measured, corrections, method, regParam);
}