Fix some documentation issues
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / UnfoldMultDists.C
CommitLineData
671df6c9 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 */
bfab35d9 11#include <TFile.h>
12#include <TList.h>
13#include <TH1.h>
14#include <TH2.h>
281a2bf8 15#include <THStack.h>
997ba0f4 16#include <TLegend.h>
17#include <TLegendEntry.h>
281a2bf8 18#include <TClass.h>
19#include <TRegexp.h>
20#include <TMath.h>
997ba0f4 21#include <TParameter.h>
997ba0f4 22#include <TMultiGraph.h>
281a2bf8 23#include <TGraphAsymmErrors.h>
997ba0f4 24#include "RooUnfold.h"
25#include "RooUnfoldResponse.h"
281a2bf8 26#include <fstream>
bfab35d9 27
281a2bf8 28/**
29 * Class to do unfolding of raw histograms produced by AliForwardMultDists
bfab35d9 30 *
671df6c9 31 * @ingroup pwglf_forward_multdist
bfab35d9 32 */
281a2bf8 33struct Unfolder
bfab35d9 34{
281a2bf8 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
bfab35d9 49 };
281a2bf8 50
bfab35d9 51 /**
281a2bf8 52 * Constructor
bfab35d9 53 */
54 Unfolder() {}
55 /**
281a2bf8 56 * Get a top collection from a file
bfab35d9 57 *
281a2bf8 58 * @param fileName Name of file
59 * @param results Wheter it's the results collection or not
60 *
61 * @return Collection or null
bfab35d9 62 */
281a2bf8 63 static TCollection* GetTop(const TString& fileName, Bool_t results=false)
bfab35d9 64 {
281a2bf8 65 TFile* file = TFile::Open(fileName, "READ");
66 if (!file) {
67 Error("GetTop", "Failed to open %s", fileName.Data());
68 return 0;
bfab35d9 69 }
281a2bf8 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;
bfab35d9 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
c8b1a7db 85 * @param verbose Be verbose
bfab35d9 86 *
87 * @return Pointer to object or null
88 */
281a2bf8 89 static TObject* GetObject(TCollection* c, const TString& name,
90 TClass* cl, Bool_t verbose=true)
bfab35d9 91 {
92 TObject* o = c->FindObject(name);
93 if (!o) {
281a2bf8 94 if (verbose)
95 Warning("GetObject", "%s not found in %s", name.Data(), c->GetName());
bfab35d9 96 return 0;
97 }
98 if (cl && !o->IsA()->InheritsFrom(cl)) {
281a2bf8 99 if (verbose)
100 Warning("GetCollection", "%s is not a %s but a %s",
101 name.Data(), cl->GetName(), o->ClassName());
bfab35d9 102 return 0;
103 }
104 return o;
105 }
bfab35d9 106 /**
107 * Get a collection
108 *
c8b1a7db 109 * @param c Parent collection
110 * @param name Name of object to findf
111 * @param verbose Be verbose
bfab35d9 112 *
113 * @return
114 */
281a2bf8 115 static TCollection* GetCollection(TCollection* c,
116 const TString& name,
117 Bool_t verbose=-true)
bfab35d9 118 {
281a2bf8 119 return static_cast<TCollection*>(GetObject(c, name,
120 TCollection::Class(),
121 verbose));
bfab35d9 122 }
123 /**
124 * Get a 1D histogram from a collection
125 *
126 * @param c Collection
127 * @param name Nanme of histogram
c8b1a7db 128 * @param verbose Be verbose
bfab35d9 129 *
130 * @return Pointer to object or null
131 */
281a2bf8 132 static TH1* GetH1(TCollection* c, const TString& name, Bool_t verbose=true)
bfab35d9 133 {
281a2bf8 134 return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verbose));
bfab35d9 135 }
136 /**
137 * Get a 2D histogram from a collection
138 *
139 * @param c Collection
140 * @param name Nanme of histogram
c8b1a7db 141 * @param verbose Be verbose
bfab35d9 142 *
143 * @return Pointer to object or null
144 */
281a2bf8 145 static TH2* GetH2(TCollection* c, const TString& name, Bool_t verbose=true)
bfab35d9 146 {
281a2bf8 147 return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verbose));
bfab35d9 148 }
281a2bf8 149 /**
150 * Get an unsigned short parameter from the collection
151 *
152 * @param c Collection
153 * @param name Parameter name
c8b1a7db 154 * @param v Value
281a2bf8 155 *
156 * @return Value
157 */
158 static void GetParameter(TCollection* c, const TString& name, UShort_t& v)
bfab35d9 159 {
281a2bf8 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
c8b1a7db 168 * @param v Value
281a2bf8 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
c8b1a7db 182 * @param v Value
281a2bf8 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 *
c8b1a7db 194 * @param method Method
281a2bf8 195 *
c8b1a7db 196 * @return Method identifier
281a2bf8 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;
bfab35d9 217 }
281a2bf8 218 pMethod++;
bfab35d9 219 }
281a2bf8 220 if (pMethod->id == 0xDeadBeef)
221 Error("MethodId", "Unknown unfolding method: %s", method.Data());
222
223 return pMethod->id;
bfab35d9 224 }
281a2bf8 225
bfab35d9 226 /**
281a2bf8 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
bfab35d9 241 *
281a2bf8 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 *
c8b1a7db 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>
281a2bf8 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
c8b1a7db 277 * has the expected @a corrFile @e in @e addition to the
281a2bf8 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)
bfab35d9 287 {
281a2bf8 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;
281a2bf8 301 ULong_t trig;
302 Double_t minZ;
303 Double_t maxZ;
304 GetParameter(mTop, "sys", sys);
c8b1a7db 305 GetParameter(mTop, "sNN", sNN);
281a2bf8 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);
bfab35d9 312
281a2bf8 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 }
bfab35d9 319
281a2bf8 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
52b36573 326 SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
327 corrFile.IsNull());
281a2bf8 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 }
52b36573 356 /**
357 * Append an & to a string and the next term.
358 *
359 * @param trg Output string
360 * @param what Term
361 */
281a2bf8 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
c8b1a7db 379 * @param self Self-consistency check
281a2bf8 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,
52b36573 390 Bool_t self) const
281a2bf8 391 {
392 dir->cd();
bfab35d9 393
52b36573 394 TParameter<bool>* pM = new TParameter<bool>("self", self);
395 pM->SetBit(BIT(19));
396 pM->Write();
397
281a2bf8 398 TNamed* outMeth = new TNamed("method", method.Data());
399 outMeth->SetUniqueID(mId);
400 outMeth->Write();
bfab35d9 401
281a2bf8 402 TParameter<double>* pR = new TParameter<double>("regParam", regParam);
52b36573 403 pR->SetBit(BIT(19));
281a2bf8 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);
52b36573 466 pY->SetBit(BIT(19));
281a2bf8 467 pY->Write();
bfab35d9 468
281a2bf8 469 TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
52b36573 470 pZ->SetBit(BIT(19));
281a2bf8 471 pZ->Write();
281a2bf8 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.
c8b1a7db 500 * @param sys Collision system
501 * @param sNN Collision energy
281a2bf8 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})");
52b36573 525 THStack* allRatio = (sys != 1 ? 0 :
526 new THStack("ratios", "Ratios to other"));
281a2bf8 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.
bfab35d9 533 static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
281a2bf8 534 TIter next(measured);
bfab35d9 535 TObject* o = 0;
281a2bf8 536 Int_t i = 0;
997ba0f4 537 Double_t r = regParam;
bfab35d9 538 while ((o = next())) {
281a2bf8 539 // Go back to where we where
540 dir->cd();
541
bfab35d9 542 // if not a collection, don't bother
543 if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
997ba0f4 544
bfab35d9 545 // If it doesn't match our regular expression, don't bother
546 TString n(o->GetName());
547 if (n.Index(regex) == kNPOS) {
997ba0f4 548 // Warning("ScanType", "%s in %s doesn't match eta range regexp",
549 // n.Data(), real->GetName());
bfab35d9 550 continue;
551 }
281a2bf8 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
52b36573 559 TH1* result = 0;
281a2bf8 560 Bin2Stack(binS, i, allMeasured, allTruth, allTruthA,
52b36573 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);
281a2bf8 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;
bfab35d9 585 }
52b36573 586 if (allRatio && allRatio->GetHists()) {
587 if (allRatio->GetHists()->GetEntries() > 0)
588 dir->Add(allRatio);
589 else
590 delete allRatio;
591 }
bfab35d9 592 }
593 /**
281a2bf8 594 * Process a single eta bin
bfab35d9 595 *
281a2bf8 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)
bfab35d9 609 {
281a2bf8 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)
bfab35d9 618 return 0;
bfab35d9 619
281a2bf8 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);
bfab35d9 633
281a2bf8 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");
bfab35d9 665 }
281a2bf8 666 phist++;
bfab35d9 667 }
bfab35d9 668
281a2bf8 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})");
bfab35d9 683
bfab35d9 684
281a2bf8 685 // Make a stack
686 tit = measured->GetName();
bfab35d9 687 tit.ReplaceAll("m", "-");
688 tit.ReplaceAll("p", "+");
689 tit.ReplaceAll("d", ".");
690 tit.ReplaceAll("_", "<#it{#eta}<");
bfab35d9 691 THStack* stack = new THStack("all", tit);
281a2bf8 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");
bfab35d9 697 dir->Add(stack);
698
281a2bf8 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
bfab35d9 758 return stack;
759 }
281a2bf8 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 }
bfab35d9 778 /**
281a2bf8 779 * Add the bin histograms to our summary stacks
bfab35d9 780 *
281a2bf8 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$
c8b1a7db 788 * @param result The result in this bin
bfab35d9 789 */
281a2bf8 790 void Bin2Stack(const THStack* bin, Int_t i,
791 THStack* measured,
792 THStack* truth,
793 THStack* accepted,
794 THStack* unfolded,
52b36573 795 THStack* corrected,
796 TH1*& result)
bfab35d9 797 {
281a2bf8 798 Int_t open, closed;
799 Double_t factor;
800 Float_t size;
801 BinAttributes(i, open, closed, size, factor);
bfab35d9 802
281a2bf8 803 TIter next(bin->GetHists());
bfab35d9 804 TH1* h = 0;
281a2bf8 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;
bfab35d9 816 }
281a2bf8 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
52b36573 823 if (col == kColorCorrected) result = cln;
bfab35d9 824
281a2bf8 825 // Make sure we do not get the old legend
826 TObject* tst = cln->FindObject("legend");
827 if (tst) cln->GetListOfFunctions()->Remove(tst);
bfab35d9 828
281a2bf8 829 tmp->Add(cln, next.GetOption());
bfab35d9 830 }
281a2bf8 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;
997ba0f4 843
281a2bf8 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);
997ba0f4 852 }
bfab35d9 853 }
281a2bf8 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
c8b1a7db 862 * @param alice Possible ALICE result on return
863 * @param cms Possible CMS result on return
281a2bf8 864 */
865 void Other2Stack(const TString& name, Int_t i,
52b36573 866 UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
867 TGraph*& alice, TGraph*& cms)
bfab35d9 868 {
281a2bf8 869 if (!allALICE && !allCMS) return;
870
871 TString tmp(name);
872 tmp.ReplaceAll("p", "+");
873 tmp.ReplaceAll("m", "-");
874 tmp.ReplaceAll("_", " ");
52b36573 875 tmp.ReplaceAll("d", ".");
281a2bf8 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;
bfab35d9 881 }
281a2bf8 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();
bfab35d9 885
281a2bf8 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");
52b36573 903 alice = g;
281a2bf8 904 }
bfab35d9 905 }
281a2bf8 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");
52b36573 913 cms = g;
914 }
915 }
916 }
917 /**
918 * Create ratios to other data
919 *
c8b1a7db 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
52b36573 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);
281a2bf8 958 }
52b36573 959 all->Add(r);
960 pg++;
bfab35d9 961 }
52b36573 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));
281a2bf8 970
52b36573 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);
bfab35d9 979 }
52b36573 980
bfab35d9 981 /**
281a2bf8 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.
bfab35d9 985 *
281a2bf8 986 * @param stack Stack to get the legend from/modify
987 *
988 * @return The legend object or null
bfab35d9 989 */
281a2bf8 990 TLegend* StackLegend(THStack* stack)
bfab35d9 991 {
281a2bf8 992 TList* hists = stack->GetHists();
993 if (!hists) return 0;
994
995 TObject* first = hists->First();
996 if (!first) return 0;
bfab35d9 997
281a2bf8 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;
bfab35d9 1012 }
281a2bf8 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)
bfab35d9 1020 {
52b36573 1021 TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);",
281a2bf8 1022 type, eta, sNN);
1023 TGraphAsymmErrors* g =
1024 reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
1025 if (!g) {
52b36573 1026 // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
1027 // type, eta, sNN);
281a2bf8 1028 return 0;
bfab35d9 1029 }
52b36573 1030
281a2bf8 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);
bfab35d9 1036 }
1037 return g;
281a2bf8 1038 }
bfab35d9 1039};
281a2bf8 1040
1041void
1042UnfoldMultDists(const TString& method="Bayes",
1043 Double_t regParam=1e-30,
1044 const TString& measured="forward_multdists.root",
1045 const TString& corrections="")
bfab35d9 1046{
281a2bf8 1047 Unfolder m;
1048 m.Run(measured, corrections, method, regParam);
bfab35d9 1049}