AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
UnfoldingHF.C
Go to the documentation of this file.
1 //-----------------------------------------------------//
2 // //
3 // Base Macro to evaluate pt-eta unfolding //
4 // for charm mesons //
5 // //
6 // Usage: //
7 // //
8 // 1) By default need output.root file from AliCF... //
9 // task //
10 // 2) ./output.root - used to evaluate the eff //
11 // 3) To correct, the data sample should be spitted in //
12 // 2, in this case in GetMeasuredSpectrum() connect //
13 // the second container ./output_data.root //
14 // 4) Set the number of iterations for Baysian unf //
15 // 5) Set min. Chi2 //
16 // 6) Select if unfold at Acceptance level or //
17 // after PPR cuts //
18 // 7) Gaussian error propagation assumed, may be //
19 // errors overestimated - to be studied //
20 // //
21 //-----------------------------------------------------//
22 // //
23 // a.grelli@uu.nl- Utrecht for R.Vernaut example //
24 // //
25 //-----------------------------------------------------//
26 // //
27 // NOTE: by default the macro is computing the //
28 // unfolding after acceptance step. To evaluate //
29 // after PPR switch on the option in the CF task! //
30 // and change in this macro the selection steps //
31 // //
32 //-----------------------------------------------------//
33 
34 #include "TString.h"
35 #include <iostream>
36 #include <fstream>
37 
38 void UnfoldingHF() {
39 
40  // not jet impemented the "systematics" mode!
41 
42  //TString smode(analysis_mode);
43  //TString hist(method);
44 
45  Load();
46 
47  printf("==================================================================\n");
48  printf("=========== RUNNING D MESONS UNFOLDING ==========\n");
49  printf("==================================================================\n");
50 
51  // --------------------- get variables from container ----------------------//
52 
53  AliCFDataGrid *measured = (AliCFDataGrid*) GetMeasuredSpectrum();
54  AliCFDataGrid *reconstructed = (AliCFDataGrid*) GetReconstructedSpectrum();
55  AliCFDataGrid *generated = (AliCFDataGrid*) GetGeneratedSpectrum();
56  AliCFEffGrid *efficiency = (AliCFEffGrid*) GetEfficiency();
57 
58  // the response matrix
59 
60  THnSparse *response = (THnSparse*) GetResponseMatrix();
61 
62  // -------------------------------------------------------------------------//
63  // //
64  // The HF correction framework has 12 dimensions. This is bad for //
65  // pt-eta unfolding so reduce the dimesions from 12 to 2. //
66  // //
67  //--------------------------------------------------------------------------//
68 
69 
70  THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
71  THnSparse* data = CreateGuessed(((AliCFGridSparse*)measured->GetData())->GetGrid());
72  THnSparse* Reco = CreateGuessed(((AliCFGridSparse*)reconstructed->GetData())->GetGrid()) ;
73 
74 
75  // best guess, traslate the 12 dimentions container over 2 (eta and pt)
76 
77  Int_t dim[2];
78 
79  dim[0] = 0;
80  dim[1] = 1;
81 
82  THnSparse* BestGuess = guessed ->Projection(2,dim);
83  THnSparse* DataSample = data ->Projection(2,dim);
84  THnSparse* Reconstructed = Reco ->Projection(2,dim);
85  THnSparse* SimpleEff = Reco ->Projection(2,dim);
86 
87  SimpleEff->Divide(Reconstructed,BestGuess,1.,1.);
88 
89 
90  // here I do the unfolding and I set the number of iterations and the chi2
91 
92  AliCFUnfolding unfolding("unfolding","",2,response,SimpleEff,DataSample,BestGuess);
93  unfolding.SetMaxNumberOfIterations(100);
94  unfolding.SetMaxChi2PerDOF(0.00000000005);
95  unfolding.Unfold();
96 
97  THnSparse* result = unfolding.GetUnfolded();
98 
99 
100  TH2D* h_gen = generated ->Project(1,0);
101  TH2D* h_meas = DataSample ->Projection(1,0);
102  TH2D* h_reco = Reconstructed ->Projection(1,0);
103  TH2D* h_guessed = BestGuess ->Projection(1,0);
104  TH2D* h_eff = SimpleEff ->Projection(1,0);
105  TH2D* h_unf = result ->Projection(1,0);
106 
107  // Apply simple efficiency correction
108 
109 
110  TH2D* simpleC = (TH2D*)h_meas->Clone();
111  simpleC->Divide(h_eff);
112 
113  TCanvas * canvas3 = new TCanvas("Unfolded efficiency map","canvas 3",1000,700);
114 
115  canvas3->cd();
116 
117  TH2D* ratio = (TH2D*)h_unf->Clone();
118  ratio->SetTitle("corrected using unfolding");
119  ratio->Divide(h_unf,h_guessed,1,1);
120 
121  ratio->Draw("cont4z");
122 
123  TCanvas * canvas4 = new TCanvas("Simple Efficiency map","",1000,700);
124 
125  canvas4->cd();
126  h_meas->SetTitle("simple correction");
127  h_meas->Divide(h_eff);
128  h_meas->Divide(h_guessed);
129  h_meas->Draw("cont4z");
130 
131  TCanvas * distribution2 = new TCanvas("dist2","",1000,700);
132 
133  distribution2->cd();
134 
135  TH1D* gen2 = generated->Project(1); // generated pt
136  gen2->SetLineColor(1);
137  gen2->SetTitle("D0 #eta distribution");
138  gen2->GetXaxis()->SetTitle("eta");
139  gen2->SetMarkerStyle(21);
140  gen2->SetMarkerSize(1.);
141  gen2->SetMarkerColor(1);
142  gen2->Draw();
143 
144  TH1D* meas2 = measured->Project(1); // generated pt
145  meas2->SetTitle("measurements");
146  meas2->SetLineColor(2);
147  meas2->SetMarkerStyle(22);
148  meas2->SetMarkerSize(1.);
149  meas2->SetMarkerColor(2);
150  meas2->DrawCopy("same");
151 
152  TH1D* unf2 = result->Projection(1); // generated pt
153  unf2->SetTitle("measurements");
154  unf2->SetLineColor(3);
155  unf2->SetMarkerStyle(23);
156  unf2->SetMarkerSize(1.);
157  unf2->SetMarkerColor(3);
158  unf2->DrawCopy("SAME");
159 
160  corr2 = new TH1D();
161  corr2->Sumw2();
162  corr2 = simpleC->ProjectionY();
163  corr2->SetTitle("measurements");
164  corr2->SetLineColor(4);
165  corr2->SetMarkerStyle(24);
166  corr2->SetMarkerSize(1.);
167  corr2->SetMarkerColor(4);
168  corr2->Sumw2();
169  corr2->DrawCopy("SAME");
170 
171  leg = new TLegend(0.1,0.7,0.48,0.9);
172  leg->SetHeader("Distributions");
173  leg->AddEntry(gen2,"Generated PYTHIA","l");
174  leg->AddEntry(meas2,"Reconstructed","l");
175  leg->AddEntry(corr2,"Corrected for eff","l");
176  leg->AddEntry(unf2,"Unfolded","l");
177 
178  leg->Draw();
179 
180  TCanvas * distribution3 = new TCanvas("dist3","",1000,700);
181 
182  distribution3->cd();
183 
184 
185  TH1D* gen3 = generated->Project(0); // generated pt
186  Int_t Entries = gen3->GetEntries();
187  gen3->SetLineColor(1);
188  gen3->SetTitle("D0 pt distribution");
189  gen3->GetXaxis()->SetTitle("pt (GeV/c)");
190  gen3->SetMarkerStyle(21);
191  gen3->SetMarkerSize(1.);
192  gen3->SetMarkerColor(1);
193  gen3->Draw();
194 
195 
196  distribution3->cd();
197  TH1D* meas3 = measured->Project(0); // generated pt
198  meas3->SetTitle("measurements");
199  meas3->SetLineColor(2);
200  meas3->SetMarkerStyle(22);
201  meas3->SetMarkerSize(1.);
202  meas3->SetMarkerColor(2);
203  meas3->DrawCopy("SAME");
204 
205  TH1D* unf3 = result->Projection(0); //
206  unf3->SetTitle("unfolded");
207  unf3->SetLineColor(3);
208  unf3->SetMarkerStyle(23);
209  unf3->SetMarkerSize(1.);
210  unf3->SetMarkerColor(3);
211  unf3->Sumw2();
212  unf3->DrawCopy("SAME");
213 
214  TH1D* corr3 = simpleC->ProjectionX();
215  corr3->SetTitle("corrected");
216  corr3->SetLineColor(4);
217  corr3->SetMarkerStyle(24);
218  corr3->SetMarkerSize(1.);
219  corr3->SetMarkerColor(4);
220  corr3->Sumw2();
221  corr3->DrawCopy("SAME");
222 
223  leg = new TLegend(0.1,0.7,0.48,0.9);
224  leg -> SetHeader("Distributions");
225  leg -> AddEntry(gen3,"Generated PYTHIA","l");
226  leg -> AddEntry(meas3,"Reconstructed","l");
227  leg -> AddEntry(corr3,"Corrected for eff","l");
228  leg -> AddEntry(unf3,"Unfolded","l");
229 
230  leg->Draw();
231 
232  return;
233 }
234 
235 // ====================================================================
236 
237 
238 void Load(){
239  gSystem->Load("libANALYSIS");
240  gSystem->Load("libCORRFW");
241  gStyle->SetPalette(1);
242  gStyle->SetOptStat(0);
243 }
244 
245 //___________________________________________________________-
246 
248  TFile * f = TFile::Open("output.root","read");
249  AliCFContainer* c = (AliCFContainer*)f->Get("container");
250  AliCFDataGrid* data1 = new AliCFDataGrid("data1","",*c);
251  data1->SetMeasured(3);
252  return data1;
253 }
255  TFile * f = TFile::Open("output.root","read");
256  AliCFContainer* c = (AliCFContainer*)f->Get("container");
257  AliCFDataGrid* data2 = new AliCFDataGrid("data2","",*c);
258  data2->SetMeasured(3);
259 
260  return data2;
261 }
263  TFile * f = TFile::Open("output.root","read");
264  AliCFContainer* c = (AliCFContainer*)f->Get("container");
265  AliCFDataGrid* data3 = new AliCFDataGrid("data3","",*c);
266  data3->SetMeasured(1);
267  return data3;
268 }
270  TFile * f = TFile::Open("output.root","read");
271  AliCFContainer* c = (AliCFContainer*)f->Get("container");
272  AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
273  eff->CalculateEfficiency(3,1);
274  return eff;
275 }
276 THnSparse* GetResponseMatrix() {
277  TFile * f = TFile::Open("output.root","read");
278  THnSparse* h = (THnSparse*)f->Get("correlation");
279  return h;
280 }
281 
282 THnSparse* GetResponseMatrixPPR() {
283  TFile * f = TFile::Open("output.root","read");
284  THnSparse* h = (THnSparse*)f->Get("correlationPPR");
285  return h;
286 }
287 
288 THnSparse* CreateGuessed(const THnSparse* h) {
289  THnSparse* guessed = (THnSparse*) h->Clone();
290  return guessed ;
291 }
292 
293 
TObject * GetMeasuredSpectrum()
Definition: UnfoldingHF.C:247
TObject * GetGeneratedSpectrum()
Definition: UnfoldingHF.C:262
TSystem * gSystem
TCanvas * c
Definition: TestFitELoss.C:172
THnSparse * CreateGuessed(const THnSparse *h)
Definition: UnfoldingHF.C:288
void Load()
Definition: UnfoldingHF.C:238
int Int_t
Definition: External.C:63
TObject * GetReconstructedSpectrum()
Definition: UnfoldingHF.C:254
THnSparse * GetResponseMatrix()
Definition: UnfoldingHF.C:276
Definition: External.C:228
Definition: External.C:212
TObject * GetEfficiency()
Definition: UnfoldingHF.C:269
THnSparse * GetResponseMatrixPPR()
Definition: UnfoldingHF.C:282
void UnfoldingHF()
Definition: UnfoldingHF.C:38