4ceb8853bf218d85f035fea8d08d6c842e3dc0ed
[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     UShort_t maxN; 
284     ULong_t  trig; 
285     Double_t minZ; 
286     Double_t maxZ;
287     GetParameter(mTop, "sys",     sys);   
288     GetParameter(mTop, "snn",     sNN);   
289     GetParameter(mTop, "maxN",    maxN);          
290     GetParameter(mTop, "trigger", trig);
291     GetParameter(mTop, "minIpZ",  minZ); 
292     GetParameter(mTop, "maxIpZ",  maxZ); 
293     if (sys == 0 || sNN == 0) 
294       Warning("Run", "System (%d) and/or collision energy (%d) unknown", 
295               sys, sNN);
296     
297     // Open the output file 
298     TFile* out = TFile::Open("forward_unfolded.root", "RECREATE");
299     if (!out) { 
300       Error("Run", "Failed to open output file");
301       return;
302     }
303
304     // Decode method option and store in file 
305     TString meth(method);
306     UInt_t  mId = MethodId(meth);
307     if (mId == 0xDeadBeef) return;
308
309     // Store information 
310     SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,maxN);
311
312     // Load other data 
313     TString savPath(gROOT->GetMacroPath());
314     gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2/scripts",
315                              gROOT->GetMacroPath()));
316     // Always recompile 
317     if (!gROOT->GetClass("OtherPNch"))
318       gROOT->LoadMacro("OtherPNchData.C++");
319     gROOT->SetMacroPath(savPath);
320
321     // Loop over the input 
322     const char*  inputs[] = { "symmetric", "positive", "negative", 0 };
323     const char** pinput   = inputs;
324     while (*pinput) { 
325       TCollection* mInput = GetCollection(mTop, *pinput, false);
326       TCollection* cInput = GetCollection(cTop, *pinput, false);
327       if (mInput && cInput)
328         ProcessType(mInput, cInput, mId, regParam, out,
329                     sys, sNN);
330       pinput++;
331     }      
332
333     out->Write();
334     // out->ls();
335     out->Close();
336
337     SaveSummarize();
338   }
339   static void AppendAnd(TString& trg, const TString& what)
340   {
341     if (!trg.IsNull()) trg.Append(" & ");
342     trg.Append(what);
343   }
344   /** 
345    * Store some information on the output
346    * 
347    * @param dir      Where to store
348    * @param method   Method used
349    * @param mId      Method identifier 
350    * @param regParam Regularization parameter 
351    * @param sys      Collision system
352    * @param sNN      Center of mass energy 
353    * @param trigger  Trigger mask 
354    * @param minIpZ   Least z coordinate of interaction point
355    * @param maxIpZ   Largest z coordinate of interaction point
356    * @param maxN     Largest @f$N_{ch}@f$ to consider 
357    */
358   void SaveInformation(TDirectory* dir, 
359                        const TString& method,
360                        Int_t          mId,
361                        Double_t       regParam, 
362                        UShort_t       sys, 
363                        UShort_t       sNN,
364                        UInt_t         trigger,
365                        Double_t       minIpZ, 
366                        Double_t       maxIpZ, 
367                        UShort_t       maxN) const
368   {
369     dir->cd();
370
371     TNamed* outMeth = new TNamed("method", method.Data());
372     outMeth->SetUniqueID(mId);
373     outMeth->Write();
374
375     TParameter<double>* pR = new TParameter<double>("regParam", regParam);
376     pR->Write();
377     
378     TString tS = (sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "?");
379     TNamed* pS = new TNamed("sys", tS.Data()); pS->SetUniqueID(sys);
380     pS->Write();
381     
382     TString tE;
383     if      (sNN <  1000) tE = Form("%dGeV", sNN);
384     else if (sNN <  3000) tE = Form("%4.2fTeV", float(sNN)/1000);
385     else                  tE = Form("%dTeV", sNN/1000);
386     TNamed* pE = new TNamed("sNN", tE.Data()); pE->SetUniqueID(sNN);
387     pE->Write();
388     
389     TString tT;
390     /** 
391      * Bits of the trigger pattern
392      */
393     enum { 
394       /** In-elastic collision */
395       kInel        = 0x0001, 
396       /** In-elastic collision with at least one SPD tracklet */
397       kInelGt0     = 0x0002, 
398       /** Non-single diffractive collision */
399       kNSD         = 0x0004, 
400       /** Empty bunch crossing */
401       kEmpty       = 0x0008, 
402       /** A-side trigger */
403       kA           = 0x0010, 
404       /** B(arrel) trigger */
405       kB           = 0x0020, 
406       /** C-side trigger */
407       kC           = 0x0080,  
408       /** Empty trigger */
409       kE           = 0x0100,
410       /** pileup from SPD */
411       kPileUp      = 0x0200,    
412       /** true NSD from MC */
413       kMCNSD       = 0x0400,    
414       /** Offline MB triggered */
415       kOffline     = 0x0800,
416       /** At least one SPD cluster */ 
417       kNClusterGt0 = 0x1000,
418       /** V0-AND trigger */
419       kV0AND       = 0x2000, 
420       /** Satellite event */
421       kSatellite   = 0x4000
422     };
423     if ((trigger & kInel)        != 0x0) AppendAnd(tT, "INEL");
424     if ((trigger & kInelGt0)     != 0x0) AppendAnd(tT, "INEL>0");
425     if ((trigger & kNSD)         != 0x0) AppendAnd(tT, "NSD");
426     if ((trigger & kV0AND)       != 0x0) AppendAnd(tT, "V0AND");
427     if ((trigger & kA)           != 0x0) AppendAnd(tT, "A");
428     if ((trigger & kB)           != 0x0) AppendAnd(tT, "B");
429     if ((trigger & kC)           != 0x0) AppendAnd(tT, "C");
430     if ((trigger & kE)           != 0x0) AppendAnd(tT, "E");
431     if ((trigger & kMCNSD)       != 0x0) AppendAnd(tT, "MCNSD");
432     if ((trigger & kNClusterGt0) != 0x0) AppendAnd(tT, "NCluster>0");
433     if ((trigger & kSatellite)   != 0x0) AppendAnd(tT, "Satellite");
434     TNamed* pT = new TNamed("trigger", tT.Data()); pT->SetUniqueID(trigger);
435     pT->Write();
436     
437     TParameter<double>* pY = new TParameter<double>("minIpZ", minIpZ);
438     pY->Write();
439
440     TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
441     pZ->Write();
442
443     TParameter<int>* pN = new TParameter<int>("maxN", maxN);
444     pN->Write();
445   }
446   /** 
447    * Save a script to do a summary of this step 
448    * 
449    */                  
450   void SaveSummarize()
451   {
452     std::ofstream f("SummarizeUnfold.C");
453     f << "void SummarizeUnfold(const char* file=\"forward_unfolded.root\")\n"
454       << "{\n"
455       << "  const char* fwd=\"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n"
456       << "  gROOT->LoadMacro(Form(\"%s/DrawUnfoldedSummary.C\",fwd));\n"
457       << "  DrawUnfoldedSummary(file);\n"
458       << "}\n"
459       << "// EOF" << std::endl;
460     f.close();
461   }
462   /** 
463    * Process a single type - i.e., one of <i>symmetric</i>,
464    * <i>positive</i>, <i>negative</i>, or <i>other</i> - by looping
465    * over all contained objects and call ProcessBin for each found
466    * bin.
467    * 
468    * @param measured     Input collection of measured data
469    * @param corrections  Input collection of correction data
470    * @param method       Unfolding method to use 
471    * @param regParam     Regularisation parameter
472    * @param out          Output directory. 
473    */
474   void ProcessType(TCollection* measured, 
475                    TCollection* corrections,
476                    UInt_t       method,
477                    Double_t     regParam,
478                    TDirectory*  out,
479                    UShort_t     sys, 
480                    UShort_t     sNN)
481   {
482     Printf("  Processing %s ...", measured->GetName());
483     TDirectory* dir = out->mkdir(measured->GetName());
484     
485     // Make some summary stacks 
486     THStack*  allMeasured  = new THStack("measured",      
487                                          "Measured P(#it{N}_{ch})");
488     THStack*  allTruth     = new THStack("truth",        
489                                          "MC 'truth' P(#it{N}_{ch})");
490     THStack*  allTruthA    = new THStack("truthAccepted",
491                                          "MC 'truth' accepted P(#it{N}_{ch})");
492     THStack*  allUnfolded  = new THStack("unfolded",
493                                          "Unfolded P(#it{N}_{ch})");
494     THStack*  allCorrected = new THStack("corrected",
495                                          "Corrected P(#it{N}_{ch})");
496     TMultiGraph* allALICE  = (sys != 1 ? 0 : 
497                               new TMultiGraph("alice", "ALICE Published"));
498     TMultiGraph* allCMS    = (sys != 1 ? 0 : 
499                               new TMultiGraph("cms", "CMS Published"));
500
501     // Loop over the list of objects. 
502     static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
503     TIter          next(measured);
504     TObject*       o = 0;
505     Int_t          i = 0;
506     Double_t       r = regParam;
507     while ((o = next())) {
508       // Go back to where we where 
509       dir->cd();
510       
511       // if not a collection, don't bother 
512       if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
513     
514       // If it doesn't match our regular expression, don't bother 
515       TString n(o->GetName());
516       if (n.Index(regex) == kNPOS) { 
517         // Warning("ScanType", "%s in %s doesn't match eta range regexp", 
518         //         n.Data(), real->GetName());
519         continue;
520       }
521       TCollection* mBin = static_cast<TCollection*>(o);
522       TCollection* cBin = GetCollection(corrections, n.Data());
523       if (!cBin) continue;
524
525       THStack* binS = ProcessBin(mBin, cBin, method, r, dir);
526       if (!binS) continue;
527
528       Bin2Stack(binS, i, allMeasured, allTruth, allTruthA, 
529                 allUnfolded, allCorrected);
530       Other2Stack(o->GetName(), i, sNN, allALICE, allCMS);
531       i++;
532     }
533     dir->Add(allMeasured);
534     dir->Add(allTruth);
535     dir->Add(allTruthA);
536     dir->Add(allUnfolded);
537     dir->Add(allCorrected);
538     if (allALICE && allALICE->GetListOfGraphs()) {
539       if (allALICE->GetListOfGraphs()->GetEntries() > 0)
540         dir->Add(allALICE);
541       else 
542         delete allALICE;
543     }
544     if (allCMS && allCMS->GetListOfGraphs()) {
545       if (allCMS->GetListOfGraphs()->GetEntries() > 0) 
546         dir->Add(allCMS);
547       else 
548         delete allCMS;
549     }
550   }
551   /** 
552    * Process a single eta bin 
553    * 
554    * @param measured     Input collection of measured data
555    * @param corrections  Input collection of correction data
556    * @param method       Unfolding method to use 
557    * @param regParam     Regularisation parameter
558    * @param out          Output directory. 
559    *
560    * @return Stack of histograms or null 
561    */
562   THStack* ProcessBin(TCollection* measured, 
563                       TCollection* corrections, 
564                       UInt_t       method,
565                       Double_t     regParam, 
566                       TDirectory*  out)
567   {
568     Printf("   Processing %s ...", measured->GetName());
569     // Try to get the data 
570     TH1* inRaw    = GetH1(measured, "rawDist");
571     TH1* inTruth  = GetH1(corrections, "truth");
572     TH1* inTruthA = GetH1(corrections, "truthAccepted");
573     TH1* inTrgVtx = GetH1(corrections, "triggerVertex");
574     TH2* inResp   = GetH2(corrections, "response");
575     if (!inRaw || !inTruth || !inTruthA || !inTrgVtx || !inResp) 
576       return 0;
577     
578     // Make output directory
579     TDirectory* dir = out->mkdir(measured->GetName());
580     dir->cd();
581
582     // Copy the input to the output 
583     TH1* outRaw    = static_cast<TH1*>(inRaw    ->Clone("measured"));
584     TH1* outTruth  = static_cast<TH1*>(inTruth  ->Clone("truth"));
585     TH1* outTruthA = static_cast<TH1*>(inTruthA ->Clone("truthAccepted"));
586     TH1* outTrgVtx = static_cast<TH1*>(inTrgVtx ->Clone("triggerVertex"));
587     TH2* outResp   = static_cast<TH2*>(inResp   ->Clone("response"));
588
589     // Make our response matrix 
590     RooUnfoldResponse matrix(0, 0, inResp);
591     
592     // Store regularization parameter 
593     Double_t             r        = regParam;
594     RooUnfold::Algorithm algo     = (RooUnfold::Algorithm)method;
595     RooUnfold*           unfolder = RooUnfold::New(algo, &matrix, inRaw, r);
596     unfolder->SetVerbose(1);
597
598     // Do the unfolding and get the result
599     TH1* res = unfolder->Hreco();
600     res->SetDirectory(0);
601
602     // Make a copy to store on the output 
603     TH1* outUnfold = static_cast<TH1*>(res->Clone("unfolded"));
604     TString tit(outUnfold->GetTitle());
605     tit.ReplaceAll("Unfold Reponse matrix", "Unfolded P(#it{N}_{ch})");
606     outUnfold->SetTitle(tit);
607
608     // Clone the unfolded results and divide by the trigger/vertex
609     // bias correction
610     TH1* outCorr   = static_cast<TH1*>(outUnfold->Clone("corrected"));
611     outCorr->Divide(inTrgVtx);
612     tit.ReplaceAll("Unfolded", "Corrected");
613     outCorr->SetTitle(tit);
614
615     // Now normalize the output to integral=1 
616     TH1*  hists[] = { outRaw, outUnfold, outCorr, 0 };
617     TH1** phist   = hists;
618     while (*phist) { 
619       TH1* h = *phist;
620       if (h) { 
621         Double_t intg = h->Integral(1, h->GetXaxis()->GetXmax());
622         h->Scale(1. / intg, "width");
623       }
624       phist++;
625     }
626     
627     // And make ratios
628     TH1* ratioTrue = static_cast<TH1*>(outCorr->Clone("ratioCorrTruth"));
629     tit = ratioTrue->GetTitle();
630     tit.ReplaceAll("Corrected", "Corrected/MC 'truth'");
631     ratioTrue->SetTitle(tit);
632     ratioTrue->Divide(outTruth);
633     ratioTrue->SetYTitle("P_{corrected}(#it{N}_{ch})/P_{truth}(#it{N}_{ch})");
634
635     TH1* ratioAcc  = static_cast<TH1*>(outUnfold->Clone("ratioUnfAcc"));
636     tit = ratioAcc->GetTitle();
637     tit.ReplaceAll("Unfolded", "Unfolded/MC selected");
638     ratioAcc->SetTitle(tit);
639     ratioAcc->Divide(outTruthA);
640     ratioAcc->SetYTitle("P_{unfolded}(#it{N}_{ch})/P_{MC}(#it{N}_{ch})");
641     
642
643     // Make a stack 
644     tit = measured->GetName();
645     tit.ReplaceAll("m", "-");
646     tit.ReplaceAll("p", "+");
647     tit.ReplaceAll("d", ".");
648     tit.ReplaceAll("_", "<#it{#eta}<");
649     THStack* stack = new THStack("all", tit);
650     stack->Add(outTruth,  "E2");
651     stack->Add(outTruthA, "E2");
652     stack->Add(outRaw,    "E1");
653     stack->Add(outUnfold, "E1");
654     stack->Add(outCorr,   "E1");
655     dir->Add(stack);
656
657     // Rest of the function is devoted to making the output look nice 
658     outRaw   ->SetDirectory(dir); 
659     outTruth ->SetDirectory(dir);  
660     outTruthA->SetDirectory(dir);  
661     outTrgVtx->SetDirectory(dir);  
662     outResp  ->SetDirectory(dir);  
663     outUnfold->SetDirectory(dir);   
664     outCorr  ->SetDirectory(dir); 
665
666     outRaw   ->SetMarkerStyle(20);  // Measured is closed
667     outTruth ->SetMarkerStyle(24);  // MC is open
668     outTruthA->SetMarkerStyle(24);  // MC is open
669     outTrgVtx->SetMarkerStyle(20);  // Derived is closed
670     outUnfold->SetMarkerStyle(20);  // Derived is closed   
671     outCorr  ->SetMarkerStyle(20);  // Derived is closed 
672
673     outRaw   ->SetMarkerSize(0.9); 
674     outTruth ->SetMarkerSize(1.6);  
675     outTruthA->SetMarkerSize(1.4);  
676     outTrgVtx->SetMarkerSize(1.0);  
677     outUnfold->SetMarkerSize(0.9);   
678     outCorr  ->SetMarkerSize(1.0);
679  
680     outRaw   ->SetMarkerColor(kColorMeasured); 
681     outTruth ->SetMarkerColor(kColorTruth);  
682     outTruthA->SetMarkerColor(kColorAccepted);  
683     outTrgVtx->SetMarkerColor(kColorTrgVtx);  
684     outUnfold->SetMarkerColor(kColorUnfolded);   
685     outCorr  ->SetMarkerColor(kColorCorrected); 
686
687     outRaw   ->SetFillColor(kColorError);     
688     outTruth ->SetFillColor(kColorError);  
689     outTruthA->SetFillColor(kColorError);  
690     outTrgVtx->SetFillColor(kColorError);  
691     outUnfold->SetFillColor(kColorError);   
692     outCorr  ->SetFillColor(kColorError); 
693
694     outRaw   ->SetFillStyle(0); 
695     outTruth ->SetFillStyle(1001);  
696     outTruthA->SetFillStyle(1001);  
697     outTrgVtx->SetFillStyle(0);  
698     outUnfold->SetFillStyle(0);   
699     outCorr  ->SetFillStyle(0);
700
701     outRaw   ->SetLineColor(kBlack); 
702     outTruth ->SetLineColor(kBlack);  
703     outTruthA->SetLineColor(kBlack);  
704     outTrgVtx->SetLineColor(kBlack);  
705     outUnfold->SetLineColor(kBlack);   
706     outCorr  ->SetLineColor(kBlack); 
707
708     // Legend 
709     TLegend* l = StackLegend(stack);
710     l->AddEntry(outRaw,     "Raw",                 "lp");
711     l->AddEntry(outTruth,   "MC 'truth'",          "fp");
712     l->AddEntry(outTruthA,  "MC 'truth' accepted", "fp");
713     l->AddEntry(outUnfold,  "Unfolded",            "lp");
714     l->AddEntry(outCorr,    "Corrected",           "lp");
715
716     return stack;
717   }
718   static void BinAttributes(Int_t i,
719                             Int_t&    open, 
720                             Int_t&    closed, 
721                             Float_t&  size,
722                             Double_t& factor) 
723   {
724     // --- Setup for markers -----------------------------------------
725     const Int_t   nMarkers = 7;
726     const Int_t   aClosed[] = { 20,  21,  22,  23,  29,  33,  34  };
727     const Int_t   aOpen[]   = { 24,  25,  26,  32,  30,  27,  28  };
728     const Float_t aSize[]   = { 1.1, 1.0, 1.2, 1.2, 1.2, 1.2, 1.0 };
729     Int_t         j         = i % nMarkers;
730
731     open   = aOpen[j];
732     closed = aClosed[j];
733     size   = aSize[j];
734     factor = TMath::Power(10, i);
735   }
736   /** 
737    * Add the bin histograms to our summary stacks 
738    * 
739    * @param bin       Bin stack
740    * @param i         Current off-set in the stacks 
741    * @param measured  All measured @f$ P(N_{ch})@f$ 
742    * @param truth     All MC truth @f$ P(N_{ch})@f$ 
743    * @param accepted  All MC accepted @f$ P(N_{ch})@f$ 
744    * @param unfolded  All unfolded @f$ P(N_{ch})@f$ 
745    * @param corrected All corrected @f$ P(N_{ch})@f$ 
746    */
747   void Bin2Stack(const THStack* bin, Int_t i, 
748                  THStack* measured, 
749                  THStack* truth, 
750                  THStack* accepted, 
751                  THStack* unfolded,
752                  THStack* corrected)
753   {
754     Int_t open, closed;
755     Double_t factor; 
756     Float_t  size;
757     BinAttributes(i, open, closed, size, factor);
758
759     TIter next(bin->GetHists());
760     TH1*  h = 0;
761     while ((h = static_cast<TH1*>(next()))) {
762       THStack* tmp = 0;
763       Int_t    col = h->GetMarkerColor();
764       Int_t    sty = 0;
765       switch (col) { 
766       case kColorMeasured:  tmp = measured;   sty = closed;  break;
767       case kColorTruth:     tmp = truth;      sty = open;    break;
768       case kColorAccepted:  tmp = accepted;   sty = open;    break;
769       case kColorUnfolded:  tmp = unfolded;   sty = closed;  break;
770       case kColorCorrected: tmp = corrected;  sty = closed;  break;
771       default: continue; 
772       }
773       // Now clone, and add to the appropriate stack 
774       TH1* cln = static_cast<TH1*>(h->Clone(h->GetName()));
775       cln->SetDirectory(0);
776       cln->SetMarkerStyle(sty);
777       cln->SetMarkerSize(size);
778       cln->Scale(factor); // Scale by 10^i
779
780       // Make sure we do not get the old legend 
781       TObject* tst = cln->FindObject("legend");
782       if (tst) cln->GetListOfFunctions()->Remove(tst);
783
784       tmp->Add(cln, next.GetOption());
785     }
786     
787     // Add entries to our stacks 
788     TString   txt      = bin->GetTitle();
789     if      (i == 0) txt.Append(" (#times1)");
790     else if (i == 1) txt.Append(" (#times10)");
791     else             txt.Append(Form(" (#times10^{%d})", i));
792     THStack*  stacks[] = { measured, truth, accepted, unfolded, corrected, 0 };
793     THStack** pstack   = stacks;
794     while (*pstack) { 
795       TLegend* leg = StackLegend(*pstack);
796       pstack++;
797       if (!leg) continue;
798       
799       TObject* dummy = 0;
800       TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
801       e->SetMarkerStyle(closed);
802       e->SetMarkerSize(1.2*size);
803       e->SetMarkerColor(kBlack);
804       e->SetFillColor(0);
805       e->SetFillStyle(0);
806       e->SetLineColor(kBlack);
807     }
808   }
809   /** 
810    * Add distributions from other experiments to stacks 
811    * 
812    * @param name     Name of current bin 
813    * @param i        Index of current bin
814    * @param sNN      Center of mass energy 
815    * @param allALICE Stack of ALICE data 
816    * @param allCMS   Stack of CMS data 
817    */
818   void Other2Stack(const TString& name, Int_t i,
819                    UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS) 
820   {
821     if (!allALICE && !allCMS) return;
822
823     TString tmp(name);
824     tmp.ReplaceAll("p", "+");
825     tmp.ReplaceAll("m", "-");
826     tmp.ReplaceAll("_", " ");
827     TObjArray* tokens = tmp.Tokenize(" ");
828     if (!tokens || tokens->GetEntriesFast() < 2) { 
829       Error("Other2Stack", "Failed to decode eta range from %s", name.Data());
830       if (tokens) tokens->Delete();
831       return;
832     }
833     Double_t eta1 = static_cast<TObjString*>(tokens->At(0))->String().Atof();
834     Double_t eta2 = static_cast<TObjString*>(tokens->At(1))->String().Atof();
835     tokens->Delete();
836     
837     if (TMath::Abs(eta2-eta1) > 1e3) 
838       // Not symmetric bin 
839       return;
840     Double_t aEta = TMath::Abs(eta1);
841
842     Int_t open, closed;
843     Double_t factor; 
844     Float_t  size;
845     BinAttributes(i, open, closed, size, factor);
846
847     if (allALICE) {
848       TGraphAsymmErrors* g = GetOther(1, aEta, sNN, factor);
849       if (g) {
850         g->SetMarkerStyle(closed);
851         g->SetMarkerColor(kColorALICE);
852         g->SetMarkerSize(size);
853         allALICE->Add(g, "p same");
854       }
855     }
856     if (allCMS) {
857       TGraphAsymmErrors* g = GetOther(0, aEta, sNN, factor);
858       if (g) {
859         g->SetMarkerStyle(closed);
860         g->SetMarkerColor(kColorCMS);
861         g->SetMarkerSize(size);
862         allCMS->Add(g, "p same");
863       }
864     }
865
866     
867   }
868   /** 
869    * Get or create a stack legend.  This is done by adding a TLegend
870    * object to the list of functions for the first histogram in the
871    * stack.
872    * 
873    * @param stack Stack to get the legend from/modify 
874    * 
875    * @return The legend object or null
876    */
877   TLegend* StackLegend(THStack* stack) 
878   {
879     TList* hists = stack->GetHists();
880     if (!hists) return 0;
881     
882     TObject* first = hists->First();
883     if (!first) return 0;
884
885     TH1*    hist = static_cast<TH1*>(first);
886     TList*  list = hist->GetListOfFunctions();
887     TObject* o   = list->FindObject("legend");
888     if (o) return static_cast<TLegend*>(o);
889     
890     TLegend* l = new TLegend(0.65, 0.65, 0.9, 0.9, "", "NDC");
891     l->SetName("legend");
892     l->SetBorderSize(0);
893     l->SetFillColor(0);
894     l->SetFillStyle(0);
895     l->SetTextFont(42);
896     list->Add(l);
897     
898     return l;
899   }
900   
901   /* =================================================================
902    *
903    * Measurements from other sources, such as published ALICE, CMS, ...
904    */
905   TGraphAsymmErrors* GetOther(UShort_t type, Double_t eta, UShort_t sNN,
906                               Int_t factor)
907   {
908     TString oC = Form("OtherPNch::GetData(%hu,%f,%hu)", 
909                       type, eta, sNN);
910     TGraphAsymmErrors* g = 
911       reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
912     if (!g) { 
913       Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
914               type, eta, sNN);
915       return 0;
916     }
917
918
919     for (Int_t j = 0; j < g->GetN(); j++) { 
920       g->SetPoint(j, g->GetX()[j], g->GetY()[j]*factor);
921       g->SetPointError(j, g->GetEXlow()[j], g->GetEXhigh()[j], 
922                        g->GetEYlow()[j]*factor, g->GetEYhigh()[j]*factor);
923     }
924     return g;
925   }    
926 };
927
928 void
929 UnfoldMultDists(const TString& method="Bayes", 
930                 Double_t       regParam=1e-30,
931                 const TString& measured="forward_multdists.root", 
932                 const TString& corrections="")
933 {
934   Unfolder m;
935   m.Run(measured, corrections, method, regParam);
936 }