Fix some documentation issues
[u/mrichter/AliRoot.git] / 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(1);
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) > 1e3) 
887       // Not symmetric bin 
888       return;
889     Double_t aEta = TMath::Abs(eta1);
890
891     Int_t open, closed;
892     Double_t factor; 
893     Float_t  size;
894     BinAttributes(i, open, closed, size, factor);
895
896     if (allALICE) {
897       TGraphAsymmErrors* g = GetOther(1, aEta, sNN, factor);
898       if (g) {
899         g->SetMarkerStyle(closed);
900         g->SetMarkerColor(kColorALICE);
901         g->SetMarkerSize(size);
902         allALICE->Add(g, "p same");
903         alice = g;
904       }
905     }
906     if (allCMS) {
907       TGraphAsymmErrors* g = GetOther(0, aEta, sNN, factor);
908       if (g) {
909         g->SetMarkerStyle(closed);
910         g->SetMarkerColor(kColorCMS);
911         g->SetMarkerSize(size);
912         allCMS->Add(g, "p same");
913         cms = g;
914       }
915     }
916   }
917   /** 
918    * Create ratios to other data 
919    * 
920    * @param ib      Bin number  
921    * @param res     Result
922    * @param alice   ALICE result if any
923    * @param cms     CMS result if any
924    * @param all     Stack to add ratio to 
925    */
926   void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
927   {
928     if (!all || !res || !(alice || cms)) return;
929
930     Int_t        off  = 5*ib;
931     TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
932     TGraph**     pg   = gs;
933     while (*pg) { 
934       TGraph*     g = *pg;
935       const char* n = (g == alice ? "ALICE" : "CMS");
936
937       TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
938       TString tit(r->GetTitle());
939       tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
940       r->SetTitle(tit);
941       r->SetMarkerColor(g->GetMarkerColor());
942       r->SetLineColor(g->GetLineColor());
943
944       TObject* tst = r->FindObject("legend");
945       if (tst) r->GetListOfFunctions()->Remove(tst);
946
947       for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
948         Double_t c = r->GetBinContent(i);
949         Double_t e = r->GetBinError(i);
950         Double_t o = g->Eval(r->GetBinCenter(i));
951         if (o < 1e-12) { 
952           r->SetBinContent(i, 0);
953           r->SetBinError(i, 0);
954           continue;
955         }
956         r->SetBinContent(i, (c - o) / o + off);
957         r->SetBinError(i, e / o);
958       }
959       all->Add(r);
960       pg++;
961     }
962     TLegend* leg = StackLegend(all);
963     if (!leg) return;
964       
965     TString   txt      = res->GetTitle();
966     txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
967     if      (ib == 0) txt.Append(" "); // (#times1)");
968     // else if (ib == 1) txt.Append(" (#times10)");
969     else              txt.Append(Form(" (+%d)", off));
970
971     TObject* dummy = 0;
972     TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
973     e->SetMarkerStyle(res->GetMarkerStyle());
974     e->SetMarkerSize(res->GetMarkerSize());
975     e->SetMarkerColor(kBlack);
976     e->SetFillColor(0);
977     e->SetFillStyle(0);
978     e->SetLineColor(kBlack);
979   }
980
981   /** 
982    * Get or create a stack legend.  This is done by adding a TLegend
983    * object to the list of functions for the first histogram in the
984    * stack.
985    * 
986    * @param stack Stack to get the legend from/modify 
987    * 
988    * @return The legend object or null
989    */
990   TLegend* StackLegend(THStack* stack) 
991   {
992     TList* hists = stack->GetHists();
993     if (!hists) return 0;
994     
995     TObject* first = hists->First();
996     if (!first) return 0;
997
998     TH1*    hist = static_cast<TH1*>(first);
999     TList*  list = hist->GetListOfFunctions();
1000     TObject* o   = list->FindObject("legend");
1001     if (o) return static_cast<TLegend*>(o);
1002     
1003     TLegend* l = new TLegend(0.65, 0.65, 0.9, 0.9, "", "NDC");
1004     l->SetName("legend");
1005     l->SetBorderSize(0);
1006     l->SetFillColor(0);
1007     l->SetFillStyle(0);
1008     l->SetTextFont(42);
1009     list->Add(l);
1010     
1011     return l;
1012   }
1013   
1014   /* =================================================================
1015    *
1016    * Measurements from other sources, such as published ALICE, CMS, ...
1017    */
1018   TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
1019                               Int_t factor)
1020   {
1021     TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);", 
1022                       type, eta, sNN);
1023     TGraphAsymmErrors* g = 
1024       reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
1025     if (!g) { 
1026       // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
1027       //         type, eta, sNN);
1028       return 0;
1029     }
1030   
1031
1032     for (Int_t j = 0; j < g->GetN(); j++) { 
1033       g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
1034       g->SetPointError(j, g->GetEXlow()[j], g->GetEXhigh()[j], 
1035                        g->GetEYlow()[j]*factor, g->GetEYhigh()[j]*factor);
1036     }
1037     return g;
1038   }    
1039 };
1040
1041 void
1042 UnfoldMultDists(const TString& method="Bayes", 
1043                 Double_t       regParam=1e-30,
1044                 const TString& measured="forward_multdists.root", 
1045                 const TString& corrections="")
1046 {
1047   Unfolder m;
1048   m.Run(measured, corrections, method, regParam);
1049 }