Attachment 'CalibrationFunctions.C'
Download 1 //###################################################################
2 // Determine Reconstruction efficiency and Energy calibration
3 // functions to be applied to theoretical background expectation (RecEff)
4 // and Data samples (Energy calibration)
5 //
6 // Reconstruction efficiency: Look if a ProcessedReco photon
7 // connected to MC photon via MCTruth is found within 0.1 rad (~5.7 deg)
8 //
9 // ErecVsEgen: Ratio of highest Processed Photon energy to MC energy.
10 // Processed photon has to be connected to MC photon
11 // and within 0.1 rad. Only MC photons not converted before
12 // ECAL (or before end of tracking) are taken into account
13 // Due to merging there should be only one Processed photon
14 // connected to MC photon in 0.1 rad
15 //###################################################################
16 //###################################################################
17
18 # include "GlobalVariables.C"
19
20 // ########### CUTS
21
22 double energy_cut = 2.0;
23 double cosTheta_cut[2] = {0.0,0.995};
24 int conv_cut = 1;
25
26 void CalibrationFunctions(double e_min, double cos_min, double cos_max, int conv){
27
28 energy_cut = e_min;
29 cosTheta_cut[0] = cos_min;
30 cosTheta_cut[1] = cos_max;
31 conv_cut = conv;
32
33 GammaReconstructonEfficiency();
34 ErecVsEgen();
35
36 gROOT->ProcessLine(".q");
37
38 // return;
39 }
40
41
42 void GammaReconstructonEfficiency(){
43
44 printf("Reconstruction efficiencies\n");
45
46 // ############ OPEN OUTPUTFILE
47 TFile* f_out = new TFile("../data/Auxiliary/ReconstructionEfficiencies.root","RECREATE");
48
49 // ############ RECONSTRUCTION EFFICIENCY HISTOGRAMS IN COS(THETA), THETA, E AND PT
50
51
52 TCanvas* c1 = new TCanvas("c1","Reconstruction efficiency in cos(#Theta)");
53 TCanvas* c2 = new TCanvas("c2","Reconstruction efficiency in #Theta");
54 TCanvas* c3 = new TCanvas("c3","Reconstruction efficiency in Energy");
55 TCanvas* c4 = new TCanvas("c4","Reconstruction efficiency in PT");
56 TCanvas* c5 = new TCanvas("c5","Reconstruction efficiency in E and cos(#Theta)");
57
58 TProfile* p1 = new TProfile("RecEff_Cos","Reconstruction efficiency;cos(#Theta_{gen,#gamma});Efficiency",100,cosTheta_cut[0],cosTheta_cut[1]);
59 TProfile* p2 = new TProfile("RecEff_Theta","Reconstruction efficiency;#Theta_{gen,#gamma} [rad];Efficiency",100,0.0,3.2);
60 TProfile* p3 = new TProfile("RecEff_E","Reconstruction efficiency;E_{gen,#gamma} [GeV];Efficiency",250,0,250);
61 TProfile* p4 = new TProfile("RecEff_PT","Reconstruction efficiency;P_{T,#gamma} [GeV];Efficiency",250,0,250);
62 TProfile2D* p5 = new TProfile2D("RecEff_E_Cos","Reconstruction efficiency;E_{gen,#gamma} [GeV]; cos(#Theta_{gen,#gamma});Efficiency",25,0,250,10,cosTheta_cut[0],cosTheta_cut[1]);
63
64
65 p1->GetYaxis()->SetRangeUser(0.8,1.05);
66 p2->GetYaxis()->SetRangeUser(0.8,1.05);
67 p3->GetYaxis()->SetRangeUser(0.8,1.05);
68 p4->GetYaxis()->SetRangeUser(0.8,1.05);
69 p5->GetZaxis()->SetRangeUser(0.8,1.05);
70
71
72 // ########## OPEN INPUT FILES
73
74 TChain* tree = new TChain("Tree");
75 tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
76 tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root");
77
78 LoadTree(tree);
79
80 TChain* proctree = new TChain("ProcTree");
81 proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
82 proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root");
83 LoadPreProcTree(proctree);
84
85
86 // ########## FILL HISTOGRAMS
87
88 int nEntries = proctree->GetEntries();
89 // nEntries = 1e5;
90
91 for ( int i(0); i < nEntries; i++ ){
92 tree->GetEntry(i);
93 proctree->GetEntry(i);
94 for ( int j(0); j < nMCParticles; j++ ){
95 // Check if MC is photon and apply cuts
96 if ( MC_type[j] == 22 && MC_isReal[j] && MC_energy[j] > energy_cut &&
97 TMath::Abs(MC_cosTheta[j]) > cosTheta_cut[0] && TMath::Abs(MC_cosTheta[j]) < cosTheta_cut[1] ){
98 int counter(0);
99 for ( int k(0); k < nProcessedRecos; k++ ){
100 if ( ProcessedReco_type[k] == 22 ){
101 double angle = TMath::ACos((MC_px[j]*ProcessedReco_px[k]+MC_py[j]*ProcessedReco_py[k]
102 +MC_pz[j]*ProcessedReco_pz[k])/(MC_absP[j]*ProcessedReco_absP[k]));
103 if (angle < 0.0) cout << "ANGLE < 0.0" << endl;
104
105 if ( MC_id[j] == ProcessedRecoRelatedToId[k] && angle < 0.1){
106 counter++;
107 break; // break as soon as one PFlow photon is found
108 }
109 }
110 }
111
112 double PT = MC_energy[j]*TMath::Sin(TMath::ACos(MC_cosTheta[j]));
113
114 // Fill for each "real" MC photon, counter is either 0 or 1
115 p1->Fill(TMath::Abs(MC_cosTheta[j]),counter);
116 p2->Fill(TMath::ACos(MC_cosTheta[j]),counter);
117 p3->Fill(MC_energy[j],counter);
118 p4->Fill(PT,counter);
119 p5->Fill(MC_energy[j],TMath::Abs(MC_cosTheta[j]),counter);
120 }
121 }
122 }
123
124
125 // ########## DRAW
126
127
128
129 TLatex latex;
130 latex.SetNDC();
131 latex.SetTextSize(0.05);
132 latex.SetTextAlign(13);
133
134 c1->cd();
135 p1->DrawCopy();
136 latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p1->GetMean(2)));
137
138 c2->cd();
139 p2->DrawCopy();
140 c2->cd();
141 latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p2->GetMean(2)));
142
143 c3->cd();
144 p3->DrawCopy();
145 latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p3->GetMean(2)));
146
147 c4->cd();
148 p4->DrawCopy();
149 latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p4->GetMean(2)));
150
151 c5->cd();
152 p5->DrawCopy();
153
154
155 // ############ PRINT TO FILE
156
157 f_out->cd();
158
159 p1->Write();
160 p2->Write();
161 p3->Write();
162 p4->Write();
163 p5->Write();
164
165 c1->Print(Form("Reconstruction_Efficiency_cosTheta.eps"));
166 c1->Print(Form("Reconstruction_Efficiency_cosTheta.png"));
167
168
169 c2->Print(Form("Reconstruction_Efficiency_Theta.eps"));
170 c2->Print(Form("Reconstruction_Efficiency_Theta.png"));
171
172 c3->Print(Form("Reconstruction_Efficiency_Energy.eps"));
173 c3->Print(Form("Reconstruction_Efficiency_Energy.png"));
174
175 c4->Print(Form("Reconstruction_Efficiency_PT.eps"));
176 c4->Print(Form("Reconstruction_Efficiency_PT.png"));
177
178 c5->Print(Form("Reconstruction_Efficiency_Energy_cosTheta.eps"));
179 c5->Print(Form("Reconstruction_Efficiency_Energy_cosTheta.png"));
180
181 //
182 f_out->Close();
183 return;
184 }
185
186
187 // #########################################################################################################
188 // #########################################################################################################
189
190
191
192 void ErecVsEgen(){
193
194 printf("Energy calibration\n");
195
196 // ######## OPEN OUTPUTFILE
197
198 TFile* f_out = new TFile("../data/Auxiliary/Energycalibration.root","RECREATE");
199
200
201 char str[128];
202
203 // ######## CUTS
204
205 double z_cut,r_cut;
206
207 if (conv_cut == 1){
208 sprintf(str,"InECAL");
209 z_cut = 2450.0;
210 r_cut = 1843.3;
211 }
212 if (conv_cut == 2){
213 sprintf(str,"AfterTracking");
214 z_cut = 2246.0; // cut at tpc dimensions
215 r_cut = 1739.0;
216 }
217
218 // ############ ENERGY RATIO HISTOGRAMS IN COS(THETA), THETA, E AND PT
219
220 TCanvas* c1 = new TCanvas("c2","Ratio E_rec/E_gen vs cos(Theta_gen)");
221 TCanvas* c2 = new TCanvas("c3","Ratio E_rec/E_gen vs Theta_gen");
222 TCanvas* c3 = new TCanvas("c1","Ratio E_rec/E_gen vs E_gen");
223 TCanvas* c4 = new TCanvas("c4","Ratio E_rec/E_gen vs PT_gen");
224 TCanvas* c5 = new TCanvas("c5","Ratio E_rec/E_gen vs PT_gen");
225
226 TProfile* p1 = new TProfile("ErecEgen_Cos","Ratio E_rec/E_gen vs cos(Theta_gen);cos(#Theta_{gen});#bar{E_{rec}/E_{gen}}",150,0,1);
227 TProfile* p2 = new TProfile("ErecEgen_Theta","Ratio E_rec/E_gen vs Theta_gen;#Theta_{gen};#bar{E_{rec}/E_{gen}}",150,0,3.2);
228 TProfile* p3 = new TProfile("ErecEgen_E","Ratio E_rec/E_gen vs E_gen;E_{gen} [GeV];#bar{E_{rec}/E_{gen}}",250,0,250);
229 TProfile* p4 = new TProfile("ErecEgen_PT","Ratio E_rec/E_gen vs PT_gen;P_{T} [GeV];#bar{E_{rec}/E_{gen}}",250,0,250);
230 TProfile2D* p5 = new TProfile2D("ErecEgen_E_Cos","Ratio E_rec/E_gen vs E and Cos;E_{gen} [GeV];cos(#Theta_{gen});#bar{E_{rec}/E_{gen}}",25,0,250,10,0.0,1.0);
231
232 p1->GetYaxis()->SetRangeUser(0.8,1.05);
233 p2->GetYaxis()->SetRangeUser(0.8,1.05);
234 p3->GetYaxis()->SetRangeUser(0.8,1.05);
235 p4->GetYaxis()->SetRangeUser(0.8,1.05);
236 p5->GetZaxis()->SetRangeUser(0.8,1.05);
237
238
239 // ########## OPEN INPUT FILES
240
241 TChain* tree = new TChain("Tree");
242 tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
243 tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root");
244
245 LoadTree(tree);
246
247 TChain* proctree = new TChain("ProcTree");
248 proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
249 proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root");
250 LoadPreProcTree(proctree);
251
252
253
254 // ####### FILL HISTOGRAMS
255
256 int nEntries = tree->GetEntries();
257 //nEntries = 1e5;
258
259 for ( int i(0); i < nEntries; i++ ){
260 tree->GetEntry(i);
261 proctree->GetEntry(i);
262 for ( int j(0); j < nMCParticles; j++ ){
263 if ( MC_type[j] == 22 && MC_isReal[j] && MC_energy[j] > energy_cut &&
264 TMath::Abs(MC_cosTheta[j]) > cosTheta_cut[0] && TMath::Abs(MC_cosTheta[j]) < cosTheta_cut[1] ){
265
266 // cut conversions
267
268 double z = TMath::Abs(MC_endpoint_z[j]);
269 double r = TMath::Sqrt(MC_endpoint_x[j]*MC_endpoint_x[j] + MC_endpoint_y[j]*MC_endpoint_y[j]);
270
271 if ( conv_cut && (z < z_cut && r < r_cut) )
272 continue;
273
274 // find max energy reco gamma
275 double maxEnergy(0.0);
276 int Recoid(-1);
277 for ( int k(0); k < nProcessedRecos; k++ ){ // Find PFlow photon with highest energy for MC photon
278 if ( ProcessedReco_type[k] == 22 && MC_id[j] == ProcessedRecoRelatedToId[k] && ProcessedReco_energy[k] > maxEnergy){
279 maxEnergy = ProcessedReco_energy[k];
280 Recoid = k;
281 }
282 }
283 if (Recoid != -1){
284 double angle = TMath::ACos((MC_px[j]*ProcessedReco_px[Recoid]+MC_py[j]*ProcessedReco_py[Recoid]
285 +MC_pz[j]*ProcessedReco_pz[Recoid])/(MC_absP[j]*ProcessedReco_absP[Recoid]));
286 if (angle < 0.1){
287 double ratio = ProcessedReco_energy[Recoid]/MC_energy[j];
288 p1->Fill(TMath::Abs(MC_cosTheta[j]),ratio);
289 p2->Fill(TMath::ACos(MC_cosTheta[j]),ratio);
290 p3->Fill(MC_energy[j],ratio);
291 p4->Fill(TMath::Sin(TMath::ACos(MC_cosTheta[j]))*MC_energy[j],ratio);
292 p5->Fill(MC_energy[j],MC_cosTheta[j],ratio);
293 }
294 }
295 }
296 }
297 }
298
299 // ######## FIT
300
301 TF1* func = new TF1("EnergyCalibration_Func",fitfunc,0.0,1.0,13);
302
303 func->SetParameter(0,1.0);
304
305 func->SetParameter(1,-0.1);
306 func->SetParameter(2,0.23);
307 func->SetParameter(3,0.01);
308
309 func->SetParameter(4,-0.1);
310 func->SetParameter(5,0.58);
311 func->SetParameter(6,0.01);
312
313 func->SetParameter(7,-0.1);
314 func->SetParameter(8,0.77);
315 func->SetParameter(9,0.01);
316
317 func->SetParameter(10,-0.1);
318 func->SetParameter(11,0.98);
319 func->SetParameter(12,0.01);
320
321 // ########## DRAW
322
323 c1->cd();
324 p1->Fit(func,"R","");
325 p1->DrawCopy();
326
327 c2->cd();
328 p2->DrawCopy();
329
330 c3->cd();
331 p3->DrawCopy();
332
333 c4->cd();
334 p4->DrawCopy();
335
336 c5->cd();
337 p5->DrawCopy();
338 // ########### PRINT
339
340 f_out->cd();
341
342 p1->Write();
343 p2->Write();
344 p3->Write();
345 p4->Write();
346 p5->Write();
347 func->Write();
348
349 f_out->Close();
350
351
352 c1->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_cosTheta_%s.eps",str));
353 c1->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_cosTheta_%s.png",str));
354 //c1->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_cosTheta_%s.root",str));
355
356 c2->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Theta_%s.eps",str));
357 c2->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Theta_%s.png",str));
358
359 c3->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_%s.eps",str));
360 c3->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_%s.png",str));
361
362 c4->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_PT_%s.eps",str));
363 c4->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_PT_%s.png",str));
364
365 c5->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_cosTheta_%s.eps",str));
366 c5->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_cosTheta_%s.png",str));
367 return;
368 }
369
370
371 // ########### FITFUNCTION
372
373 double constant(double* x, double* par){
374 return par[0];
375 }
376
377 double gaussdip(double* x, double* par){
378 return par[0]/(sqrt(2*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2*par[2]*par[2]));
379 }
380
381 double fitfunc(double* x, double* par){
382 return constant(x,par) + gaussdip(x,&par[1]) + gaussdip(x,&par[4]) + gaussdip(x,&par[7]) + gaussdip(x,&par[10]);
383 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.