1 //-----------------------------------------------------//
3 // Base Macro to evaluate pt-eta unfolding //
8 // 1) By default need output.root file from AliCF... //
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 //
18 // 7) Gaussian error propagation assumed, may be //
19 // errors overestimated - to be studied //
21 //-----------------------------------------------------//
23 // a.grelli@uu.nl- Utrecht for R.Vernaut example //
25 //-----------------------------------------------------//
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 //
32 //-----------------------------------------------------//
40 // not jet impemented the "systematics" mode!
42 //TString smode(analysis_mode);
43 //TString hist(method);
47 printf("==================================================================\n");
48 printf("=========== RUNNING D MESONS UNFOLDING ==========\n");
49 printf("==================================================================\n");
51 // --------------------- get variables from container ----------------------//
53 AliCFDataGrid *measured = (AliCFDataGrid*) GetMeasuredSpectrum();
54 AliCFDataGrid *reconstructed = (AliCFDataGrid*) GetReconstructedSpectrum();
55 AliCFDataGrid *generated = (AliCFDataGrid*) GetGeneratedSpectrum();
56 AliCFEffGrid *efficiency = (AliCFEffGrid*) GetEfficiency();
58 // the response matrix
60 THnSparse *response = (THnSparse*) GetResponseMatrix();
62 // -------------------------------------------------------------------------//
64 // The HF correction framework has 12 dimensions. This is bad for //
65 // pt-eta unfolding so reduce the dimesions from 12 to 2. //
67 //--------------------------------------------------------------------------//
70 THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
71 THnSparse* data = CreateGuessed(((AliCFGridSparse*)measured->GetData())->GetGrid());
72 THnSparse* Reco = CreateGuessed(((AliCFGridSparse*)reconstructed->GetData())->GetGrid()) ;
75 // best guess, traslate the 12 dimentions container over 2 (eta and pt)
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);
87 SimpleEff->Divide(Reconstructed,BestGuess,1.,1.);
90 // here I do the unfolding and I set the number of iterations and the chi2
92 AliCFUnfolding unfolding("unfolding","",2,response,SimpleEff,DataSample,BestGuess);
93 unfolding.SetMaxNumberOfIterations(100);
94 unfolding.SetMaxChi2PerDOF(0.00000000005);
97 THnSparse* result = unfolding.GetUnfolded();
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);
107 // Apply simple efficiency correction
110 TH2D* simpleC = (TH2D*)h_meas->Clone();
111 simpleC->Divide(h_eff);
113 TCanvas * canvas3 = new TCanvas("Unfolded efficiency map","canvas 3",1000,700);
117 TH2D* ratio = (TH2D*)h_unf->Clone();
118 ratio->SetTitle("corrected using unfolding");
119 ratio->Divide(h_unf,h_guessed,1,1);
121 ratio->Draw("cont4z");
123 TCanvas * canvas4 = new TCanvas("Simple Efficiency map","",1000,700);
126 h_meas->SetTitle("simple correction");
127 h_meas->Divide(h_eff);
128 h_meas->Divide(h_guessed);
129 h_meas->Draw("cont4z");
131 TCanvas * distribution2 = new TCanvas("dist2","",1000,700);
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);
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");
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");
162 corr2 = simpleC->ProjectionY();
163 corr2->SetTitle("measurements");
164 corr2->SetLineColor(4);
165 corr2->SetMarkerStyle(24);
166 corr2->SetMarkerSize(1.);
167 corr2->SetMarkerColor(4);
169 corr2->DrawCopy("SAME");
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");
180 TCanvas * distribution3 = new TCanvas("dist3","",1000,700);
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);
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");
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);
212 unf3->DrawCopy("SAME");
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);
221 corr3->DrawCopy("SAME");
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");
235 // ====================================================================
239 gSystem->Load("libANALYSIS");
240 gSystem->Load("libCORRFW");
241 gStyle->SetPalette(1);
242 gStyle->SetOptStat(0);
245 //___________________________________________________________-
247 TObject* GetMeasuredSpectrum() {
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);
254 TObject* GetReconstructedSpectrum() {
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);
262 TObject* GetGeneratedSpectrum() {
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);
269 TObject* GetEfficiency() {
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);
276 THnSparse* GetResponseMatrix() {
277 TFile * f = TFile::Open("output.root","read");
278 THnSparse* h = (THnSparse*)f->Get("correlation");
282 THnSparse* GetResponseMatrixPPR() {
283 TFile * f = TFile::Open("output.root","read");
284 THnSparse* h = (THnSparse*)f->Get("correlationPPR");
288 THnSparse* CreateGuessed(const THnSparse* h) {
289 THnSparse* guessed = (THnSparse*) h->Clone();