Fix some documentation issues
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / UnfoldMult.C
CommitLineData
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 */
38TObject* 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 */
57TObject* 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 */
77TCollection* 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 */
91TCollection* 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 */
105TH1* 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 */
119TH1* 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 */
131TH2* 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 */
145TH2* 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 */
159TFile* 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 173typedef 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 */
187void 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 202void 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 */
213TH1* 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 */
221void 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 */
233void 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//____________________________________________________________________
261void 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//____________________________________________________________________
300void 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//____________________________________________________________________
339TH1* 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//____________________________________________________________________
358void 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