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