UnfoldMultDists.C

Go to the documentation of this file.

#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>

struct Unfolder
{
  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
  };

  Unfolder() {}
  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;
  }
  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;
  }
  static TCollection* GetCollection(TCollection*   c, 
            const TString& name, 
            Bool_t         verbose=-true)
  {
    return static_cast<TCollection*>(GetObject(c, name, 
                 TCollection::Class(),
                 verbose));
  }
  static TH1* GetH1(TCollection* c, const TString& name, Bool_t verbose=true) 
  {
    return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verbose));
  }
  static TH2* GetH2(TCollection* c, const TString& name, Bool_t verbose=true) 
  {
    return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verbose));
  }
  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());
  }
  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());
  }
  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());
  }
  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;
  }
  
  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_PHYSICS)/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();
  }
  static void AppendAnd(TString& trg, const TString& what)
  {
    if (!trg.IsNull()) trg.Append(" & ");
    trg.Append(what);
  }
  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;
    enum { 
      kInel        = 0x0001, 
      kInelGt0     = 0x0002, 
      kNSD         = 0x0004, 
      kEmpty       = 0x0008, 
      kA           = 0x0010, 
      kB           = 0x0020, 
      kC           = 0x0080,  
      kE           = 0x0100,
      kPileUp      = 0x0200,    
      kMCNSD       = 0x0400,    
      kOffline     = 0x0800,
      kNClusterGt0 = 0x1000,
      kV0AND       = 0x2000, 
      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();
  }
  void SaveSummarize()
  {
    std::ofstream f("SummarizeUnfold.C");
    f << "void SummarizeUnfold(const char* file=\"forward_unfolded.root\")\n"
      << "{\n"
      << "  const char* fwd=\"$ALICE_PHYSICS/PWGLF/FORWARD/analysis2\";\n"
      << "  gROOT->LoadMacro(Form(\"%s/DrawUnfoldedSummary.C\",fwd));\n"
      << "  DrawUnfoldedSummary(file);\n"
      << "}\n"
      << "// EOF" << std::endl;
    f.close();
  }
  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;
    }
  }
  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(0);

    // 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);
  }
  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);
    }
  }
  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) > 1e-3) {
      // Not symmetric bin 
      // Info("Other2Stack", "bin [%f,%f] is not symmetric (%f)",
      //      eta1, eta2, TMath::Abs(eta2-eta1));
      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;
      }
    }
  }
  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);
  }

  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);
}