Attachment 'Selection.C'
Download 1 # include <string>
2 # include "Functions.C"
3
4
5
6
7
8 void Selection(int Ng_min, double e_min, double e_max, double cos_min, double cos_max, double pt_max, double evis_max, double echar_max){
9
10
11 printf("Selection\n");
12
13 // ############ CUTFLOW OUTPUT FILE
14 TFile* f_cutflow = new TFile("../data/Auxiliary/SelectionFlow.root","RECREATE");
15
16
17 // ###############################
18 TFile* f = 0;
19 ifstream txtfile("../steeringfiles/ILD00_Calibrated_RootFiles.txt");
20 char name[256];
21 string newname;
22
23 // ############################# FOR EACH ROOT FILE
24
25 while (!txtfile.eof()){
26
27 txtfile.getline(name,256,'\n');
28 printf("%s\n",name);
29
30 // ############################# OPEN ORIGINAL FILE, LOAD TREES
31
32 f = new TFile(name,"OPEN");
33 if (f->IsZombie()) continue;
34
35 TTree* tree = (TTree*)f->Get("Tree");
36 TTree* globaltree = (TTree*)f->Get("GlobalParameters");
37 TTree* calibratedtree = (TTree*)f->Get("CalibratedTree");
38 TTree* bcaltree = (TTree*)f->Get("BcalTree");
39
40 LoadTree(tree);
41 LoadGlobalTree(globaltree);
42 LoadCalibratedTree(calibratedtree);
43 LoadBcalTree(bcaltree);
44
45 tree->GetEntry(0);
46 globaltree->GetEntry(0);
47
48 TTree* xsectree = 0;
49 // if ( !strcmp(Process,"n1n1a") || !strcmp(Process,"n2n2a") || !strcmp(Process,"n3n3a")){
50 // xsectree = (TTree*)f->Get("XSecTree");
51 // LoadXSecTree(xsectree);
52 // cout << Process <<endl;
53 // }
54
55 // ########################### CREATE CUTFLOW HISTOGRAMS
56
57 char histoname[128];
58 double ft = (double)NumberOfEventsInFile/NumberOfEventsInRun;
59 f_cutflow->cd();
60
61
62 // Ngamma min
63
64 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut1_E",Process,P_e,P_p,MCLumi*ft);
65 TH1D* h1_E = new TH1D(histoname,histoname,300,0.0,300);
66 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut1_Cos",Process,P_e,P_p,MCLumi*ft);
67 TH1D* h1_cos = new TH1D(histoname,histoname,100,-1.0,1.0);
68
69 // phase space
70
71 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut2_E",Process,P_e,P_p,MCLumi*ft);
72 TH1D* h2_E = new TH1D(histoname,histoname,300,0.0,300);
73 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut2_Cos",Process,P_e,P_p,MCLumi*ft);
74 TH1D* h2_cos = new TH1D(histoname,histoname,100,-1.0,1.0);
75
76 // pt
77
78 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut5_E",Process,P_e,P_p,MCLumi*ft);
79 TH1D* h5_E = new TH1D(histoname,histoname,300,0.0,300);
80 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut5_Cos",Process,P_e,P_p,MCLumi*ft);
81 TH1D* h5_cos = new TH1D(histoname,histoname,100,-1.0,1.0);
82
83
84 // E vis
85
86 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut6_E",Process,P_e,P_p,MCLumi*ft);
87 TH1D* h6_E = new TH1D(histoname,histoname,300,0.0,300);
88 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut6_Cos",Process,P_e,P_p,MCLumi*ft);
89 TH1D* h6_cos = new TH1D(histoname,histoname,100,-1.0,1.0);
90
91 // bhabha cut
92
93 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut8_E",Process,P_e,P_p,MCLumi*ft);
94 TH1D* h8_E = new TH1D(histoname,histoname,300,0.0,300);
95 sprintf(histoname,"%s_Pe%.0fPp%.0f_lumi%.1f_Cut8_Cos",Process,P_e,P_p,MCLumi*ft);
96 TH1D* h8_cos = new TH1D(histoname,histoname,100,-1.0,1.0);
97
98 // ########################## CREATE PROCESS SPECIFIC OUTPUT FILES,
99 // ########## ADD SELECTION TREE
100
101
102 newname = name;
103 newname = (newname.substr(newname.find("ILD"))).replace(6,10,"Selection");
104 newname.insert(0,"../data/Selection/");
105
106
107 printf("%s\n",newname.c_str());
108
109 TFile* sel = new TFile(newname.c_str(),"RECREATE");
110 TTree* sel_tree = new TTree("Tree","Tree");
111 TTree* sel_calibratedtree = new TTree("CalibratedTree","CalibratedTree");
112 TTree* sel_gammatree = new TTree("SelPhotonTree","SelPhotonTree");
113 TTree* sel_globaltree = globaltree->CloneTree();
114 TTree* sel_bcaltree = new TTree("BcalTree","BcalTree");
115
116 SetTreeBranches(sel_tree);
117 SetCalibratedTreeBranches(sel_calibratedtree);
118 SetSelPhotonTreeBranches(sel_gammatree);
119 SetBcalTreeBranches(sel_bcaltree);
120
121 // TTree* sel_xsectree = 0;
122 // if ( !strcmp(Process,"n1n1a") || !strcmp(Process,"n2n2a") || !strcmp(Process,"n3n3a")){
123 // sel_xsectree = new TTree("XSecTree","XSecTree");
124 // SetXSecTreeBranches(sel_xsectree);
125 // }
126
127
128 // ############# BOOK KEEPING
129
130
131 double n_before = 0.0;
132 double n_reco = 0.0;
133 double n_after_kinetic_cuts = 0.0;
134 double n_after_all_cuts = 0.0;
135 double n_empty = 0.0; // no reco photons in event
136
137 // ################ LOOP OVER EVENTS
138
139 int n_entries = tree->GetEntries();
140 //n_entries = 1e1;
141
142 for (int j(0); j < n_entries; j++){
143 if (j%10000 == 0) printf("event %i\n",j);
144
145 int n_photons = 0; // number of photons in event
146
147 tree->GetEntry(j);
148 calibratedtree->GetEntry(j);
149 bcaltree->GetEntry(j);
150 // if (!strcmp(Process,"n1n1a") || !strcmp(Process,"n2n2a") || !strcmp(Process,"n3n3a")){
151 // xsectree->GetEntry(j);
152 // }
153
154
155 // ######### INDEX OF MOST ENERGETIC PHOTON
156 int index = GetHighEnergyCalibratedPhotonWOBeamCal(n_photons);
157
158 n_before++;
159
160 // #################### APPLY CUTS
161
162 if (index == -1){ // just make sure not to run into unallocated mem in next test
163 continue;
164 n_empty++;
165 }
166 n_reco++;
167
168 double Evis = 0.0;
169 double p_tracks = 0.0;
170 double pT_tracks = 0.0;
171 double E_neut = 0.0;
172 double E_neut_m_E_gamma = 0.0;
173 double E_char = 0.0;
174
175 // -----------------------------------------------
176 double track_pt;
177 double track_p ;
178 double track_px;
179 double track_py;
180 double track_pz;
181 // -----------------------------------------------
182
183 // calibrated tree
184
185 for (int k(0); k < nCalibratedRecos; k++){
186 if (!CalibratedReco_isBCal[k]){
187 Evis += CalibratedReco_energy[k];
188 if (CalibratedReco_charge[k] != 0.0){
189 E_char += CalibratedReco_energy[k];
190 p_tracks += CalibratedReco_absP[k];
191 pT_tracks += sqrt(CalibratedReco_px[k]*CalibratedReco_px[k] + CalibratedReco_py[k]*CalibratedReco_py[k]);
192 }
193 else {
194 E_neut += CalibratedReco_energy[k];
195 }
196 }
197 }
198 if (index != -1){
199 Evis = Evis - CalibratedReco_energy[index];
200 }
201
202 if(Evis == 0.0) Evis = -3.0;
203 if(p_tracks == 0.0) p_tracks = -3.0;
204 if(pT_tracks == 0.0) pT_tracks = -3.0;
205 if(E_char == 0.0) E_char = -3.0;
206 if(E_neut == 0.0) E_neut = -3.0;
207 E_neut_m_E_gamma = E_neut;
208 if (index != -1){
209 E_neut_m_E_gamma = E_neut - CalibratedReco_energy[index];
210 }
211 if(E_neut_m_E_gamma == 0.0) E_neut_m_E_gamma = -3.0;
212
213 int highpttrack = 0;
214 for (int l(0); l < nTracks; l++){
215 track_pt = 3*4*0.0001*1/Track_omega[l];
216 if(track_pt<0.) track_pt = -1*track_pt;
217 if (track_pt > pt_max) highpttrack = 1;
218 }
219
220
221 // CUT ONE (N_GAMMA_MIN)
222 h1_E->Fill(CalibratedReco_energy[index]); // one photon required, reco eff
223 h1_cos->Fill(CalibratedReco_cosTheta[index]);
224
225 bool pass = false;
226
227 // PHASE SPACE
228 if (CalibratedReco_energy[index] > e_min && CalibratedReco_energy[index] < e_max &&
229 TMath::Abs(CalibratedReco_cosTheta[index]) < cos_max && TMath::Abs(CalibratedReco_cosTheta[index]) > cos_min){
230
231 h2_E->Fill(CalibratedReco_energy[index]); // events after phase space cuts
232 h2_cos->Fill(CalibratedReco_cosTheta[index]);
233
234 n_after_kinetic_cuts++;
235
236 // PT
237 if (!highpttrack){
238 // if (pT_tracks <= pt_max){
239
240 h5_E->Fill(CalibratedReco_energy[index],1.0); // events after pt
241 h5_cos->Fill(CalibratedReco_cosTheta[index],1.0);
242
243 // Evis - Egamma
244 if (Evis < evis_max){
245
246 h6_E->Fill(CalibratedReco_energy[index],1.0); // events after Evis
247 h6_cos->Fill(CalibratedReco_cosTheta[index],1.0);
248
249 // E_char
250 if (E_char < echar_max){
251
252 // h7_E->Fill(CalibratedReco_energy[index],1.0); // events after Echar
253 // h7_cos->Fill(CalibratedReco_cosTheta[index],1.0);
254
255
256 // Bhabha rejection in forward calorimeters
257 h8_E->Fill(CalibratedReco_energy[index],EventWeight);
258 h8_cos->Fill(CalibratedReco_cosTheta[index],EventWeight);
259
260 n_after_all_cuts += EventWeight;
261 pass = true;
262
263 }
264 }
265 }
266 }
267
268 // ################# FILL SELECTION TREE VARIABLES WITH SELECTED PHOTON
269
270 if (pass){
271 nSelPhotons = 0;
272 SelPhoton_id[nSelPhotons] = 1;
273 SelPhoton_type[nSelPhotons] = CalibratedReco_type[index];
274 SelPhoton_px[nSelPhotons] = CalibratedReco_px[index];
275 SelPhoton_py[nSelPhotons] = CalibratedReco_py[index];
276 SelPhoton_pz[nSelPhotons] = CalibratedReco_pz[index];
277 SelPhoton_energy[nSelPhotons] = CalibratedReco_energy[index];
278 SelPhoton_mass[nSelPhotons] = CalibratedReco_mass[index];
279 SelPhoton_charge[nSelPhotons] = CalibratedReco_charge[index];
280 SelPhoton_position_x[nSelPhotons] = CalibratedReco_position_x[index];
281 SelPhoton_position_y[nSelPhotons] = CalibratedReco_position_y[index];
282 SelPhoton_position_z[nSelPhotons] = CalibratedReco_position_z[index];
283
284 SelPhoton_absP[nSelPhotons] = CalibratedReco_absP[index];
285 SelPhoton_cosTheta[nSelPhotons] = CalibratedReco_cosTheta[index];
286 SelPhoton_phi[nSelPhotons] = CalibratedReco_phi[index];
287
288 SelPhotonToCalibratedId[nSelPhotons] = CalibratedReco_id[index];
289 SelPhotonRelatedToId[nSelPhotons] = CalibratedRecoRelatedToId[index];
290
291 nSelPhotons++;
292
293 sel_tree->Fill();
294 sel_calibratedtree->Fill();
295 sel_gammatree->Fill();
296 // if (!strcmp(Process,"n1n1a") || !strcmp(Process,"n2n2a") || !strcmp(Process,"n3n3a")){
297 // sel_xsectree->Fill();
298 // }
299 sel_bcaltree->Fill();
300 }
301 }
302
303 // ############### END OF EVENT LOOP
304
305 printf("%i empty events in file: %s\n",n_empty,name);
306
307 sel->cd();
308
309 Entries_BeforeCuts = n_reco;
310 Entries_AfterKinCuts = n_after_kinetic_cuts;
311 Entries_AfterAllCuts = n_after_all_cuts;
312 CrossSection_AfterKinCuts = CrossSection*((double)Entries_AfterKinCuts/Entries_BeforeCuts);
313
314 Efficiency_KinCuts = Entries_AfterKinCuts/Entries_BeforeCuts;
315 Efficiency_RecCuts = Entries_AfterAllCuts/Entries_AfterKinCuts;
316 Efficiency_AllCuts = Entries_AfterAllCuts/Entries_BeforeCuts;
317
318 TTree* sel_efficiencies = new TTree("SelectionEfficiencies","SelectionEfficiencies");
319 SetSelectionEfficienciesTreeBranches(sel_efficiencies);
320
321 sel_efficiencies->Fill();
322
323 // ######### WRITE TREES TO SELECTION FILE
324
325
326 sel_tree->Write();
327 sel_globaltree->Write();
328 sel_calibratedtree->Write();
329 sel_gammatree->Write();
330 // if (!strcmp(Process,"n1n1a") || !strcmp(Process,"n2n2a") || !strcmp(Process,"n3n3a")){
331 // sel_xsectree->Write();
332 // }
333 sel_bcaltree->Write();
334 sel_efficiencies->Write();
335
336 f->Close();
337 sel->Close();
338 delete f;
339
340 f_cutflow->cd();
341
342 h1_E->Write();
343 h1_cos->Write();
344 h2_E->Write();
345 h2_cos->Write();
346 h5_E->Write();
347 h5_cos->Write();
348 h6_E->Write();
349 h6_cos->Write();
350 // h7_E->Write();
351 // h7_cos->Write();
352 h8_E->Write();
353 h8_cos->Write();
354 }
355 f_cutflow->Close();
356 // // ################# END OF FILE LOOP
357 // CalculateEfficiencies();
358
359 gROOT->ProcessLine(".q");
360 return;
361 }
362
363
364 // #############################################################
365 // #############################################################
366
367
368 void CalculateEfficiencies(){
369
370 double n_before = 0;
371 double n_afterkin = 0;
372 double n_afterall = 0;
373 //double avXSec = .0;
374
375 double N_after_norm[7] = {0.0};
376 double N_before_norm[7] = {0.0};
377 double polarisedXSec[7] = {0.0}; // 0/0, 80/00, 80/30, 80/60, -80/00, -80/30, -80/60
378 double polarisedXSec_recoCut[7] = {0.0};
379
380 printf("Calculate polarisation dependent selection efficiencies for signal-like processes\n");
381
382 // ###############################
383 TFile* f = 0;
384 ifstream txtfile("../steeringfiles/ILD00_Selection_RootFiles.txt");
385 ifstream txtfile2("../steeringfiles/ILD00_Selection_RootFiles.txt");
386 char name[256];
387 txtfile2.getline(name,256,'\n');
388 string file;
389 string nextfile;
390
391 while (!txtfile.eof()){
392 txtfile.getline(name,256,'\n');
393 file = name;
394
395 if (txtfile2.eof()) nextfile = "eof";
396 else{
397 txtfile2.getline(name,256,'\n');
398 nextfile = name;
399 }
400 printf("File: %s\nnext file: %s\n\n",file.c_str(),nextfile.c_str());
401
402 // ############################# OPEN ORIGINAL FILE, LOAD TREES
403
404 f = new TFile(file.c_str(),"OPEN");
405 if (f->IsZombie()) continue;
406
407 TTree* eff_tree = (TTree*)f->Get("SelectionEfficiencies");
408 TTree* tree = (TTree*)f->Get("Tree");
409 TTree* globaltree = (TTree*)f->Get("GlobalParameters");
410 LoadTree(tree);
411 LoadGlobalTree(globaltree);
412 LoadSelectionEfficienciesTree(eff_tree);
413
414
415 eff_tree->GetEntry(0);
416 tree->GetEntry(0);
417 globaltree->GetEntry(0);
418
419 n_before += Entries_BeforeCuts; // sum up for each process and polarisation (-1/1) (1/-1)
420 n_afterkin += Entries_AfterKinCuts;
421 n_afterall += Entries_AfterAllCuts;
422
423 if ((file.substr(0,52)).compare(nextfile.substr(0,52))){ // if next file is different process
424 // for single gamma processes
425 if (!strcmp(Process,"n1n1a") || !strcmp(Process,"n2n2a") || !strcmp(Process,"n3n3a")){
426 CreatePolarisationWeights(0.0,0.0);
427 // Lumi normalised events after selection cuts
428 N_after_norm[0] += n_afterall/MCLumi* GetPolarisationWeight(P_e,P_p);
429 // Lumi normalised events before selection cuts
430 N_before_norm[0] += n_afterkin/MCLumi* GetPolarisationWeight(P_e,P_p);
431
432 CreatePolarisationWeights(0.8,0.0);
433 // Lumi normalised events after selection cuts
434 N_after_norm[1] += n_afterall/MCLumi* GetPolarisationWeight(P_e,P_p);
435 // Lumi normalised events before selection cuts
436 N_before_norm[1] += n_afterkin/MCLumi* GetPolarisationWeight(P_e,P_p);
437
438
439 CreatePolarisationWeights(0.8,-0.3);
440 // Lumi normalised events after selection cuts
441 N_after_norm[2] += n_afterall/MCLumi* GetPolarisationWeight(P_e,P_p);
442 // Lumi normalised events before selection cuts
443 N_before_norm[2] += n_afterkin/MCLumi* GetPolarisationWeight(P_e,P_p);
444
445 CreatePolarisationWeights(0.8,-0.6);
446 // Lumi normalised events after selection cuts
447 N_after_norm[3] += n_afterall/MCLumi* GetPolarisationWeight(P_e,P_p);
448 // Lumi normalised events before selection cuts
449 N_before_norm[3] += n_afterkin/MCLumi* GetPolarisationWeight(P_e,P_p);
450
451 CreatePolarisationWeights(-0.8,0.0);
452 // Lumi normalised events after selection cuts
453 N_after_norm[4] += n_afterall/MCLumi* GetPolarisationWeight(P_e,P_p);
454 // Lumi normalised events before selection cuts
455 N_before_norm[4] += n_afterkin/MCLumi* GetPolarisationWeight(P_e,P_p);
456
457
458 CreatePolarisationWeights(-0.8,0.3);
459 // Lumi normalised events after selection cuts
460 N_after_norm[5] += n_afterall/MCLumi* GetPolarisationWeight(P_e,P_p);
461 // Lumi normalised events before selection cuts
462 N_before_norm[5] += n_afterkin/MCLumi* GetPolarisationWeight(P_e,P_p);
463
464
465 CreatePolarisationWeights(-0.8,0.6);
466 // Lumi normalised events after selection cuts
467 N_after_norm[6] += n_afterall/MCLumi* GetPolarisationWeight(P_e,P_p);
468 // Lumi normalised events before selection cuts
469 N_before_norm[6] += n_afterkin/MCLumi* GetPolarisationWeight(P_e,P_p);
470 }
471 n_before = 0;
472 n_afterkin = 0;
473 n_afterall = 0;
474
475 }
476 }
477
478
479 TFile* f_out = new TFile("../data/Auxiliary/SelectionScalingFactors.root","RECREATE");
480 TTree* t = new TTree("Scalingfactors","Scalingfactors");
481 SetScalingFactorsTreeBranches(t);
482
483
484 ScalingFactorp0p0 = N_after_norm[0]/N_before_norm[0];
485 ScalingFactorp8p0 = N_after_norm[1]/N_before_norm[1];
486 ScalingFactorp8m3 = N_after_norm[2]/N_before_norm[2];
487 ScalingFactorp8m6 = N_after_norm[3]/N_before_norm[3];
488 ScalingFactorm8p0 = N_after_norm[4]/N_before_norm[4];
489 ScalingFactorm8p3 = N_after_norm[5]/N_before_norm[5];
490 ScalingFactorm8p6 = N_after_norm[6]/N_before_norm[6];
491
492
493 PolXSecp0p0 = polarisedXSec_recoCut[0];
494 PolXSecp8p0 = polarisedXSec_recoCut[1];
495 PolXSecp8m3 = polarisedXSec_recoCut[2];
496 PolXSecp8m6 = polarisedXSec_recoCut[3];
497 PolXSecm8p0 = polarisedXSec_recoCut[4];
498 PolXSecm8p3 = polarisedXSec_recoCut[5];
499 PolXSecm8p6 = polarisedXSec_recoCut[6];
500
501 // RECO CUT SCALING FACTORS FOR DIFFERENT POLARISATIONS
502 printf("Cross section scaling factors\n");
503 printf("(0.0/ 0.0): %.4f\n",ScalingFactorp0p0);
504 printf("(0.8/ 0.0): %.4f\n",ScalingFactorp8p0);
505 printf("(0.8/-0.3): %.4f\n",ScalingFactorp8m3);
506 printf("(0.8/-0.6): %.4f\n",ScalingFactorp8m6);
507 printf("(-0.8/0.0): %.4f\n",ScalingFactorm8p0);
508 printf("(-0.8/0.3): %.4f\n",ScalingFactorm8p3);
509 printf("(-0.8/0.6): %.4f\n",ScalingFactorm8p6);
510
511 //
512 // printf("\nCross sections\n");
513 // printf("(0.0/ 0.0): %.4f\n",PolXSecp0p0);
514 // printf("(0.8/ 0.0): %.4f\n",PolXSecp8p0);
515 // printf("(0.8/-0.3): %.4f\n",PolXSecp8m3);
516 // printf("(0.8/-0.6): %.4f\n",PolXSecp8m6);
517 // printf("(-0.8/0.3): %.4f\n",PolXSecm8p3);
518
519 t->Fill();
520 f_out->cd();
521 t->Write();
522
523 f_out->Close();
524 return;
525 }
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.