Fix some documentation issues
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / scripts / UnfoldMult.C
1 /**
2  * @file   UnfoldMult.C
3  * @author Valentina Zaccolo
4  * @date   Tue Mar 05 12:56:56 2013
5  * 
6  * @brief  Bayesian Method Unfolding Script
7  *
8  * @ingroup pwglf_forward_scripts
9  * @ingroup pwglf_forward_multdist
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>
23 /** 
24  * @defgroup pwglf_forward_scripts_unfoldmult UnfoldMult stuff 
25  *
26  * @ingroup pwglf_forward_multdist
27  */
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 
35  *
36  * @ingroup pwglf_forward_scripts_unfoldmult
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 
54  *
55  * @ingroup pwglf_forward_scripts_unfoldmult
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 
74  *
75  * @ingroup pwglf_forward_scripts_unfoldmult
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 
88  *
89  * @ingroup pwglf_forward_scripts_unfoldmult
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 
102  *
103  * @ingroup pwglf_forward_scripts_unfoldmult
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 
116  *
117  * @ingroup pwglf_forward_scripts_unfoldmult
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 
142  *
143  * @ingroup pwglf_forward_scripts_unfoldmult
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
156  *
157  * @ingroup pwglf_forward_scripts_unfoldmult
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 }
168 /** 
169  * A type definition to make it easier to switch to TDirectory structure 
170  *
171  * @ingroup pwglf_forward_scripts_unfoldmult
172  */
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  
184  *
185  * @ingroup pwglf_forward_scripts_unfoldmult
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
199  *
200  * @ingroup pwglf_forward_scripts_unfoldmult
201  */
202 void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist); 
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
210  *
211  * @ingroup pwglf_forward_scripts_unfoldmult
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 
230  *
231  * @ingroup pwglf_forward_scripts_unfoldmult
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