4f09623b |
1 | /** |
671df6c9 |
2 | * @file UnfoldMult.C |
4f09623b |
3 | * @author Valentina Zaccolo |
4 | * @date Tue Mar 05 12:56:56 2013 |
5 | * |
6 | * @brief Bayesian Method Unfolding Script |
7 | * |
671df6c9 |
8 | * @ingroup pwglf_forward_scripts |
9 | * @ingroup pwglf_forward_multdist |
4f09623b |
10 | */ |
11 | #include <iostream> |
12 | #include <TList.h> |
13 | #include <TFile.h> |
14 | #include <TObject.h> |
15 | #include <TDirectory.h> |
16 | #include "TH1D.h" |
17 | #include "TH2D.h" |
18 | #include <TROOT.h> |
19 | #include <TSystem.h> |
20 | #include "RooUnfoldBayes.h" |
21 | #include "RooUnfoldResponse.h" |
22 | #include <TMath.h> |
671df6c9 |
23 | /** |
24 | * @defgroup pwglf_forward_scripts_unfoldmult UnfoldMult stuff |
25 | * |
26 | * @ingroup pwglf_forward_multdist |
27 | */ |
4f09623b |
28 | /** |
29 | * Get an object from a directory |
30 | * |
31 | * @param d Directory |
32 | * @param name Name of object |
33 | * |
34 | * @return Pointer to object or null if not found |
671df6c9 |
35 | * |
36 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
37 | */ |
38 | TObject* GetObject(TDirectory* d, const char* name) |
39 | { |
40 | TObject* o = d->Get(name); |
41 | if (!o) { |
42 | Error("GetCollection", "%s not found in %s", name, d->GetName()); |
43 | return 0; |
44 | } |
45 | return o; |
46 | } |
47 | /** |
48 | * Get an object from a collection |
49 | * |
50 | * @param d Collection |
51 | * @param name Name of object |
52 | * |
53 | * @return Pointer to object or null if not found |
671df6c9 |
54 | * |
55 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
56 | */ |
57 | TObject* GetObject(TCollection* d, const char* name) |
58 | { |
59 | TObject* o = d->FindObject(name); |
60 | if (!o) { |
61 | Error("GetCollection", "%s not found in %s", name, d->GetName()); |
62 | d->ls(); |
63 | return 0; |
64 | } |
65 | return o; |
66 | } |
67 | /** |
68 | * Get a collection from a directory |
69 | * |
70 | * @param d Parent directory |
71 | * @param name Name of collection |
72 | * |
73 | * @return Pointer to collection or null |
671df6c9 |
74 | * |
75 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
76 | */ |
77 | TCollection* GetCollection(TDirectory* d, const char* name) |
78 | { |
79 | return static_cast<TCollection*>(GetObject(d, name)); |
80 | } |
81 | /** |
82 | * Get a collection from a collection |
83 | * |
84 | * @param d Parent collection |
85 | * @param name Name of collection |
86 | * |
87 | * @return Pointer to collection or null |
671df6c9 |
88 | * |
89 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
90 | */ |
91 | TCollection* GetCollection(TCollection* d, const char* name) |
92 | { |
93 | return static_cast<TCollection*>(GetObject(d, name)); |
94 | } |
95 | /** |
96 | * Get a 1D-histogram from a directory |
97 | * |
98 | * @param d Parent directory |
99 | * @param name Name of histogram |
100 | * |
101 | * @return Pointer to histogram or null |
671df6c9 |
102 | * |
103 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
104 | */ |
105 | TH1* GetH1(TDirectory* d, const char* name) |
106 | { |
107 | return static_cast<TH1*>(GetObject(d, name)); |
108 | } |
109 | /** |
110 | * Get a 1D-histogram from a collection |
111 | * |
112 | * @param d Parent collectoin |
113 | * @param name Name of histogram |
114 | * |
115 | * @return Pointer to histogram or null |
671df6c9 |
116 | * |
117 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
118 | */ |
119 | TH1* GetH1(TCollection* d, const char* name) |
120 | { |
121 | return static_cast<TH1*>(GetObject(d, name)); |
122 | } |
123 | /** |
124 | * Get a 2D-histogram from a directory |
125 | * |
126 | * @param d Parent directory |
127 | * @param name Name of histogram |
128 | * |
129 | * @return Pointer to histogram or null |
130 | */ |
131 | TH2* GetH2(TDirectory* d, const char* name) |
132 | { |
133 | return static_cast<TH2*>(GetObject(d, name)); |
134 | } |
135 | /** |
136 | * Get a 2D-histogram from a collection |
137 | * |
138 | * @param d Parent collection |
139 | * @param name Name of histogram |
140 | * |
141 | * @return Pointer to histogram or null |
671df6c9 |
142 | * |
143 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
144 | */ |
145 | TH2* GetH2(TCollection* d, const char* name) |
146 | { |
147 | return static_cast<TH2*>(GetObject(d, name)); |
148 | } |
149 | /** |
150 | * Open a file |
151 | * |
152 | * @param filename File name |
153 | * @param readNotWrite If true, open read-only |
154 | * |
155 | * @return Pointer to file, or null |
671df6c9 |
156 | * |
157 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
158 | */ |
159 | TFile* OpenFile(const char* filename, Bool_t readNotWrite) |
160 | { |
161 | TFile* ret = TFile::Open(filename, readNotWrite ? "READ" : "RECREATE"); |
162 | if (!ret) { |
163 | Error("OpenFile", "Failed to open the file \"%s\"", filename); |
164 | return 0; |
165 | } |
166 | return ret; |
167 | } |
671df6c9 |
168 | /** |
169 | * A type definition to make it easier to switch to TDirectory structure |
170 | * |
171 | * @ingroup pwglf_forward_scripts_unfoldmult |
172 | */ |
4f09623b |
173 | typedef TCollection Dir; |
174 | |
175 | /** |
176 | * Get histograms for @f$\eta@f$-bin from @a lowLim to @a highLim and |
177 | * process it . |
178 | * |
179 | * @param lowLim Lower @f$\eta@f$ limit |
180 | * @param highLim Upper @f$\eta@f$ limit |
181 | * @param resp Directory containing response matrices |
182 | * @param data Directory containing data histogram |
183 | * @param out Output directory |
671df6c9 |
184 | * |
185 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
186 | */ |
187 | void GetHistos(Double_t lowLim, |
188 | Double_t highLim, |
189 | Dir* resp, |
190 | Dir* data, |
191 | TDirectory* out); |
192 | /** |
193 | * Do the actual unfolding and corrections |
194 | * |
195 | * @param response Response matrix |
196 | * @param measured Measured distribution |
197 | * @param mcHist MC truth distribution |
198 | * @param esdHist MC distribution after ESD event selection |
671df6c9 |
199 | * |
200 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
201 | */ |
671df6c9 |
202 | void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist); |
4f09623b |
203 | /** |
204 | * Caclulate the trigger bias correction |
205 | * |
206 | * @param mcHist MC truth distribution |
207 | * @param esdHist MC distribution after ESD event selection |
208 | * |
209 | * @return Histogram of trigger bias correction |
671df6c9 |
210 | * |
211 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
212 | */ |
213 | TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist); |
214 | /** |
215 | * Apply the trigger bias correction to the unfolded result |
216 | * |
217 | * @param unfoldBefore Unfolded measured distribution |
218 | * @param triggerBias Trigger bias correction |
219 | * @param mcHist MC truth distribution |
220 | */ |
221 | void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1* mcHist); |
222 | |
223 | |
224 | /** |
225 | * Main Loop to call eta range |
226 | * |
227 | * @param respFile File containing response matrices |
228 | * @param dataFile File containing measured distributions |
229 | * @param outFile Output file |
671df6c9 |
230 | * |
231 | * @ingroup pwglf_forward_scripts_unfoldmult |
4f09623b |
232 | */ |
233 | void UnfoldMult(const char* respFile="forward_response.root", |
234 | const char* dataFile="forward_multiplicy.root", |
235 | const char* outFile="unfolded.root") |
236 | { |
237 | TFile* output = OpenFile(outFile,false); |
238 | TFile* resp = OpenFile(respFile, true); |
239 | TFile* data = OpenFile(dataFile, true); |
240 | if (!output || !resp || !data) return; |
241 | |
242 | Dir* topResp = GetCollection(resp, "ResponseMatricesSums"); |
243 | Dir* topData = GetCollection(data, "MultSums"); |
244 | if (!topData || !topResp) return; |
245 | |
246 | Double_t limits[] = {5.1, 3.4, 3.0, 2.4, 2.0, 1.5, 1.0, 0.5, 0. }; |
247 | Double_t* limit = limits; |
248 | |
249 | while ((*limit) > 0.1) { |
250 | if ((*limit) >5.) GetHistos(-3.4, (*limit), topResp, topData, output); |
251 | else GetHistos(-(*limit), (*limit), topResp, topData, output); |
252 | limit++; |
253 | output->cd(); |
254 | } |
255 | output->Close(); |
256 | resp->Close(); |
257 | data->Close(); |
258 | } |
259 | |
260 | //____________________________________________________________________ |
261 | void GetHistos(Double_t lowLim, |
262 | Double_t highLim, |
263 | Dir* topResp, |
264 | Dir* topData, |
265 | TDirectory* out) |
266 | { |
267 | // Format the name of the bin |
268 | TString sLow( TString::Format("%+03d",int(10*lowLim))); |
269 | TString sHigh(TString::Format("%+03d",int(10*highLim))); |
270 | sLow .ReplaceAll("+", "plus"); |
271 | sLow .ReplaceAll("-", "minus"); |
272 | sHigh.ReplaceAll("+", "plus"); |
273 | sHigh.ReplaceAll("-", "minus"); |
274 | TString name = TString::Format("%s_%s", sLow.Data(), sHigh.Data()); |
275 | TDirectory* dir = out->mkdir(name); |
276 | dir->cd(); |
277 | |
278 | // Get our collections |
279 | Dir* listResp = GetCollection(topResp, name); |
280 | Dir* listData = GetCollection(topData,name); |
281 | if(!listResp || !listData) return; |
282 | |
283 | // Get the MC histograms |
284 | TH2* response = GetH2(listResp, "responseMatrix"); |
285 | TH1* mcHist = GetH1(listResp, "fMCNSD"); |
286 | TH1* esdHist = GetH1(listResp, "fESDNSD"); |
287 | if (!response || !mcHist || !esdHist) return; |
288 | |
289 | // Get the data histogram |
290 | TH1* measured = GetH1(listData, "mult"); |
291 | if(!measured) return; |
292 | |
293 | // Now do the unfolding |
294 | ProcessUnfold(response, measured, mcHist, esdHist); |
295 | } |
296 | |
297 | |
298 | |
299 | //____________________________________________________________________ |
300 | void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist) |
301 | { |
302 | Int_t nX = response->GetNbinsX(); |
303 | Int_t nY = response->GetNbinsY(); |
304 | TH1* projY = static_cast<TH1*>(response->ProjectionY("projY",1,nX,"")); |
305 | TH1* projX = static_cast<TH1*>(response->ProjectionX("projX",1,nY,"")); |
306 | projX->SetDirectory(0); |
307 | projY->SetDirectory(0); |
308 | |
309 | TH2* responseTranspose = static_cast<TH2*>(response->Clone("response")); |
310 | for(Int_t i = 1; i <= nX; i++) { |
311 | for(Int_t j = 1; j <= nY; j++) { |
312 | responseTranspose->SetBinContent(j,i, response->GetBinContent(i,j)); |
313 | responseTranspose->SetBinError(j,i, response->GetBinError(i,j)); |
314 | } |
315 | } |
316 | |
317 | RooUnfoldResponse* responseObject = new RooUnfoldResponse(projY, projX, |
318 | responseTranspose, |
319 | "", ""); |
320 | RooUnfold* unfold = new RooUnfoldBayes(responseObject, |
321 | measured, 5); |
322 | TH1* unfolded = static_cast<TH1*>(unfold->Hreco()); |
323 | TH1* triggerBias = GetTriggerBias(mcHist, esdHist); |
324 | DoCorrection(unfolded, TriggerBias, mcHist); |
325 | |
326 | Double_t scale_unf = 1/unfolded->Integral(); |
327 | unfolded->Scale(scale_unf); |
328 | unfolded->Write("Unfolded"); |
329 | |
330 | Double_t scale_raw = 1/measured->Integral(); |
331 | measured->Scale(scale_raw); |
332 | measured->Write("Raw"); |
333 | triggerBias->Write("TriggerBiasCorr"); |
334 | } |
335 | |
336 | |
337 | |
338 | //____________________________________________________________________ |
339 | TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist) |
340 | { |
341 | mcHist->Sumw2(); |
342 | esdHist->Sumw2(); |
343 | |
344 | TH1* corr = new TH1D("corr", "corr", 40, -0.5, 39.5); |
345 | for (Int_t j=1; j<corr->GetNbinsX(); j++) { |
346 | Double_t errSqr = 0.; |
347 | Double_t cESD = esdHist->GetBinContent(j); |
348 | Double_t cMS = mcHist->GetBinContent(j); |
349 | Double_t c = (cMC == 0 ? 1 : |
350 | cESD == 0 ? 1/cMC : cESD/cMC); |
351 | corr->SetBinContent(j, c); |
352 | corr->SetBinError(j, c*TMath::Sqrt(errSqr)); |
353 | } |
354 | return corr; |
355 | } |
356 | |
357 | //____________________________________________________________________ |
358 | void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1*) |
359 | { |
360 | for (Int_t k=1; k<=35; k++) { |
361 | Double_t tb = triggerBias->GetBinContent(k); |
362 | Double_t ub = unfoldBefore->GetBinContent(k); |
363 | if (tb > 1e-5 && ub > 0) { |
364 | unfoldBefore->SetBinContent(k, ub / tb); |
365 | Double_t eub = UnfoldBefore->GetBinError(k); |
366 | Double_t etb = TriggerBias->GetBinError(k); |
367 | unfoldBefore->SetBinError(k, TMath::Sqrt(TMath::Power(eub/ub)+ |
368 | TMath::Power(etb/tb))*ub/tb); |
369 | } |
370 | } |
371 | } |
372 | // |
373 | // EOF |
374 | // |
375 | |
376 | |
377 | |