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