A better way to specify the Nch axis for the MultDists task.
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / UnfoldMultDists.C
CommitLineData
bfab35d9 1#include <TFile.h>
2#include <TList.h>
3#include <TH1.h>
4#include <TH2.h>
281a2bf8 5#include <THStack.h>
997ba0f4 6#include <TLegend.h>
7#include <TLegendEntry.h>
281a2bf8 8#include <TClass.h>
9#include <TRegexp.h>
10#include <TMath.h>
997ba0f4 11#include <TParameter.h>
997ba0f4 12#include <TMultiGraph.h>
281a2bf8 13#include <TGraphAsymmErrors.h>
997ba0f4 14#include "RooUnfold.h"
15#include "RooUnfoldResponse.h"
281a2bf8 16#include <fstream>
bfab35d9 17
281a2bf8 18/**
19 * Class to do unfolding of raw histograms produced by AliForwardMultDists
bfab35d9 20 *
bfab35d9 21 */
281a2bf8 22struct Unfolder
bfab35d9 23{
281a2bf8 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
bfab35d9 38 };
281a2bf8 39
bfab35d9 40 /**
281a2bf8 41 * Constructor
bfab35d9 42 */
43 Unfolder() {}
44 /**
281a2bf8 45 * Get a top collection from a file
bfab35d9 46 *
281a2bf8 47 * @param fileName Name of file
48 * @param results Wheter it's the results collection or not
49 *
50 * @return Collection or null
bfab35d9 51 */
281a2bf8 52 static TCollection* GetTop(const TString& fileName, Bool_t results=false)
bfab35d9 53 {
281a2bf8 54 TFile* file = TFile::Open(fileName, "READ");
55 if (!file) {
56 Error("GetTop", "Failed to open %s", fileName.Data());
57 return 0;
bfab35d9 58 }
281a2bf8 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;
bfab35d9 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 */
281a2bf8 77 static TObject* GetObject(TCollection* c, const TString& name,
78 TClass* cl, Bool_t verbose=true)
bfab35d9 79 {
80 TObject* o = c->FindObject(name);
81 if (!o) {
281a2bf8 82 if (verbose)
83 Warning("GetObject", "%s not found in %s", name.Data(), c->GetName());
bfab35d9 84 return 0;
85 }
86 if (cl && !o->IsA()->InheritsFrom(cl)) {
281a2bf8 87 if (verbose)
88 Warning("GetCollection", "%s is not a %s but a %s",
89 name.Data(), cl->GetName(), o->ClassName());
bfab35d9 90 return 0;
91 }
92 return o;
93 }
bfab35d9 94 /**
95 * Get a collection
96 *
97 * @param c
98 * @param name
99 *
100 * @return
101 */
281a2bf8 102 static TCollection* GetCollection(TCollection* c,
103 const TString& name,
104 Bool_t verbose=-true)
bfab35d9 105 {
281a2bf8 106 return static_cast<TCollection*>(GetObject(c, name,
107 TCollection::Class(),
108 verbose));
bfab35d9 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 */
281a2bf8 118 static TH1* GetH1(TCollection* c, const TString& name, Bool_t verbose=true)
bfab35d9 119 {
281a2bf8 120 return static_cast<TH1*>(GetObject(c, name, TH1::Class(), verbose));
bfab35d9 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 */
281a2bf8 130 static TH2* GetH2(TCollection* c, const TString& name, Bool_t verbose=true)
bfab35d9 131 {
281a2bf8 132 return static_cast<TH2*>(GetObject(c, name, TH2::Class(), verbose));
bfab35d9 133 }
281a2bf8 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)
bfab35d9 143 {
281a2bf8 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;
bfab35d9 199 }
281a2bf8 200 pMethod++;
bfab35d9 201 }
281a2bf8 202 if (pMethod->id == 0xDeadBeef)
203 Error("MethodId", "Unknown unfolding method: %s", method.Data());
204
205 return pMethod->id;
bfab35d9 206 }
281a2bf8 207
bfab35d9 208 /**
281a2bf8 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
bfab35d9 223 *
281a2bf8 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)
bfab35d9 269 {
281a2bf8 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;
281a2bf8 283 ULong_t trig;
284 Double_t minZ;
285 Double_t maxZ;
286 GetParameter(mTop, "sys", sys);
287 GetParameter(mTop, "snn", sNN);
281a2bf8 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);
bfab35d9 294
281a2bf8 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 }
bfab35d9 301
281a2bf8 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
52b36573 308 SaveInformation(out,meth,mId,regParam,sys,sNN,trig,minZ,maxZ,
309 corrFile.IsNull());
281a2bf8 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 }
52b36573 338 /**
339 * Append an & to a string and the next term.
340 *
341 * @param trg Output string
342 * @param what Term
343 */
281a2bf8 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
281a2bf8 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,
52b36573 371 Bool_t self) const
281a2bf8 372 {
373 dir->cd();
bfab35d9 374
52b36573 375 TParameter<bool>* pM = new TParameter<bool>("self", self);
376 pM->SetBit(BIT(19));
377 pM->Write();
378
281a2bf8 379 TNamed* outMeth = new TNamed("method", method.Data());
380 outMeth->SetUniqueID(mId);
381 outMeth->Write();
bfab35d9 382
281a2bf8 383 TParameter<double>* pR = new TParameter<double>("regParam", regParam);
52b36573 384 pR->SetBit(BIT(19));
281a2bf8 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);
52b36573 447 pY->SetBit(BIT(19));
281a2bf8 448 pY->Write();
bfab35d9 449
281a2bf8 450 TParameter<double>* pZ = new TParameter<double>("maxIpZ", maxIpZ);
52b36573 451 pZ->SetBit(BIT(19));
281a2bf8 452 pZ->Write();
281a2bf8 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})");
52b36573 504 THStack* allRatio = (sys != 1 ? 0 :
505 new THStack("ratios", "Ratios to other"));
281a2bf8 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.
bfab35d9 512 static TRegexp regex("[pm][0-9]d[0-9]*_[pm][0-9]d[0-9]*");
281a2bf8 513 TIter next(measured);
bfab35d9 514 TObject* o = 0;
281a2bf8 515 Int_t i = 0;
997ba0f4 516 Double_t r = regParam;
bfab35d9 517 while ((o = next())) {
281a2bf8 518 // Go back to where we where
519 dir->cd();
520
bfab35d9 521 // if not a collection, don't bother
522 if (!o->IsA()->InheritsFrom(TCollection::Class())) continue;
997ba0f4 523
bfab35d9 524 // If it doesn't match our regular expression, don't bother
525 TString n(o->GetName());
526 if (n.Index(regex) == kNPOS) {
997ba0f4 527 // Warning("ScanType", "%s in %s doesn't match eta range regexp",
528 // n.Data(), real->GetName());
bfab35d9 529 continue;
530 }
281a2bf8 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
52b36573 538 TH1* result = 0;
281a2bf8 539 Bin2Stack(binS, i, allMeasured, allTruth, allTruthA,
52b36573 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);
281a2bf8 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;
bfab35d9 564 }
52b36573 565 if (allRatio && allRatio->GetHists()) {
566 if (allRatio->GetHists()->GetEntries() > 0)
567 dir->Add(allRatio);
568 else
569 delete allRatio;
570 }
bfab35d9 571 }
572 /**
281a2bf8 573 * Process a single eta bin
bfab35d9 574 *
281a2bf8 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)
bfab35d9 588 {
281a2bf8 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)
bfab35d9 597 return 0;
bfab35d9 598
281a2bf8 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);
bfab35d9 612
281a2bf8 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");
bfab35d9 644 }
281a2bf8 645 phist++;
bfab35d9 646 }
bfab35d9 647
281a2bf8 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})");
bfab35d9 662
bfab35d9 663
281a2bf8 664 // Make a stack
665 tit = measured->GetName();
bfab35d9 666 tit.ReplaceAll("m", "-");
667 tit.ReplaceAll("p", "+");
668 tit.ReplaceAll("d", ".");
669 tit.ReplaceAll("_", "<#it{#eta}<");
bfab35d9 670 THStack* stack = new THStack("all", tit);
281a2bf8 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");
bfab35d9 676 dir->Add(stack);
677
281a2bf8 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
bfab35d9 737 return stack;
738 }
281a2bf8 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 }
bfab35d9 757 /**
281a2bf8 758 * Add the bin histograms to our summary stacks
bfab35d9 759 *
281a2bf8 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$
bfab35d9 767 */
281a2bf8 768 void Bin2Stack(const THStack* bin, Int_t i,
769 THStack* measured,
770 THStack* truth,
771 THStack* accepted,
772 THStack* unfolded,
52b36573 773 THStack* corrected,
774 TH1*& result)
bfab35d9 775 {
281a2bf8 776 Int_t open, closed;
777 Double_t factor;
778 Float_t size;
779 BinAttributes(i, open, closed, size, factor);
bfab35d9 780
281a2bf8 781 TIter next(bin->GetHists());
bfab35d9 782 TH1* h = 0;
281a2bf8 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;
bfab35d9 794 }
281a2bf8 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
52b36573 801 if (col == kColorCorrected) result = cln;
bfab35d9 802
281a2bf8 803 // Make sure we do not get the old legend
804 TObject* tst = cln->FindObject("legend");
805 if (tst) cln->GetListOfFunctions()->Remove(tst);
bfab35d9 806
281a2bf8 807 tmp->Add(cln, next.GetOption());
bfab35d9 808 }
281a2bf8 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;
997ba0f4 821
281a2bf8 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);
997ba0f4 830 }
bfab35d9 831 }
281a2bf8 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,
52b36573 842 UShort_t sNN, TMultiGraph* allALICE, TMultiGraph* allCMS,
843 TGraph*& alice, TGraph*& cms)
bfab35d9 844 {
281a2bf8 845 if (!allALICE && !allCMS) return;
846
847 TString tmp(name);
848 tmp.ReplaceAll("p", "+");
849 tmp.ReplaceAll("m", "-");
850 tmp.ReplaceAll("_", " ");
52b36573 851 tmp.ReplaceAll("d", ".");
281a2bf8 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;
bfab35d9 857 }
281a2bf8 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();
bfab35d9 861
281a2bf8 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");
52b36573 879 alice = g;
281a2bf8 880 }
bfab35d9 881 }
281a2bf8 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");
52b36573 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);
281a2bf8 934 }
52b36573 935 all->Add(r);
936 pg++;
bfab35d9 937 }
52b36573 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));
281a2bf8 946
52b36573 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);
bfab35d9 955 }
52b36573 956
bfab35d9 957 /**
281a2bf8 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.
bfab35d9 961 *
281a2bf8 962 * @param stack Stack to get the legend from/modify
963 *
964 * @return The legend object or null
bfab35d9 965 */
281a2bf8 966 TLegend* StackLegend(THStack* stack)
bfab35d9 967 {
281a2bf8 968 TList* hists = stack->GetHists();
969 if (!hists) return 0;
970
971 TObject* first = hists->First();
972 if (!first) return 0;
bfab35d9 973
281a2bf8 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;
bfab35d9 988 }
281a2bf8 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)
bfab35d9 996 {
52b36573 997 TString oC = Form("OtherPNch::GetData(%hu,%f,%hu);",
281a2bf8 998 type, eta, sNN);
999 TGraphAsymmErrors* g =
1000 reinterpret_cast<TGraphAsymmErrors*>(gROOT->ProcessLine(oC));
1001 if (!g) {
52b36573 1002 // Warning("GetOther", "No other data found for type=%d eta=%f sNN=%d",
1003 // type, eta, sNN);
281a2bf8 1004 return 0;
bfab35d9 1005 }
52b36573 1006
281a2bf8 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);
bfab35d9 1012 }
1013 return g;
281a2bf8 1014 }
bfab35d9 1015};
281a2bf8 1016
1017void
1018UnfoldMultDists(const TString& method="Bayes",
1019 Double_t regParam=1e-30,
1020 const TString& measured="forward_multdists.root",
1021 const TString& corrections="")
bfab35d9 1022{
281a2bf8 1023 Unfolder m;
1024 m.Run(measured, corrections, method, regParam);
bfab35d9 1025}