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