PWGLF/FORWARD/analysis2/scripts/UnfoldMultDists.C

   1 /**
   2  * @file   UnfoldMultDists.C
   3  * @author Christian Holm Christensen <cholm@nbi.dk>
   4  * @date   Tue Nov 12 09:25:52 2013
   5  *
   6  * @brief  A class to do unfolding
   7  *
   8  *
   9  * @ingroup pwglf_forward_multdist
  10  */
  11 #include <TFile.h>
  12 #include <TList.h>
  13 #include <TH1.h>
  14 #include <TH2.h>
  15 #include <THStack.h>
  16 #include <TLegend.h>
  17 #include <TLegendEntry.h>
  18 #include <TClass.h>
  19 #include <TRegexp.h>
  20 #include <TMath.h>
  21 #include <TParameter.h>
  22 #include <TMultiGraph.h>
  23 #include <TGraphAsymmErrors.h>
  24 #include "RooUnfold.h"
  25 #include "RooUnfoldResponse.h"
  26 #include <fstream>
  27
  28 /**
  29  * Class to do unfolding of raw histograms produced by AliForwardMultDists
  30  *
  31  * @ingroup pwglf_forward_multdist
  32  */
  33 struct Unfolder
  34 {
  35   /**
  36    * Colours used
  37    *
  38    */
  39   enum {
  40     kColorMeasured = kOrange-2,
  41     kColorTruth    = kBlue-3,
  42     kColorAccepted = kMagenta-3,
  43     kColorTrgVtx   = kBlack,
  44     kColorUnfolded = kOrange+2,
  45     kColorCorrected= kRed+2,
  46     kColorError    = kBlue-10,
  47     kColorALICE    = kPink+1,
  48     kColorCMS      = kGreen+2
  49   };
  50
  51   /**
  52    * Constructor
  53    */
  54   Unfolder() {}
  55   /**
  56    * Get a top collection from a file
  57    *
  58    * @param fileName Name of file
  59    * @param results  Wheter it's the results collection or not
  60    *
  61    * @return Collection or null
  62    */
  63   static TCollection* GetTop(const TString& fileName, Bool_t results=false)
  64   {
  65     TFile* file = TFile::Open(fileName, "READ");
  66     if (!file) {
  67       Error("GetTop", "Failed to open %s", fileName.Data());
  68       return 0;
  69     }
  70     TCollection* ret = 0;
  71     TString cName(Form("ForwardMult%s", results ? "Results" : "Sums"));
  72     file->GetObject(cName, ret);
  73     if (!ret)
  74       Error("GetTop", "Failed to get collection %s from %s",
  75             cName.Data(), fileName.Data());
  76     file->Close();
  77     return ret;
  78   }
  79   /**
  80    * Get an object from a collection
  81    *
  82    * @param c    Collection
  83    * @param name Name of object
  84    * @param cl   Possible class to check against
  85    * @param verbose  Be verbose
  86    *
  87    * @return Pointer to object or null
  88    */
  89   static TObject* GetObject(TCollection* c, const TString& name,
  90                             TClass* cl, Bool_t verbose=true)
  91   {
  92     TObject* o = c->FindObject(name);
  93     if (!o) {
  94       if (verbose)
  95         Warning("GetObject", "%s not found in %s", name.Data(), c->GetName());
  96       return 0;
  97     }
  98     if (cl && !o->IsA()->InheritsFrom(cl)) {
  99       if (verbose)
 100         Warning("GetCollection", "%s is not a %s but a %s",
 101                 name.Data(), cl->GetName(), o->ClassName());
 102       return 0;
 103     }
 104     return o;
 105   }
 106   /**
 107    * Get a collection
 108    *
 109    * @param c        Parent collection
 110    * @param name     Name of object to findf
 111    * @param verbose  Be verbose
 112    *
 113    * @return
 114    */
 115   static TCollection* GetCollection(TCollection*   c,
 116                                     const TString& name,
 117                                     Bool_t         verbose=-true)
 118   {
 119     return static_cast<TCollection*>(GetObject(c, name,
 120                                                TCollection::Class(),
 121                                                verbose));
 122   }
 123   /**
 124    * Get a 1D histogram from a collection
 125    *
 126    * @param c    Collection
 127    * @param name Nanme of histogram
 128    * @param verbose  Be verbose
 129    *
 130    * @return Pointer to object or null
 131    */
 132   static TH1* GetH1(TCollection* c, const TString& name, Bool_t verbose=true)
 133   {
 134     return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verbose));
 135   }
 136   /**
 137    * Get a 2D histogram from a collection
 138    *
 139    * @param c    Collection
 140    * @param name Nanme of histogram
 141    * @param verbose  Be verbose
 142    *
 143    * @return Pointer to object or null
 144    */
 145   static TH2* GetH2(TCollection* c, const TString& name, Bool_t verbose=true)
 146   {
 147     return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verbose));
 148   }
 149   /**
 150    * Get an unsigned short parameter from the collection
 151    *
 152    * @param c    Collection
 153    * @param name Parameter name
 154    * @param v    Value
 155    *
 156    * @return Value
 157    */
 158   static void GetParameter(TCollection* c, const TString& name, UShort_t& v)
 159   {
 160     TObject* o = GetObject(c, name, TParameter<int>::Class(), true);
 161     v = (!o ? 0 : o->GetUniqueID());
 162   }
 163   /**
 164    * Get an unsigned short parameter from the collection
 165    *
 166    * @param c    Collection
 167    * @param name Parameter name
 168    * @param v    Value
 169    *
 170    * @return Value
 171    */
 172   static void GetParameter(TCollection* c, const TString& name, ULong_t& v)
 173   {
 174     TObject* o = GetObject(c, name, TParameter<long>::Class(), true);
 175     v = (!o ? 0 : o->GetUniqueID());
 176   }
 177   /**
 178    * Get an unsigned short parameter from the collection
 179    *
 180    * @param c    Collection
 181    * @param name Parameter name
 182    * @param v    Value
 183    *
 184    * @return Value
 185    */
 186   static void GetParameter(TCollection* c, const TString& name, Double_t& v)
 187   {
 188     TObject* o = GetObject(c, name, TParameter<double>::Class(), true);
 189     v = (!o ? 0 : static_cast<TParameter<double>*>(o)->GetVal());
 190   }
 191   /**
 192    * Get the method identifier
 193    *
 194    * @param method Method
 195    *
 196    * @return Method identifier
 197    */
 198   static UInt_t MethodId(TString& method)
 199   {
 200     struct Method {
 201       UInt_t  id;
 202       TString name;
 203     };
 204     const Method methods[] = { {RooUnfold::kNone,    "None"},
 205                                {RooUnfold::kBayes,   "Bayes"},
 206                                {RooUnfold::kSVD,     "SVD"},
 207                                {RooUnfold::kBinByBin,"BinByBin"},
 208                                {RooUnfold::kTUnfold, "TUnfold"},
 209                                {RooUnfold::kInvert,  "Invert"},
 210                                {RooUnfold::kDagostini,"Dagostini"},
 211                                {0xDeadBeef,           "unknown"} };
 212     const Method* pMethod = methods;
 213     while (pMethod->id != 0xDeadBeef) {
 214       if (method.BeginsWith(pMethod->name, TString::kIgnoreCase)) {
 215         method = pMethod->name;
 216         break;
 217       }
 218       pMethod++;
 219     }
 220     if (pMethod->id == 0xDeadBeef)
 221       Error("MethodId", "Unknown unfolding method: %s", method.Data());
 222
 223     return pMethod->id;
 224   }
 225
 226   /**
 227    * Run the unfolding and correction task.
 228    *
 229    * The @a measuredFile is assumed to have the structure
 230    *
 231    * @verbatim
 232    *     /- ForwardMultSums         (TCollection)
 233    *          |- [type]             (TCollection)
 234    *          |    |- [bin]         (TCollection)
 235    *          |    |    `- rawDist  (TH1)
 236    *          |    |- [bin]
 237    *          |    ...
 238    *          |- [type]
 239    *          ...
 240    * @endverbatim
 241    *
 242    * and @a corrFile is assumed to have the structure
 243    *
 244    * @verbatim
 245    *     /- ForwardMultResults            (TCollection)
 246    *          |- [type]                   (TCollection)
 247    *          |    |- [bin]               (TCollection)
 248    *          |    |    |- truth          (TH1)
 249    *          |    |    |- truthAccepted  (TH1)
 250    *          |    |    |- triggerVertex  (TH1)
 251    *          |    |    `- response       (TH2)
 252    *          |    |- [bin]
 253    *          |    ...
 254    *          |- [type]
 255    *          ...
 256    * @endverbatim
 257    *
 258    * where @c [type] is one of <i>symmetric</i>, <i>positive</i>,
 259    * <i>negative</i>, or <i>other</i>, and [bin] is the @f$ \eta@f$
 260    * bin named like
 261    *
 262    * @verbatim
 263    *   [bin]          := [eta_spec] _ [eta_spec]
 264    *   [eta_spec]     := [sign_char] [integer_part] d [decimal_part]
 265    *   [sign_part]    := p          positive eta
 266    *                  |  m          negative eta
 267    *   [integer_part] := [number]
 268    *   [decimal_part] := [number] [number]
 269    *   [number]       := 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
 270    * @endverbatim
 271    *
 272    * That is, the bin @f$ -3\le\eta\ge3@f$ is labeled
 273    * <b>m3d00_p3d00</b>, @f$ 0\le\eta\ge2.5@f$ is <b>p0d00_p2d50</b>
 274    *
 275    * @a measuredFile and @a corrFile can point to the same file.  If
 276    * @a corrFile is not specified, it is assumed that @a measuredFile
 277    * has the expected @a corrFile @e in @e addition to the
 278    * expected content of that file.
 279    *
 280    * @param measuredFile Name of file containing measured data
 281    * @param corrFile     Name of file containing correction data
 282    * @param method       Unfolding method to use
 283    * @param regParam     Regularization parameter
 284    */
 285   void Run(const TString& measuredFile, const TString& corrFile,
 286            const TString& method="Bayes", Double_t regParam=4)
 287   {
 288     // Get the input collections
 289     if (measuredFile.IsNull()) {
 290       Error("Run", "No measurements given");
 291       return;
 292     }
 293     TCollection* mTop = GetTop(measuredFile, false);
 294     TCollection* cTop = GetTop(corrFile.IsNull() ? measuredFile : corrFile,
 295                                true);
 296     if (!mTop || !cTop) return;
 297
 298     // Get some info from the input collection
 299     UShort_t sys;
 300     UShort_t sNN;
 301     ULong_t  trig;
 302     Double_t minZ;
 303     Double_t maxZ;
 304     GetParameter(mTop, "sys",     sys);
 305     GetParameter(mTop, "sNN",     sNN);
 306     GetParameter(mTop, "trigger", trig);
 307     GetParameter(mTop, "minIpZ",  minZ);
 308     GetParameter(mTop, "maxIpZ",  maxZ);
 309     if (sys == 0 || sNN == 0)
 310       Warning("Run", "System (%d) and/or collision energy (%d) unknown",
 311               sys, sNN);
 312
 313     // Open the output file
 314     TFile* out = TFile::Open("forward_unfolded.root", "RECREATE");
 315     if (!out) {
 316       Error("Run", "Failed to open output file");
 317       return;
 318     }
 319
 320     // Decode method option and store in file
 321     TString meth(method);
 322     UInt_t  mId = MethodId(meth);
 323     if (mId == 0xDeadBeef) return;
 324
 325     // Store information
 326     SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
 327                     corrFile.IsNull());
 328
 329     // Load other data
 330     TString savPath(gROOT->GetMacroPath());
 331     gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2/scripts",
 332                              gROOT->GetMacroPath()));
 333     // Always recompile
 334     if (!gROOT->GetClass("OtherPNch"))
 335       gROOT->LoadMacro("OtherPNchData.C++");
 336     gROOT->SetMacroPath(savPath);
 337
 338     // Loop over the input
 339     const char*  inputs[] = { "symmetric", "positive", "negative", 0 };
 340     const char** pinput   = inputs;
 341     while (*pinput) {
 342       TCollection* mInput = GetCollection(mTop, *pinput, false);
 343       TCollection* cInput = GetCollection(cTop, *pinput, false);
 344       if (mInput && cInput)
 345         ProcessType(mInput, cInput, mId, regParam, out,
 346                     sys, sNN);
 347       pinput++;
 348     }
 349
 350     out->Write();
 351     // out->ls();
 352     out->Close();
 353
 354     SaveSummarize();
 355   }
 356   /**
 357    * Append an & to a string and the next term.
 358    *
 359    * @param trg  Output string
 360    * @param what Term
 361    */
 362   static void AppendAnd(TString& trg, const TString& what)
 363   {
 364     if (!trg.IsNull()) trg.Append(" & ");
 365     trg.Append(what);
 366   }
 367   /**
 368    * Store some information on the output
 369    *
 370    * @param dir      Where to store
 371    * @param method   Method used
 372    * @param mId      Method identifier
 373    * @param regParam Regularization parameter
 374    * @param sys      Collision system
 375    * @param sNN      Center of mass energy
 376    * @param trigger  Trigger mask
 377    * @param minIpZ   Least z coordinate of interaction point
 378    * @param maxIpZ   Largest z coordinate of interaction point
 379    * @param self     Self-consistency check
 380    */
 381   void SaveInformation(TDirectory* dir,
 382                        const TString& method,
 383                        Int_t          mId,
 384                        Double_t       regParam,
 385                        UShort_t       sys,
 386                        UShort_t       sNN,
 387                        UInt_t         trigger,
 388                        Double_t       minIpZ,
 389                        Double_t       maxIpZ,
 390                        Bool_t         self) const
 391   {
 392     dir->cd();
 393
 394     TParameter<bool>* pM = new TParameter<bool>("self", self);
 395     pM->SetBit(BIT(19));
 396     pM->Write();
 397
 398     TNamed* outMeth = new TNamed("method", method.Data());
 399     outMeth->SetUniqueID(mId);
 400     outMeth->Write();
 401
 402     TParameter<double>* pR = new TParameter<double>("regParam", regParam);
 403     pR->SetBit(BIT(19));
 404     pR->Write();
 405
 406     TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
 407     TNamed* pS = new TNamed("sys", tS.Data()); pS->SetUniqueID(sys);
 408     pS->Write();
 409
 410     TString tE;
 411     if      (sNN <  1000) tE = Form("%dGeV", sNN);
 412     else if (sNN <  3000) tE = Form("%4.2fTeV", float(sNN)/1000);
 413     else                  tE = Form("%dTeV", sNN/1000);
 414     TNamed* pE = new TNamed("sNN", tE.Data()); pE->SetUniqueID(sNN);
 415     pE->Write();
 416
 417     TString tT;
 418     /**
 419      * Bits of the trigger pattern
 420      */
 421     enum {
 422       /** In-elastic collision */
 423       kInel        = 0x0001,
 424       /** In-elastic collision with at least one SPD tracklet */
 425       kInelGt0     = 0x0002,
 426       /** Non-single diffractive collision */
 427       kNSD         = 0x0004,
 428       /** Empty bunch crossing */
 429       kEmpty       = 0x0008,
 430       /** A-side trigger */
 431       kA           = 0x0010,
 432       /** B(arrel) trigger */
 433       kB           = 0x0020,
 434       /** C-side trigger */
 435       kC           = 0x0080,
 436       /** Empty trigger */
 437       kE           = 0x0100,
 438       /** pileup from SPD */
 439       kPileUp      = 0x0200,
 440       /** true NSD from MC */
 441       kMCNSD       = 0x0400,
 442       /** Offline MB triggered */
 443       kOffline     = 0x0800,
 444       /** At least one SPD cluster */
 445       kNClusterGt0 = 0x1000,
 446       /** V0-AND trigger */
 447       kV0AND       = 0x2000,
 448       /** Satellite event */
 449       kSatellite   = 0x4000
 450     };
 451     if ((trigger & kInel)        != 0x0) AppendAnd(tT, "INEL");
 452     if ((trigger & kInelGt0)     != 0x0) AppendAnd(tT, "INEL>0");
 453     if ((trigger & kNSD)         != 0x0) AppendAnd(tT, "NSD");
 454     if ((trigger & kV0AND)       != 0x0) AppendAnd(tT, "V0AND");
 455     if ((trigger & kA)           != 0x0) AppendAnd(tT, "A");
 456     if ((trigger & kB)           != 0x0) AppendAnd(tT, "B");
 457     if ((trigger & kC)           != 0x0) AppendAnd(tT, "C");
 458     if ((trigger & kE)           != 0x0) AppendAnd(tT, "E");
 459     if ((trigger & kMCNSD)       != 0x0) AppendAnd(tT, "MCNSD");
 460     if ((trigger & kNClusterGt0) != 0x0) AppendAnd(tT, "NCluster>0");
 461     if ((trigger & kSatellite)   != 0x0) AppendAnd(tT, "Satellite");
 462     TNamed* pT = new TNamed("trigger", tT.Data()); pT->SetUniqueID(trigger);
 463     pT->Write();
 464
 465     TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
 466     pY->SetBit(BIT(19));
 467     pY->Write();
 468
 469     TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
 470     pZ->SetBit(BIT(19));
 471     pZ->Write();
 472   }
 473   /**
 474    * Save a script to do a summary of this step
 475    *
 476    */
 477   void SaveSummarize()
 478   {
 479     std::ofstream f("SummarizeUnfold.C");
 480     f << "void SummarizeUnfold(const char* file=\"forward_unfolded.root\")\n"
 481       << "{\n"
 482       << "  const char* fwd=\"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n"
 483       << "  gROOT->LoadMacro(Form(\"%s/DrawUnfoldedSummary.C\",fwd));\n"
 484       << "  DrawUnfoldedSummary(file);\n"
 485       << "}\n"
 486       << "// EOF" << std::endl;
 487     f.close();
 488   }
 489   /**
 490    * Process a single type - i.e., one of <i>symmetric</i>,
 491    * <i>positive</i>, <i>negative</i>, or <i>other</i> - by looping
 492    * over all contained objects and call ProcessBin for each found
 493    * bin.
 494    *
 495    * @param measured     Input collection of measured data
 496    * @param corrections  Input collection of correction data
 497    * @param method       Unfolding method to use
 498    * @param regParam     Regularisation parameter
 499    * @param out          Output directory.
 500    * @param sys          Collision system
 501    * @param sNN          Collision energy
 502    */
 503   void ProcessType(TCollection* measured,
 504                    TCollection* corrections,
 505                    UInt_t       method,
 506                    Double_t     regParam,
 507                    TDirectory*  out,
 508                    UShort_t     sys,
 509                    UShort_t     sNN)
 510   {
 511     Printf("  Processing %s ...", measured->GetName());
 512     TDirectory* dir = out->mkdir(measured->GetName());
 513
 514     // Make some summary stacks
 515     THStack*  allMeasured  = new THStack("measured",
 516                                          "Measured P(#it{N}_{ch})");
 517     THStack*  allTruth     = new THStack("truth",
 518                                          "MC 'truth' P(#it{N}_{ch})");
 519     THStack*  allTruthA    = new THStack("truthAccepted",
 520                                          "MC 'truth' accepted P(#it{N}_{ch})");
 521     THStack*  allUnfolded  = new THStack("unfolded",
 522                                          "Unfolded P(#it{N}_{ch})");
 523     THStack*  allCorrected = new THStack("corrected",
 524                                          "Corrected P(#it{N}_{ch})");
 525     THStack*  allRatio     = (sys != 1 ? 0 :
 526                               new THStack("ratios", "Ratios to other"));
 527     TMultiGraph* allALICE  = (sys != 1 ? 0 :
 528                               new TMultiGraph("alice", "ALICE Published"));
 529     TMultiGraph* allCMS    = (sys != 1 ? 0 :
 530                               new TMultiGraph("cms", "CMS Published"));
 531
 532     // Loop over the list of objects.
 533     static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
 534     TIter          next(measured);
 535     TObject*       o = 0;
 536     Int_t          i = 0;
 537     Double_t       r = regParam;
 538     while ((o = next())) {
 539       // Go back to where we where
 540       dir->cd();
 541
 542       // if not a collection, don't bother
 543       if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
 544
 545       // If it doesn't match our regular expression, don't bother
 546       TString n(o->GetName());
 547       if (n.Index(regex) == kNPOS) {
 548         // Warning("ScanType", "%s in %s doesn't match eta range regexp",
 549         //         n.Data(), real->GetName());
 550         continue;
 551       }
 552       TCollection* mBin = static_cast<TCollection*>(o);
 553       TCollection* cBin = GetCollection(corrections, n.Data());
 554       if (!cBin) continue;
 555
 556       THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
 557       if (!binS) continue;
 558
 559       TH1* result = 0;
 560       Bin2Stack(binS, i, allMeasured, allTruth, allTruthA,
 561                 allUnfolded, allCorrected, result);
 562
 563       TGraph* alice = 0;
 564       TGraph* cms   = 0;
 565       Other2Stack(o->GetName(), i, sNN, allALICE, allCMS, alice, cms);
 566       Ratio2Stack(i, result, alice, cms, allRatio);
 567       i++;
 568     }
 569     dir->Add(allMeasured);
 570     dir->Add(allTruth);
 571     dir->Add(allTruthA);
 572     dir->Add(allUnfolded);
 573     dir->Add(allCorrected);
 574     if (allALICE && allALICE->GetListOfGraphs()) {
 575       if (allALICE->GetListOfGraphs()->GetEntries() > 0)
 576         dir->Add(allALICE);
 577       else
 578         delete allALICE;
 579     }
 580     if (allCMS && allCMS->GetListOfGraphs()) {
 581       if (allCMS->GetListOfGraphs()->GetEntries() > 0)
 582         dir->Add(allCMS);
 583       else
 584         delete allCMS;
 585     }
 586     if (allRatio && allRatio->GetHists()) {
 587       if (allRatio->GetHists()->GetEntries() > 0)
 588         dir->Add(allRatio);
 589       else
 590         delete allRatio;
 591     }
 592   }
 593   /**
 594    * Process a single eta bin
 595    *
 596    * @param measured     Input collection of measured data
 597    * @param corrections  Input collection of correction data
 598    * @param method       Unfolding method to use
 599    * @param regParam     Regularisation parameter
 600    * @param out          Output directory.
 601    *
 602    * @return Stack of histograms or null
 603    */
 604   THStack* ProcessBin(TCollection* measured,
 605                       TCollection* corrections,
 606                       UInt_t       method,
 607                       Double_t     regParam,
 608                       TDirectory*  out)
 609   {
 610     Printf("   Processing %s ...", measured->GetName());
 611     // Try to get the data
 612     TH1* inRaw    = GetH1(measured,    "rawDist");
 613     TH1* inTruth  = GetH1(corrections, "truth");
 614     TH1* inTruthA = GetH1(corrections, "truthAccepted");
 615     TH1* inTrgVtx = GetH1(corrections, "triggerVertex");
 616     TH2* inResp   = GetH2(corrections, "response");
 617     if (!inRaw || !inTruth || !inTruthA || !inTrgVtx || !inResp)
 618       return 0;
 619
 620     // Make output directory
 621     TDirectory* dir = out->mkdir(measured->GetName());
 622     dir->cd();
 623
 624     // Copy the input to the output
 625     TH1* outRaw    = static_cast<TH1*>(inRaw    ->Clone("measured"));
 626     TH1* outTruth  = static_cast<TH1*>(inTruth  ->Clone("truth"));
 627     TH1* outTruthA = static_cast<TH1*>(inTruthA ->Clone("truthAccepted"));
 628     TH1* outTrgVtx = static_cast<TH1*>(inTrgVtx ->Clone("triggerVertex"));
 629     TH2* outResp   = static_cast<TH2*>(inResp   ->Clone("response"));
 630
 631     // Make our response matrix
 632     RooUnfoldResponse matrix(0, 0, inResp);
 633
 634     // Store regularization parameter
 635     Double_t             r        = regParam;
 636     RooUnfold::Algorithm algo     = (RooUnfold::Algorithm)method;
 637     RooUnfold*           unfolder = RooUnfold::New(algo, &matrix, inRaw, r);
 638     unfolder->SetVerbose(0);
 639
 640     // Do the unfolding and get the result
 641     TH1* res = unfolder->Hreco();
 642     res->SetDirectory(0);
 643
 644     // Make a copy to store on the output
 645     TH1* outUnfold = static_cast<TH1*>(res->Clone("unfolded"));
 646     TString tit(outUnfold->GetTitle());
 647     tit.ReplaceAll("Unfold Reponse matrix", "Unfolded P(#it{N}_{ch})");
 648     outUnfold->SetTitle(tit);
 649
 650     // Clone the unfolded results and divide by the trigger/vertex
 651     // bias correction
 652     TH1* outCorr   = static_cast<TH1*>(outUnfold->Clone("corrected"));
 653     outCorr->Divide(inTrgVtx);
 654     tit.ReplaceAll("Unfolded", "Corrected");
 655     outCorr->SetTitle(tit);
 656
 657     // Now normalize the output to integral=1
 658     TH1*  hists[] = { outRaw, outUnfold, outCorr, 0 };
 659     TH1** phist   = hists;
 660     while (*phist) {
 661       TH1* h = *phist;
 662       if (h) {
 663         Double_t intg = h->Integral(1, h->GetXaxis()->GetXmax());
 664         h->Scale(1. / intg, "width");
 665       }
 666       phist++;
 667     }
 668
 669     // And make ratios
 670     TH1* ratioTrue = static_cast<TH1*>(outCorr->Clone("ratioCorrTruth"));
 671     tit = ratioTrue->GetTitle();
 672     tit.ReplaceAll("Corrected", "Corrected/MC 'truth'");
 673     ratioTrue->SetTitle(tit);
 674     ratioTrue->Divide(outTruth);
 675     ratioTrue->SetYTitle("P_{corrected}(#it{N}_{ch})/P_{truth}(#it{N}_{ch})");
 676
 677     TH1* ratioAcc  = static_cast<TH1*>(outUnfold->Clone("ratioUnfAcc"));
 678     tit = ratioAcc->GetTitle();
 679     tit.ReplaceAll("Unfolded", "Unfolded/MC selected");
 680     ratioAcc->SetTitle(tit);
 681     ratioAcc->Divide(outTruthA);
 682     ratioAcc->SetYTitle("P_{unfolded}(#it{N}_{ch})/P_{MC}(#it{N}_{ch})");
 683
 684
 685     // Make a stack
 686     tit = measured->GetName();
 687     tit.ReplaceAll("m", "-");
 688     tit.ReplaceAll("p", "+");
 689     tit.ReplaceAll("d", ".");
 690     tit.ReplaceAll("_", "<#it{#eta}<");
 691     THStack* stack = new THStack("all", tit);
 692     stack->Add(outTruth,  "E2");
 693     stack->Add(outTruthA, "E2");
 694     stack->Add(outRaw,    "E1");
 695     stack->Add(outUnfold, "E1");
 696     stack->Add(outCorr,   "E1");
 697     dir->Add(stack);
 698
 699     // Rest of the function is devoted to making the output look nice
 700     outRaw   ->SetDirectory(dir);
 701     outTruth ->SetDirectory(dir);
 702     outTruthA->SetDirectory(dir);
 703     outTrgVtx->SetDirectory(dir);
 704     outResp  ->SetDirectory(dir);
 705     outUnfold->SetDirectory(dir);
 706     outCorr  ->SetDirectory(dir);
 707
 708     outRaw   ->SetMarkerStyle(20);  // Measured is closed
 709     outTruth ->SetMarkerStyle(24);  // MC is open
 710     outTruthA->SetMarkerStyle(24);  // MC is open
 711     outTrgVtx->SetMarkerStyle(20);  // Derived is closed
 712     outUnfold->SetMarkerStyle(20);  // Derived is closed
 713     outCorr  ->SetMarkerStyle(20);  // Derived is closed
 714
 715     outRaw   ->SetMarkerSize(0.9);
 716     outTruth ->SetMarkerSize(1.6);
 717     outTruthA->SetMarkerSize(1.4);
 718     outTrgVtx->SetMarkerSize(1.0);
 719     outUnfold->SetMarkerSize(0.9);
 720     outCorr  ->SetMarkerSize(1.0);
 721
 722     outRaw   ->SetMarkerColor(kColorMeasured);
 723     outTruth ->SetMarkerColor(kColorTruth);
 724     outTruthA->SetMarkerColor(kColorAccepted);
 725     outTrgVtx->SetMarkerColor(kColorTrgVtx);
 726     outUnfold->SetMarkerColor(kColorUnfolded);
 727     outCorr  ->SetMarkerColor(kColorCorrected);
 728
 729     outRaw   ->SetFillColor(kColorError);
 730     outTruth ->SetFillColor(kColorError);
 731     outTruthA->SetFillColor(kColorError);
 732     outTrgVtx->SetFillColor(kColorError);
 733     outUnfold->SetFillColor(kColorError);
 734     outCorr  ->SetFillColor(kColorError);
 735
 736     outRaw   ->SetFillStyle(0);
 737     outTruth ->SetFillStyle(1001);
 738     outTruthA->SetFillStyle(1001);
 739     outTrgVtx->SetFillStyle(0);
 740     outUnfold->SetFillStyle(0);
 741     outCorr  ->SetFillStyle(0);
 742
 743     outRaw   ->SetLineColor(kBlack);
 744     outTruth ->SetLineColor(kBlack);
 745     outTruthA->SetLineColor(kBlack);
 746     outTrgVtx->SetLineColor(kBlack);
 747     outUnfold->SetLineColor(kBlack);
 748     outCorr  ->SetLineColor(kBlack);
 749
 750     // Legend
 751     TLegend* l = StackLegend(stack);
 752     l->AddEntry(outRaw,     "Raw",                 "lp");
 753     l->AddEntry(outTruth,   "MC 'truth'",          "fp");
 754     l->AddEntry(outTruthA,  "MC 'truth' accepted", "fp");
 755     l->AddEntry(outUnfold,  "Unfolded",            "lp");
 756     l->AddEntry(outCorr,    "Corrected",           "lp");
 757
 758     return stack;
 759   }
 760   static void BinAttributes(Int_t i,
 761                             Int_t&    open,
 762                             Int_t&    closed,
 763                             Float_t&  size,
 764                             Double_t& factor)
 765   {
 766     // --- Setup for markers -----------------------------------------
 767     const Int_t   nMarkers = 7;
 768     const Int_t   aClosed[] = { 20,  21,  22,  23,  29,  33,  34  };
 769     const Int_t   aOpen[]   = { 24,  25,  26,  32,  30,  27,  28  };
 770     const Float_t aSize[]   = { 1.1, 1.0, 1.2, 1.2, 1.2, 1.2, 1.0 };
 771     Int_t         j         = i % nMarkers;
 772
 773     open   = aOpen[j];
 774     closed = aClosed[j];
 775     size   = aSize[j];
 776     factor = TMath::Power(10, i);
 777   }
 778   /**
 779    * Add the bin histograms to our summary stacks
 780    *
 781    * @param bin       Bin stack
 782    * @param i         Current off-set in the stacks
 783    * @param measured  All measured @f$ P(N_{ch})@f$
 784    * @param truth     All MC truth @f$ P(N_{ch})@f$
 785    * @param accepted  All MC accepted @f$ P(N_{ch})@f$
 786    * @param unfolded  All unfolded @f$ P(N_{ch})@f$
 787    * @param corrected All corrected @f$ P(N_{ch})@f$
 788    * @param result    The result in this bin
 789    */
 790   void Bin2Stack(const THStack* bin, Int_t i,
 791                  THStack* measured,
 792                  THStack* truth,
 793                  THStack* accepted,
 794                  THStack* unfolded,
 795                  THStack* corrected,
 796                  TH1*&    result)
 797   {
 798     Int_t open, closed;
 799     Double_t factor;
 800     Float_t  size;
 801     BinAttributes(i, open, closed, size, factor);
 802
 803     TIter next(bin->GetHists());
 804     TH1*  h = 0;
 805     while ((h = static_cast<TH1*>(next()))) {
 806       THStack* tmp = 0;
 807       Int_t    col = h->GetMarkerColor();
 808       Int_t    sty = 0;
 809       switch (col) {
 810       case kColorMeasured:  tmp = measured;   sty = closed;  break;
 811       case kColorTruth:     tmp = truth;      sty = open;    break;
 812       case kColorAccepted:  tmp = accepted;   sty = open;    break;
 813       case kColorUnfolded:  tmp = unfolded;   sty = closed;  break;
 814       case kColorCorrected: tmp = corrected;  sty = closed;  break;
 815       default: continue;
 816       }
 817       // Now clone, and add to the appropriate stack
 818       TH1* cln = static_cast<TH1*>(h->Clone(h->GetName()));
 819       cln->SetDirectory(0);
 820       cln->SetMarkerStyle(sty);
 821       cln->SetMarkerSize(size);
 822       cln->Scale(factor); // Scale by 10^i
 823       if (col == kColorCorrected) result = cln;
 824
 825       // Make sure we do not get the old legend
 826       TObject* tst = cln->FindObject("legend");
 827       if (tst) cln->GetListOfFunctions()->Remove(tst);
 828
 829       tmp->Add(cln, next.GetOption());
 830     }
 831
 832     // Add entries to our stacks
 833     TString   txt      = bin->GetTitle();
 834     if      (i == 0) txt.Append(" (#times1)");
 835     else if (i == 1) txt.Append(" (#times10)");
 836     else             txt.Append(Form(" (#times10^{%d})", i));
 837     THStack*  stacks[] = { measured, truth, accepted, unfolded, corrected, 0 };
 838     THStack** pstack   = stacks;
 839     while (*pstack) {
 840       TLegend* leg = StackLegend(*pstack);
 841       pstack++;
 842       if (!leg) continue;
 843
 844       TObject* dummy = 0;
 845       TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
 846       e->SetMarkerStyle(closed);
 847       e->SetMarkerSize(1.2*size);
 848       e->SetMarkerColor(kBlack);
 849       e->SetFillColor(0);
 850       e->SetFillStyle(0);
 851       e->SetLineColor(kBlack);
 852     }
 853   }
 854   /**
 855    * Add distributions from other experiments to stacks
 856    *
 857    * @param name     Name of current bin
 858    * @param i        Index of current bin
 859    * @param sNN      Center of mass energy
 860    * @param allALICE Stack of ALICE data
 861    * @param allCMS   Stack of CMS data
 862    * @param alice    Possible ALICE result on return
 863    * @param cms      Possible CMS result on return
 864    */
 865   void Other2Stack(const TString& name, Int_t i,
 866                    UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
 867                    TGraph*& alice, TGraph*& cms)
 868   {
 869     if (!allALICE && !allCMS) return;
 870
 871     TString tmp(name);
 872     tmp.ReplaceAll("p", "+");
 873     tmp.ReplaceAll("m", "-");
 874     tmp.ReplaceAll("_", " ");
 875     tmp.ReplaceAll("d", ".");
 876     TObjArray* tokens = tmp.Tokenize(" ");
 877     if (!tokens || tokens->GetEntriesFast() < 2) {
 878       Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
 879       if (tokens) tokens->Delete();
 880       return;
 881     }
 882     Double_t eta1 = static_cast<TObjString*>(tokens->At(0))->String().Atof();
 883     Double_t eta2 = static_cast<TObjString*>(tokens->At(1))->String().Atof();
 884     tokens->Delete();
 885
 886     if (TMath::Abs(eta2+eta1) > 1e-3) {
 887       // Not symmetric bin
 888       // Info("Other2Stack", "bin [%f,%f] is not symmetric (%f)",
 889       //      eta1, eta2, TMath::Abs(eta2-eta1));
 890       return;
 891     }
 892     Double_t aEta = TMath::Abs(eta1);
 893
 894     Int_t open, closed;
 895     Double_t factor;
 896     Float_t  size;
 897     BinAttributes(i, open, closed, size, factor);
 898
 899     if (allALICE) {
 900       TGraphAsymmErrors* g = GetOther(1, aEta, sNN, factor);
 901       if (g) {
 902         g->SetMarkerStyle(closed);
 903         g->SetMarkerColor(kColorALICE);
 904         g->SetMarkerSize(size);
 905         allALICE->Add(g, "p same");
 906         alice = g;
 907       }
 908     }
 909     if (allCMS) {
 910       TGraphAsymmErrors* g = GetOther(0, aEta, sNN, factor);
 911       if (g) {
 912         g->SetMarkerStyle(closed);
 913         g->SetMarkerColor(kColorCMS);
 914         g->SetMarkerSize(size);
 915         allCMS->Add(g, "p same");
 916         cms = g;
 917       }
 918     }
 919   }
 920   /**
 921    * Create ratios to other data
 922    *
 923    * @param ib      Bin number
 924    * @param res     Result
 925    * @param alice   ALICE result if any
 926    * @param cms     CMS result if any
 927    * @param all     Stack to add ratio to
 928    */
 929   void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
 930   {
 931     if (!all || !res || !(alice || cms)) return;
 932
 933     Int_t        off  = 5*ib;
 934     TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
 935     TGraph**     pg   = gs;
 936     while (*pg) {
 937       TGraph*     g = *pg;
 938       const char* n = (g == alice ? "ALICE" : "CMS");
 939
 940       TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
 941       TString tit(r->GetTitle());
 942       tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
 943       r->SetTitle(tit);
 944       r->SetMarkerColor(g->GetMarkerColor());
 945       r->SetLineColor(g->GetLineColor());
 946
 947       TObject* tst = r->FindObject("legend");
 948       if (tst) r->GetListOfFunctions()->Remove(tst);
 949
 950       for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
 951         Double_t c = r->GetBinContent(i);
 952         Double_t e = r->GetBinError(i);
 953         Double_t o = g->Eval(r->GetBinCenter(i));
 954         if (o < 1e-12) {
 955           r->SetBinContent(i, 0);
 956           r->SetBinError(i, 0);
 957           continue;
 958         }
 959         r->SetBinContent(i, (c - o) / o + off);
 960         r->SetBinError(i, e / o);
 961       }
 962       all->Add(r);
 963       pg++;
 964     }
 965     TLegend* leg = StackLegend(all);
 966     if (!leg) return;
 967
 968     TString   txt      = res->GetTitle();
 969     txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
 970     if      (ib == 0) txt.Append(" "); // (#times1)");
 971     // else if (ib == 1) txt.Append(" (#times10)");
 972     else              txt.Append(Form(" (+%d)", off));
 973
 974     TObject* dummy = 0;
 975     TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
 976     e->SetMarkerStyle(res->GetMarkerStyle());
 977     e->SetMarkerSize(res->GetMarkerSize());
 978     e->SetMarkerColor(kBlack);
 979     e->SetFillColor(0);
 980     e->SetFillStyle(0);
 981     e->SetLineColor(kBlack);
 982   }
 983
 984   /**
 985    * Get or create a stack legend.  This is done by adding a TLegend
 986    * object to the list of functions for the first histogram in the
 987    * stack.
 988    *
 989    * @param stack Stack to get the legend from/modify
 990    *
 991    * @return The legend object or null
 992    */
 993   TLegend* StackLegend(THStack* stack)
 994   {
 995     TList* hists = stack->GetHists();
 996     if (!hists) return 0;
 997
 998     TObject* first = hists->First();
 999     if (!first) return 0;
1000
1001     TH1*    hist = static_cast<TH1*>(first);
1002     TList*  list = hist->GetListOfFunctions();
1003     TObject* o   = list->FindObject("legend");
1004     if (o) return static_cast<TLegend*>(o);
1005
1006     TLegend* l = new TLegend(0.65, 0.65, 0.9, 0.9, "", "NDC");
1007     l->SetName("legend");
1008     l->SetBorderSize(0);
1009     l->SetFillColor(0);
1010     l->SetFillStyle(0);
1011     l->SetTextFont(42);
1012     list->Add(l);
1013
1014     return l;
1015   }
1016
1017   /* =================================================================
1018    *
1019    * Measurements from other sources, such as published ALICE, CMS, ...
1020    */
1021   TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
1022                               Int_t factor)
1023   {
1024     TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);",
1025                       type, eta, sNN);
1026     TGraphAsymmErrors* g =
1027       reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
1028     if (!g) {
1029       // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
1030       //         type, eta, sNN);
1031       return 0;
1032     }
1033
1034
1035     for (Int_t j = 0; j < g->GetN(); j++) {
1036       g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
1037       g->SetPointError(j, g->GetEXlow()[j], g->GetEXhigh()[j],
1038                        g->GetEYlow()[j]*factor, g->GetEYhigh()[j]*factor);
1039     }
1040     return g;
1041   }
1042 };
1043
1044 void
1045 UnfoldMultDists(const TString& method="Bayes",
1046                 Double_t       regParam=1e-30,
1047                 const TString& measured="forward_multdists.root",
1048                 const TString& corrections="")
1049 {
1050   Unfolder m;
1051   m.Run(measured, corrections, method, regParam);
1052 }