Fork me on GitHub

source: svn/trunk/src/SmearUtil.cc@ 323

Last change on this file since 323 was 321, checked in by Xavier Rouby, 16 years ago

Isolation updated. ptiso implemented. etrat prepared but not finished
File size: 63.7 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \----------------------------------------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
12** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
13** FROG: [hep-ex/0901.2718v1] **
14** **
15** ------------------------------------------------------------------ **
16** **
17** Main authors: **
18** ------------- **
19** **
20** Severine Ovyn Xavier Rouby **
21** severine.ovyn@uclouvain.be xavier.rouby@cern **
22** **
23** Center for Particle Physics and Phenomenology (CP3) **
24** Universite catholique de Louvain (UCL) **
25** Louvain-la-Neuve, Belgium **
26** **
27** Copyright (C) 2008-2009, **
28** All rights reserved. **
29** **
30***********************************************************************/
31
32
33/// \file SmearUtil.cc
34/// \brief RESOLution class, and some generic definitions
35
36
37#include "SmearUtil.h"
38#include "TRandom.h"
39
40#include <iostream>
41#include <fstream>
42#include <sstream>
43#include <iomanip>
44using namespace std;
45
46//------------------------------------------------------------------------------
47
48RESOLution::RESOLution() {
49
50 // Detector characteristics
51 CEN_max_tracker = 2.5; // Maximum tracker coverage
52 CEN_max_calo_cen = 3.0; // central calorimeter coverage
53 CEN_max_calo_fwd = 5.0; // forward calorimeter pseudorapidity coverage
54 CEN_max_mu = 2.4; // muon chambers pseudorapidity coverage
55
56 // Energy resolution for electron/photon
57 // \sigma/E = C + N/E + S/\sqrt{E}
58 ELG_Scen = 0.05; // S term for central ECAL
59 ELG_Ncen = 0.25; // N term for central ECAL
60 ELG_Ccen = 0.005; // C term for central ECAL
61 ELG_Sfwd = 2.084; // S term for FCAL
62 ELG_Nfwd = 0.0; // N term for FCAL
63 ELG_Cfwd = 0.107; // C term for FCAL
64
65 // Energy resolution for hadrons in ecal/hcal/hf
66 // \sigma/E = C + N/E + S/\sqrt{E}
67 HAD_Shcal = 1.5; // S term for central HCAL
68 HAD_Nhcal = 0.; // N term for central HCAL
69 HAD_Chcal = 0.05; // C term for central HCAL
70 HAD_Shf = 2.7; // S term for FCAL
71 HAD_Nhf = 0.; // N term for FCAL
72 HAD_Chf = 0.13; // C term for FCAL
73
74 // Muon smearing
75 MU_SmearPt = 0.01;
76
77 // Tracking efficiencies
78 TRACK_ptmin = 0.9; // minimal pt needed to reach the calorimeter in GeV
79 TRACK_eff = 100; // efficiency associated to the tracking
80
81 // Calorimetric towers
82 TOWER_number = 40;
83 const float tower_eta_edges[41] = {
84 0., 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566,
85 1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.868, 2.950, 3.125, 3.300, 3.475, 3.650, 3.825, 4.000, 4.175,
86 4.350, 4.525, 4.700, 5.000}; // temporary object
87 TOWER_eta_edges = new float[TOWER_number+1];
88 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = tower_eta_edges[i];
89
90 const float tower_dphi[40] = {
91 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 10,
92 10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; // temporary object
93 TOWER_dphi = new float[TOWER_number];
94 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = tower_dphi[i];
95
96
97 // Thresholds for reconstructed objetcs
98 PTCUT_elec = 10.0;
99 PTCUT_muon = 10.0;
100 PTCUT_jet = 20.0;
101 PTCUT_gamma = 10.0;
102 PTCUT_taujet = 10.0;
103
104 // Isolation
105 ISOL_PT = 2.0; //minimal pt of tracks for isolation criteria
106 ISOL_Cone = 0.5; //Cone for isolation criteria
107 ISOL_Calo_ET = 1E99; //minimal tower energy for isolation criteria. Default off = 1E99
108 ISOL_Calo_Cone = 0.5; //Cone for calorimetric isolation
109 ISOL_Calo_Grid = 3; //Grid size (N x N) for calorimetric isolation
110
111 // General jet variable
112 JET_coneradius = 0.7; // generic jet radius ; not for tau's !!!
113 JET_jetalgo = 1; // 1 for Cone algorithm, 2 for MidPoint algorithm, 3 for SIScone algorithm, 4 for kt algorithm
114 JET_seed = 1.0; // minimum seed to start jet reconstruction
115
116 // Tagging definition
117 BTAG_b = 40;
118 BTAG_mistag_c = 10;
119 BTAG_mistag_l = 1;
120
121 // FLAGS
122 FLAG_bfield = 1; //1 to run the bfield propagation else 0
123 FLAG_vfd = 1; //1 to run the very forward detectors else 0
124 FLAG_RP = 1; //1 to run the zero degree calorimeter else 0
125 FLAG_trigger = 1; //1 to run the trigger selection else 0
126 FLAG_frog = 1; //1 to run the FROG event display
127 FLAG_lhco = 1;
128
129 // In case BField propagation allowed
130 TRACK_radius = 129; //radius of the BField coverage
131 TRACK_length = 300; //length of the BField coverage
132 TRACK_bfield_x = 0; //X composant of the BField
133 TRACK_bfield_y = 0; //Y composant of the BField
134 TRACK_bfield_z = 3.8; //Z composant of the BField
135
136 // In case Very forward detectors allowed
137 VFD_min_calo_vfd = 5.2; // very forward calorimeter (if any) like CASTOR
138 VFD_max_calo_vfd = 6.6;
139 VFD_min_zdc = 8.3;
140 VFD_s_zdc = 140; // distance of the Zero Degree Calorimeter, from the Interaction poin, in [m]
141
142 RP_220_s = 220; // distance of the RP to the IP, in meters
143 RP_220_x = 0.002; // distance of the RP to the beam, in meters
144 RP_420_s = 420; // distance of the RP to the IP, in meters
145 RP_420_x = 0.004; // distance of the RP to the beam, in meters
146 RP_IP_name = "IP5";
147 RP_beam1Card = "data/LHCB1IR5_v6.500.tfs";
148 RP_beam2Card = "data/LHCB1IR5_v6.500.tfs";
149
150 // In case FROG event display allowed
151 NEvents_Frog = 10;
152
153 //********************************************
154 //jet stuffs not defined in the input datacard
155 //********************************************
156
157 JET_overlap = 0.75;
158 // MidPoint algorithm definition
159 JET_M_coneareafraction = 0.25;
160 JET_M_maxpairsize = 2;
161 JET_M_maxiterations = 100;
162 // Define Cone algorithm.
163 JET_C_adjacencycut = 2;
164 JET_C_maxiterations = 100;
165 JET_C_iratch = 1;
166 //Define SISCone algorithm.
167 JET_S_npass = 0;
168 JET_S_protojet_ptmin= 0.0;
169
170 //For Tau-jet definition
171 TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
172 TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
173 TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
174 TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
175
176 PT_QUARKS_MIN = 2.0 ; // minimal pt needed by quarks to do b-tag
177
178 //for very forward detectors
179 RP_offsetEl_s = 120;
180 RP_offsetEl_x = 0.097;
181 RP_cross_x = -500;
182 RP_cross_y = 0.0;
183 RP_cross_ang = 142.5;
184
185}
186
187
188RESOLution::RESOLution(const RESOLution & DET) {
189 // Detector characteristics
190 CEN_max_tracker = DET.CEN_max_tracker;
191 CEN_max_calo_cen = DET.CEN_max_calo_cen;
192 CEN_max_calo_fwd = DET.CEN_max_calo_fwd;
193 CEN_max_mu = DET.CEN_max_mu;
194
195 // Energy resolution for electron/photon
196 ELG_Scen = DET.ELG_Scen;
197 ELG_Ncen = DET.ELG_Ncen;
198 ELG_Ccen = DET.ELG_Ccen;
199 ELG_Cfwd = DET.ELG_Cfwd;
200 ELG_Sfwd = DET.ELG_Sfwd;
201 ELG_Nfwd = DET.ELG_Nfwd;
202
203 // Energy resolution for hadrons in ecal/hcal/hf
204 HAD_Shcal = DET.HAD_Shcal;
205 HAD_Nhcal = DET.HAD_Nhcal;
206 HAD_Chcal = DET.HAD_Chcal;
207 HAD_Shf = DET.HAD_Shf;
208 HAD_Nhf = DET.HAD_Nhf;
209 HAD_Chf = DET.HAD_Chf;
210
211 // Muon smearing
212 MU_SmearPt = DET.MU_SmearPt;
213
214 // Tracking efficiencies
215 TRACK_ptmin = DET.TRACK_ptmin;
216 TRACK_eff = DET.TRACK_eff;
217
218 // Calorimetric towers
219 TOWER_number = DET.TOWER_number;
220 TOWER_eta_edges = new float[TOWER_number+1];
221 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];
222
223 TOWER_dphi = new float[TOWER_number];
224 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];
225
226 // Thresholds for reconstructed objetcs
227 PTCUT_elec = DET.PTCUT_elec;
228 PTCUT_muon = DET.PTCUT_muon;
229 PTCUT_jet = DET.PTCUT_jet;
230 PTCUT_gamma = DET.PTCUT_gamma;
231 PTCUT_taujet = DET.PTCUT_taujet;
232
233 // Isolation
234 ISOL_PT = DET.ISOL_PT; // tracking isolation
235 ISOL_Cone = DET.ISOL_Cone;
236 ISOL_Calo_ET = DET.ISOL_Calo_ET; // calorimeter isolation, defaut off
237 ISOL_Calo_Cone = DET.ISOL_Calo_Cone;
238 ISOL_Calo_Grid = DET.ISOL_Calo_Grid;
239
240
241 // General jet variable
242 JET_coneradius = DET.JET_coneradius;
243 JET_jetalgo = DET.JET_jetalgo;
244 JET_seed = DET.JET_seed;
245
246 // Tagging definition
247 BTAG_b = DET.BTAG_b;
248 BTAG_mistag_c = DET.BTAG_mistag_c;
249 BTAG_mistag_l = DET.BTAG_mistag_l;
250
251 // FLAGS
252 FLAG_bfield = DET.FLAG_bfield;
253 FLAG_vfd = DET.FLAG_vfd;
254 FLAG_RP = DET.FLAG_RP;
255 FLAG_trigger = DET.FLAG_trigger;
256 FLAG_frog = DET.FLAG_frog;
257 FLAG_lhco = DET.FLAG_lhco;
258
259 // In case BField propagation allowed
260 TRACK_radius = DET.TRACK_radius;
261 TRACK_length = DET.TRACK_length;
262 TRACK_bfield_x = DET.TRACK_bfield_x;
263 TRACK_bfield_y = DET.TRACK_bfield_y;
264 TRACK_bfield_z = DET.TRACK_bfield_z;
265
266 // In case Very forward detectors allowed
267 VFD_min_calo_vfd = DET.VFD_min_calo_vfd;
268 VFD_max_calo_vfd = DET.VFD_max_calo_vfd;
269 VFD_min_zdc = DET.VFD_min_zdc;
270 VFD_s_zdc = DET.VFD_s_zdc;
271
272 RP_220_s = DET.RP_220_s;
273 RP_220_x = DET.RP_220_x;
274 RP_420_s = DET.RP_420_s;
275 RP_420_x = DET.RP_420_x;
276 RP_beam1Card = DET.RP_beam1Card;
277 RP_beam2Card = DET.RP_beam2Card;
278 RP_offsetEl_s = DET.RP_offsetEl_s;
279 RP_offsetEl_x = DET.RP_offsetEl_x;
280 RP_cross_x = DET.RP_cross_x;
281 RP_cross_y = DET.RP_cross_y;
282 RP_cross_ang = DET.RP_cross_ang;
283 RP_IP_name = DET.RP_IP_name;
284
285 // In case FROG event display allowed
286 NEvents_Frog = DET.NEvents_Frog;
287
288 JET_overlap = DET.JET_overlap;
289 // MidPoint algorithm definition
290 JET_M_coneareafraction = DET.JET_M_coneareafraction;
291 JET_M_maxpairsize = DET.JET_M_maxpairsize;
292 JET_M_maxiterations = DET.JET_M_maxiterations;
293 // Define Cone algorithm.
294 JET_C_adjacencycut = DET.JET_C_adjacencycut;
295 JET_C_maxiterations = DET.JET_C_maxiterations;
296 JET_C_iratch = DET.JET_C_iratch;
297 //Define SISCone algorithm.
298 JET_S_npass = DET.JET_S_npass;
299 JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin;
300
301 //For Tau-jet definition
302 TAU_energy_scone = DET.TAU_energy_scone;
303 TAU_track_scone = DET.TAU_track_scone;
304 TAU_track_pt = DET.TAU_track_pt;
305 TAU_energy_frac = DET.TAU_energy_frac;
306
307 PT_QUARKS_MIN = DET.PT_QUARKS_MIN;
308}
309
310RESOLution& RESOLution::operator=(const RESOLution& DET) {
311 if(this==&DET) return *this;
312 // Detector characteristics
313 CEN_max_tracker = DET.CEN_max_tracker;
314 CEN_max_calo_cen = DET.CEN_max_calo_cen;
315 CEN_max_calo_fwd = DET.CEN_max_calo_fwd;
316 CEN_max_mu = DET.CEN_max_mu;
317
318 // Energy resolution for electron/photon
319 ELG_Scen = DET.ELG_Scen;
320 ELG_Ncen = DET.ELG_Ncen;
321 ELG_Ccen = DET.ELG_Ccen;
322 ELG_Cfwd = DET.ELG_Cfwd;
323 ELG_Sfwd = DET.ELG_Sfwd;
324 ELG_Nfwd = DET.ELG_Nfwd;
325
326 // Energy resolution for hadrons in ecal/hcal/hf
327 HAD_Shcal = DET.HAD_Shcal;
328 HAD_Nhcal = DET.HAD_Nhcal;
329 HAD_Chcal = DET.HAD_Chcal;
330 HAD_Shf = DET.HAD_Shf;
331 HAD_Nhf = DET.HAD_Nhf;
332 HAD_Chf = DET.HAD_Chf;
333
334 // Muon smearing
335 MU_SmearPt = DET.MU_SmearPt;
336
337 // Tracking efficiencies
338 TRACK_ptmin = DET.TRACK_ptmin;
339 TRACK_eff = DET.TRACK_eff;
340
341 // Calorimetric towers
342 TOWER_number = DET.TOWER_number;
343 TOWER_eta_edges = new float[TOWER_number+1];
344 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];
345
346 TOWER_dphi = new float[TOWER_number];
347 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];
348
349 // Thresholds for reconstructed objetcs
350 PTCUT_elec = DET.PTCUT_elec;
351 PTCUT_muon = DET.PTCUT_muon;
352 PTCUT_jet = DET.PTCUT_jet;
353 PTCUT_gamma = DET.PTCUT_gamma;
354 PTCUT_taujet = DET.PTCUT_taujet;
355
356 // Isolation
357 ISOL_PT = DET.ISOL_PT; // tracking isolation
358 ISOL_Cone = DET.ISOL_Cone;
359 ISOL_Calo_ET = DET.ISOL_Calo_ET; // calorimeter isolation, defaut off
360 ISOL_Calo_Cone = DET.ISOL_Calo_Cone;
361 ISOL_Calo_Grid = DET.ISOL_Calo_Grid;
362
363 // General jet variable
364 JET_coneradius = DET.JET_coneradius;
365 JET_jetalgo = DET.JET_jetalgo;
366 JET_seed = DET.JET_seed;
367
368 // Tagging definition
369 BTAG_b = DET.BTAG_b;
370 BTAG_mistag_c = DET.BTAG_mistag_c;
371 BTAG_mistag_l = DET.BTAG_mistag_l;
372
373 // FLAGS
374 FLAG_bfield = DET.FLAG_bfield;
375 FLAG_vfd = DET.FLAG_vfd;
376 FLAG_RP = DET.FLAG_RP;
377 FLAG_trigger = DET.FLAG_trigger;
378 FLAG_frog = DET.FLAG_frog;
379 FLAG_lhco = DET.FLAG_lhco;
380
381 // In case BField propagation allowed
382 TRACK_radius = DET.TRACK_radius;
383 TRACK_length = DET.TRACK_length;
384 TRACK_bfield_x = DET.TRACK_bfield_x;
385 TRACK_bfield_y = DET.TRACK_bfield_y;
386 TRACK_bfield_z = DET.TRACK_bfield_z;
387
388 // In case Very forward detectors allowed
389 VFD_min_calo_vfd = DET.VFD_min_calo_vfd;
390 VFD_max_calo_vfd = DET.VFD_max_calo_vfd;
391 VFD_min_zdc = DET.VFD_min_zdc;
392 VFD_s_zdc = DET.VFD_s_zdc;
393
394 RP_220_s = DET.RP_220_s;
395 RP_220_x = DET.RP_220_x;
396 RP_420_s = DET.RP_420_s;
397 RP_420_x = DET.RP_420_x;
398 RP_offsetEl_s = DET.RP_offsetEl_s;
399 RP_offsetEl_x = DET.RP_offsetEl_x;
400 RP_beam1Card = DET.RP_beam1Card;
401 RP_beam2Card = DET.RP_beam2Card;
402 RP_cross_x = DET.RP_cross_x;
403 RP_cross_y = DET.RP_cross_y;
404 RP_cross_ang = DET.RP_cross_ang;
405 RP_IP_name = DET.RP_IP_name;
406
407
408 // In case FROG event display allowed
409 NEvents_Frog = DET.NEvents_Frog;
410
411 JET_overlap = DET.JET_overlap;
412 // MidPoint algorithm definition
413 JET_M_coneareafraction = DET.JET_M_coneareafraction;
414 JET_M_maxpairsize = DET.JET_M_maxpairsize;
415 JET_M_maxiterations = DET.JET_M_maxiterations;
416 // Define Cone algorithm.
417 JET_C_adjacencycut = DET.JET_C_adjacencycut;
418 JET_C_maxiterations = DET.JET_C_maxiterations;
419 JET_C_iratch = DET.JET_C_iratch;
420 //Define SISCone algorithm.
421 JET_S_npass = DET.JET_S_npass;
422 JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin;
423
424 //For Tau-jet definition
425 TAU_energy_scone = DET.TAU_energy_scone;
426 TAU_track_scone = DET.TAU_track_scone;
427 TAU_track_pt = DET.TAU_track_pt;
428 TAU_energy_frac = DET.TAU_energy_frac;
429
430 PT_QUARKS_MIN = DET.PT_QUARKS_MIN;
431 return *this;
432}
433
434
435
436
437//------------------------------------------------------------------------------
438void RESOLution::ReadDataCard(const string datacard) {
439
440 string temp_string;
441 istringstream curstring;
442
443 ifstream fichier_a_lire(datacard.c_str());
444 if(!fichier_a_lire.good()) {
445 cout <<"** WARNING: Datadard not found, use default values **" << endl;
446 return;
447 }
448
449 while (getline(fichier_a_lire,temp_string)) {
450 curstring.clear(); // needed when using several times istringstream::str(string)
451 curstring.str(temp_string);
452 string varname;
453 float value; int ivalue; string svalue;
454
455 if(strstr(temp_string.c_str(),"#")) { }
456 else if(strstr(temp_string.c_str(),"CEN_max_tracker")) {curstring >> varname >> value; CEN_max_tracker = value;}
457 else if(strstr(temp_string.c_str(),"CEN_max_calo_cen")) {curstring >> varname >> value; CEN_max_calo_cen = value;}
458 else if(strstr(temp_string.c_str(),"CEN_max_calo_fwd")) {curstring >> varname >> value; CEN_max_calo_fwd = value;}
459 else if(strstr(temp_string.c_str(),"CEN_max_mu")) {curstring >> varname >> value; CEN_max_mu = value;}
460
461 else if(strstr(temp_string.c_str(),"VFD_min_calo_vfd")) {curstring >> varname >> value; VFD_min_calo_vfd = value;}
462 else if(strstr(temp_string.c_str(),"VFD_max_calo_vfd")) {curstring >> varname >> value; VFD_max_calo_vfd = value;}
463 else if(strstr(temp_string.c_str(),"VFD_min_zdc")) {curstring >> varname >> value; VFD_min_zdc = value;}
464 else if(strstr(temp_string.c_str(),"VFD_s_zdc")) {curstring >> varname >> value; VFD_s_zdc = value;}
465
466 else if(strstr(temp_string.c_str(),"RP_220_s")) {curstring >> varname >> value; RP_220_s = value;}
467 else if(strstr(temp_string.c_str(),"RP_220_x")) {curstring >> varname >> value; RP_220_x = value;}
468 else if(strstr(temp_string.c_str(),"RP_420_s")) {curstring >> varname >> value; RP_420_s = value;}
469 else if(strstr(temp_string.c_str(),"RP_420_x")) {curstring >> varname >> value; RP_420_x = value;}
470 else if(strstr(temp_string.c_str(),"RP_beam1Card")) {curstring >> varname >> svalue;RP_beam1Card = svalue;}
471 else if(strstr(temp_string.c_str(),"RP_beam2Card")) {curstring >> varname >> svalue;RP_beam2Card = svalue;}
472 else if(strstr(temp_string.c_str(),"RP_IP_name")) {curstring >> varname >> svalue;RP_IP_name = svalue;}
473
474 else if(strstr(temp_string.c_str(),"ELG_Scen")) {curstring >> varname >> value; ELG_Scen = value;}
475 else if(strstr(temp_string.c_str(),"ELG_Ncen")) {curstring >> varname >> value; ELG_Ncen = value;}
476 else if(strstr(temp_string.c_str(),"ELG_Ccen")) {curstring >> varname >> value; ELG_Ccen = value;}
477 else if(strstr(temp_string.c_str(),"ELG_Sfwd")) {curstring >> varname >> value; ELG_Sfwd = value;}
478 else if(strstr(temp_string.c_str(),"ELG_Cfwd")) {curstring >> varname >> value; ELG_Cfwd = value;}
479 else if(strstr(temp_string.c_str(),"ELG_Nfwd")) {curstring >> varname >> value; ELG_Nfwd = value;}
480 else if(strstr(temp_string.c_str(),"HAD_Shcal")) {curstring >> varname >> value; HAD_Shcal = value;}
481 else if(strstr(temp_string.c_str(),"HAD_Nhcal")) {curstring >> varname >> value; HAD_Nhcal = value;}
482 else if(strstr(temp_string.c_str(),"HAD_Chcal")) {curstring >> varname >> value; HAD_Chcal = value;}
483 else if(strstr(temp_string.c_str(),"HAD_Shf")) {curstring >> varname >> value; HAD_Shf = value;}
484 else if(strstr(temp_string.c_str(),"HAD_Nhf")) {curstring >> varname >> value; HAD_Nhf = value;}
485 else if(strstr(temp_string.c_str(),"HAD_Chf")) {curstring >> varname >> value; HAD_Chf = value;}
486 else if(strstr(temp_string.c_str(),"MU_SmearPt")) {curstring >> varname >> value; MU_SmearPt = value;}
487
488 else if(strstr(temp_string.c_str(),"TRACK_radius")) {curstring >> varname >> ivalue;TRACK_radius = ivalue;}
489 else if(strstr(temp_string.c_str(),"TRACK_length")) {curstring >> varname >> ivalue;TRACK_length = ivalue;}
490 else if(strstr(temp_string.c_str(),"TRACK_bfield_x")) {curstring >> varname >> value; TRACK_bfield_x = value;}
491 else if(strstr(temp_string.c_str(),"TRACK_bfield_y")) {curstring >> varname >> value; TRACK_bfield_y = value;}
492 else if(strstr(temp_string.c_str(),"TRACK_bfield_z")) {curstring >> varname >> value; TRACK_bfield_z = value;}
493 else if(strstr(temp_string.c_str(),"FLAG_bfield")) {curstring >> varname >> ivalue; FLAG_bfield = ivalue;}
494 else if(strstr(temp_string.c_str(),"TRACK_ptmin")) {curstring >> varname >> value; TRACK_ptmin = value;}
495 else if(strstr(temp_string.c_str(),"TRACK_eff")) {curstring >> varname >> ivalue;TRACK_eff = ivalue;}
496
497 else if(strstr(temp_string.c_str(),"TOWER_number")) {curstring >> varname >> ivalue;TOWER_number = ivalue;}
498 else if(strstr(temp_string.c_str(),"TOWER_eta_edges")){
499 curstring >> varname; for(unsigned int i=0; i<TOWER_number+1; i++) {curstring >> value; TOWER_eta_edges[i] = value;} }
500 else if(strstr(temp_string.c_str(),"TOWER_dphi")){
501 curstring >> varname; for(unsigned int i=0; i<TOWER_number; i++) {curstring >> value; TOWER_dphi[i] = value;} }
502
503 else if(strstr(temp_string.c_str(),"PTCUT_elec")) {curstring >> varname >> value; PTCUT_elec = value;}
504 else if(strstr(temp_string.c_str(),"PTCUT_muon")) {curstring >> varname >> value; PTCUT_muon = value;}
505 else if(strstr(temp_string.c_str(),"PTCUT_jet")) {curstring >> varname >> value; PTCUT_jet = value;}
506 else if(strstr(temp_string.c_str(),"PTCUT_gamma")) {curstring >> varname >> value; PTCUT_gamma = value;}
507 else if(strstr(temp_string.c_str(),"PTCUT_taujet")) {curstring >> varname >> value; PTCUT_taujet = value;}
508
509 else if(strstr(temp_string.c_str(),"ISOL_PT")) {curstring >> varname >> value; ISOL_PT = value;}
510 else if(strstr(temp_string.c_str(),"ISOL_Cone")) {curstring >> varname >> value; ISOL_Cone = value;}
511 else if(strstr(temp_string.c_str(),"ISOL_Calo_ET")) {curstring >> varname >> value; ISOL_Calo_ET = value;}
512 else if(strstr(temp_string.c_str(),"ISOL_Calo_Cone")) {curstring >> varname >> value; ISOL_Calo_Cone = value;}
513 else if(strstr(temp_string.c_str(),"ISOL_Calo_Grid")) {curstring >> varname >> ivalue; ISOL_Calo_Grid = ivalue;}
514
515 else if(strstr(temp_string.c_str(),"JET_coneradius")) {curstring >> varname >> value; JET_coneradius = value;}
516 else if(strstr(temp_string.c_str(),"JET_jetalgo")) {curstring >> varname >> ivalue;JET_jetalgo = ivalue;}
517 else if(strstr(temp_string.c_str(),"JET_seed")) {curstring >> varname >> value; JET_seed = value;}
518
519 else if(strstr(temp_string.c_str(),"BTAG_b")) {curstring >> varname >> ivalue;BTAG_b = ivalue;}
520 else if(strstr(temp_string.c_str(),"BTAG_mistag_c")) {curstring >> varname >> ivalue;BTAG_mistag_c = ivalue;}
521 else if(strstr(temp_string.c_str(),"BTAG_mistag_l")) {curstring >> varname >> ivalue;BTAG_mistag_l = ivalue;}
522
523 else if(strstr(temp_string.c_str(),"FLAG_vfd")) {curstring >> varname >> ivalue; FLAG_vfd = ivalue;}
524 else if(strstr(temp_string.c_str(),"FLAG_RP")) {curstring >> varname >> ivalue; FLAG_RP = ivalue;}
525 else if(strstr(temp_string.c_str(),"FLAG_trigger")) {curstring >> varname >> ivalue; FLAG_trigger = ivalue;}
526 else if(strstr(temp_string.c_str(),"FLAG_frog")) {curstring >> varname >> ivalue; FLAG_frog = ivalue;}
527 else if(strstr(temp_string.c_str(),"FLAG_lhco")) {curstring >> varname >> ivalue; FLAG_lhco = ivalue;}
528 else if(strstr(temp_string.c_str(),"NEvents_Frog")) {curstring >> varname >> ivalue; NEvents_Frog = ivalue;}
529 }
530
531 //jet stuffs not defined in the input datacard
532 JET_overlap = 0.75;
533 // MidPoint algorithm definition
534 JET_M_coneareafraction = 0.25;
535 JET_M_maxpairsize = 2;
536 JET_M_maxiterations = 100;
537 // Define Cone algorithm.
538 JET_C_adjacencycut = 2;
539 JET_C_maxiterations = 100;
540 JET_C_iratch = 1;
541 //Define SISCone algorithm.
542 JET_S_npass = 0;
543 JET_S_protojet_ptmin= 0.0;
544
545 //For Tau-jet definition
546 TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
547 TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
548 TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
549 TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
550
551}
552
553void RESOLution::Logfile(const string& LogName) {
554 //void RESOLution::Logfile(string outputfilename) {
555
556 ofstream f_out(LogName.c_str());
557
558 f_out <<"**********************************************************************"<< endl;
559 f_out <<"**********************************************************************"<< endl;
560 f_out <<"** **"<< endl;
561 f_out <<"** Welcome to **"<< endl;
562 f_out <<"** **"<< endl;
563 f_out <<"** **"<< endl;
564 f_out <<"** .ddddddd- lL hH **"<< endl;
565 f_out <<"** -Dd` `dD: Ll hH` **"<< endl;
566 f_out <<"** dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss **"<< endl;
567 f_out <<"** -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs **"<< endl;
568 f_out <<"** dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. **"<< endl;
569 f_out <<"** .Dd :dd eE. LlL PpppPP Hh Hh eE sSS **"<< endl;
570 f_out <<"** dddddd:. eee+: lL. pp. hh. hh eee+ sssssS **"<< endl;
571 f_out <<"** Pp **"<< endl;
572 f_out <<"** **"<< endl;
573 f_out <<"** Delphes, a framework for the fast simulation **"<< endl;
574 f_out <<"** of a generic collider experiment **"<< endl;
575 f_out <<"** **"<< endl;
576 f_out <<"** --- Version 1.4beta of Delphes --- **"<< endl;
577 f_out <<"** Last date of change: 9 February 2009 **"<< endl;
578 f_out <<"** **"<< endl;
579 f_out <<"** **"<< endl;
580 f_out <<"** This package uses: **"<< endl;
581 f_out <<"** ------------------ **"<< endl;
582 f_out <<"** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
583 f_out <<"** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **"<< endl;
584 f_out <<"** FROG: L. Quertenmont, V. Roberfroid [hep-ex/0901.2718v1] **"<< endl;
585 f_out <<"** **"<< endl;
586 f_out <<"** ---------------------------------------------------------------- **"<< endl;
587 f_out <<"** **"<< endl;
588 f_out <<"** Main authors: **"<< endl;
589 f_out <<"** ------------- **"<< endl;
590 f_out <<"** **"<< endl;
591 f_out <<"** Séverine Ovyn Xavier Rouby **"<< endl;
592 f_out <<"** severine.ovyn@uclouvain.be xavier.rouby@cern **"<< endl;
593 f_out <<"** Center for Particle Physics and Phenomenology (CP3) **"<< endl;
594 f_out <<"** Universite Catholique de Louvain (UCL) **"<< endl;
595 f_out <<"** Louvain-la-Neuve, Belgium **"<< endl;
596 f_out <<"** **"<< endl;
597 f_out <<"** ---------------------------------------------------------------- **"<< endl;
598 f_out <<"** **"<< endl;
599 f_out <<"** Former Delphes versions and documentation can be found on : **"<< endl;
600 f_out <<"** http://www.fynu.ucl.ac.be/delphes.html **"<< endl;
601 f_out <<"** **"<< endl;
602 f_out <<"** **"<< endl;
603 f_out <<"** Disclaimer: this program is a beta version of Delphes and **"<< endl;
604 f_out <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;
605 f_out <<"** give us your feedbacks about potential bugs **"<< endl;
606 f_out <<"** **"<< endl;
607 f_out <<"**********************************************************************"<< endl;
608 f_out <<"** **"<< endl;
609 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
610 f_out<<"* *"<<"\n";
611 f_out<<"#******************************** *"<<"\n";
612 f_out<<"# Central detector caracteristics *"<<"\n";
613 f_out<<"#******************************** *"<<"\n";
614 f_out<<"* *"<<"\n";
615 f_out << left << setw(30) <<"* Maximum tracking system: "<<""
616 << left << setw(10) <<CEN_max_tracker <<""<< right << setw(15)<<"*"<<"\n";
617 f_out << left << setw(30) <<"* Maximum central calorimeter: "<<""
618 << left << setw(10) <<CEN_max_calo_cen <<""<< right << setw(15)<<"*"<<"\n";
619 f_out << left << setw(30) <<"* Maximum forward calorimeter: "<<""
620 << left << setw(10) <<CEN_max_calo_fwd <<""<< right << setw(15)<<"*"<<"\n";
621 f_out << left << setw(30) <<"* Muon chambers coverage: "<<""
622 << left << setw(10) <<CEN_max_mu <<""<< right << setw(15)<<"*"<<"\n";
623 f_out<<"* *"<<"\n";
624 if(FLAG_RP==1){
625 f_out<<"#************************************ *"<<"\n";
626 f_out<<"# Very forward Roman Pots switched on *"<<"\n";
627 f_out<<"#************************************ *"<<"\n";
628 f_out<<"* *"<<"\n";
629 f_out << left << setw(55) <<"* Distance of the 220 RP to the IP in meters:"<<""
630 << left << setw(5) <<RP_220_s <<""<< right << setw(10)<<"*"<<"\n";
631 f_out << left << setw(55) <<"* Distance of the 220 RP to the beam in meters:"<<""
632 << left << setw(5) <<RP_220_x <<""<< right << setw(10)<<"*"<<"\n";
633 f_out << left << setw(55) <<"* Distance of the 420 RP to the IP in meters:"<<""
634 << left << setw(5) <<RP_420_s <<""<< right << setw(10)<<"*"<<"\n";
635 f_out << left << setw(55) <<"* Distance of the 420 RP to the beam in meters:"<<""
636 << left << setw(5) <<RP_420_x <<""<< right << setw(10)<<"*"<<"\n";
637 f_out << left << setw(55) <<"* Interaction point at the LHC named: "<<""
638 << left << setw(5) <<RP_IP_name <<""<< right << setw(10)<<"*"<<"\n";
639 f_out << left << setw(35) <<"* Datacard for beam 1: "<<""
640 << left << setw(25) <<RP_beam1Card <<""<< right << setw(10)<<"*"<<"\n";
641 f_out << left << setw(35) <<"* Datacard for beam 2: "<<""
642 << left << setw(25) <<RP_beam2Card <<""<< right << setw(10)<<"*"<<"\n";
643 f_out << left << setw(44) <<"* Beam separation, in meters: "<<""
644 << left << setw(6) << RP_offsetEl_x <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
645 f_out << left << setw(44) <<"* Distance from IP for Beam separation (m):"<<""
646 << left << setw(6) <<RP_offsetEl_s <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
647 f_out << left << setw(44) <<"* X offset of beam crossing in micrometers:"<<""
648 << left << setw(6) <<RP_cross_x <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
649 f_out << left << setw(44) <<"* Y offset of beam crossing in micrometers:"<<""
650 << left << setw(6) <<RP_cross_y <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
651 f_out << left << setw(44) <<"* Angle of beam crossing:"<<""
652 << left << setw(6) <<RP_cross_ang <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
653 f_out<<"* *"<<"\n";
654 }
655 else {
656 f_out<<"#************************************* *"<<"\n";
657 f_out<<"# Very forward Roman Pots switched off *"<<"\n";
658 f_out<<"#************************************* *"<<"\n";
659 f_out<<"* *"<<"\n";
660 }
661 if(FLAG_vfd==1){
662 f_out<<"#************************************** *"<<"\n";
663 f_out<<"# Very forward calorimeters switched on *"<<"\n";
664 f_out<<"#************************************** *"<<"\n";
665 f_out<<"* *"<<"\n";
666 f_out << left << setw(55) <<"* Minimum very forward calorimeter: "<<""
667 << left << setw(5) <<VFD_min_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
668 f_out << left << setw(55) <<"* Maximum very forward calorimeter: "<<""
669 << left << setw(5) <<VFD_max_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
670 f_out << left << setw(55) <<"* Minimum coverage zero_degree calorimeter "<<""
671 << left << setw(5) <<VFD_min_zdc <<""<< right << setw(10)<<"*"<<"\n";
672 f_out << left << setw(55) <<"* Distance of the ZDC to the IP, in meters: "<<""
673 << left << setw(5) <<VFD_s_zdc <<""<< right << setw(10)<<"*"<<"\n";
674 f_out<<"* *"<<"\n";
675 }
676 else {
677 f_out<<"#*************************************** *"<<"\n";
678 f_out<<"# Very forward calorimeters switched off *"<<"\n";
679 f_out<<"#*************************************** *"<<"\n";
680 f_out<<"* *"<<"\n";
681 }
682
683 f_out<<"#************************************ *"<<"\n";
684 f_out<<"# Electromagnetic smearing parameters *"<<"\n";
685 f_out<<"#************************************ *"<<"\n";
686 f_out<<"* *"<<"\n";
687 //# \sigma/E = C + N/E + S/\sqrt{E}
688 f_out << left << setw(30) <<"* S term for central ECAL: "<<""
689 << left << setw(30) <<ELG_Scen <<""<< right << setw(10)<<"*"<<"\n";
690 f_out << left << setw(30) <<"* N term for central ECAL: "<<""
691 << left << setw(30) <<ELG_Ncen <<""<< right << setw(10)<<"*"<<"\n";
692 f_out << left << setw(30) <<"* C term for central ECAL: "<<""
693 << left << setw(30) <<ELG_Ccen <<""<< right << setw(10)<<"*"<<"\n";
694 f_out << left << setw(30) <<"* S term for FCAL: "<<""
695 << left << setw(30) <<ELG_Sfwd <<""<< right << setw(10)<<"*"<<"\n";
696 f_out << left << setw(30) <<"* N term for FCAL: "<<""
697 << left << setw(30) <<ELG_Nfwd <<""<< right << setw(10)<<"*"<<"\n";
698 f_out << left << setw(30) <<"* C term for FCAL: "<<""
699 << left << setw(30) <<ELG_Cfwd <<""<< right << setw(10)<<"*"<<"\n";
700 f_out<<"* *"<<"\n";
701 f_out<<"#***************************** *"<<"\n";
702 f_out<<"# Hadronic smearing parameters *"<<"\n";
703 f_out<<"#***************************** *"<<"\n";
704 f_out<<"* *"<<"\n";
705 f_out << left << setw(30) <<"* S term for central HCAL: "<<""
706 << left << setw(30) <<HAD_Shcal <<""<< right << setw(10)<<"*"<<"\n";
707 f_out << left << setw(30) <<"* N term for central HCAL: "<<""
708 << left << setw(30) <<HAD_Nhcal <<""<< right << setw(10)<<"*"<<"\n";
709 f_out << left << setw(30) <<"* C term for central HCAL: "<<""
710 << left << setw(30) <<HAD_Chcal <<""<< right << setw(10)<<"*"<<"\n";
711 f_out << left << setw(30) <<"* S term for FCAL: "<<""
712 << left << setw(30) <<HAD_Shf <<""<< right << setw(10)<<"*"<<"\n";
713 f_out << left << setw(30) <<"* N term for FCAL: "<<""
714 << left << setw(30) <<HAD_Nhf <<""<< right << setw(10)<<"*"<<"\n";
715 f_out << left << setw(30) <<"* C term for FCAL: "<<""
716 << left << setw(30) <<HAD_Chf <<""<< right << setw(10)<<"*"<<"\n";
717 f_out<<"* *"<<"\n";
718 f_out<<"#************************* *"<<"\n";
719 f_out<<"# Muon smearing parameters *"<<"\n";
720 f_out<<"#************************* *"<<"\n";
721 f_out<<"* *"<<"\n";
722 f_out << left << setw(55) <<"* PT resolution for muons : "<<""
723 << left << setw(5) <<MU_SmearPt <<""<< right << setw(10)<<"*"<<"\n";
724 f_out<<"* *"<<"\n";
725 if(FLAG_bfield==1){
726 f_out<<"#*************************** *"<<"\n";
727 f_out<<"# Magnetic field switched on *"<<"\n";
728 f_out<<"#*************************** *"<<"\n";
729 f_out<<"* *"<<"\n";
730 f_out << left << setw(55) <<"* Radius of the BField coverage: "<<""
731 << left << setw(5) <<TRACK_radius <<""<< right << setw(10)<<"*"<<"\n";
732 f_out << left << setw(55) <<"* Length of the BField coverage: "<<""
733 << left << setw(5) <<TRACK_length <<""<< right << setw(10)<<"*"<<"\n";
734 f_out << left << setw(55) <<"* BField X component: "<<""
735 << left << setw(5) <<TRACK_bfield_x <<""<< right << setw(10)<<"*"<<"\n";
736 f_out << left << setw(55) <<"* BField Y component: "<<""
737 << left << setw(5) <<TRACK_bfield_y <<""<< right << setw(10)<<"*"<<"\n";
738 f_out << left << setw(55) <<"* BField Z component: "<<""
739 << left << setw(5) <<TRACK_bfield_z <<""<< right << setw(10)<<"*"<<"\n";
740 f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
741 << left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
742 f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
743 << left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
744 f_out<<"* *"<<"\n";
745 }
746 else {
747 f_out<<"#**************************** *"<<"\n";
748 f_out<<"# Magnetic field switched off *"<<"\n";
749 f_out<<"#**************************** *"<<"\n";
750 f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
751 << left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
752 f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
753 << left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
754 f_out<<"* *"<<"\n";
755 }
756 f_out<<"#******************** *"<<"\n";
757 f_out<<"# Calorimetric Towers *"<<"\n";
758 f_out<<"#******************** *"<<"\n";
759 f_out << left << setw(55) <<"* Number of calorimetric towers in eta, for eta>0: "<<""
760 << left << setw(5) << TOWER_number <<""<< right << setw(10)<<"*"<<"\n";
761 f_out << left << setw(55) <<"* Tower edges in eta, for eta>0: "<<"" << right << setw(15)<<"*"<<"\n";
762 f_out << "* ";
763 for (unsigned int i=0; i<TOWER_number+1; i++) {
764 f_out << left << setw(7) << TOWER_eta_edges[i];
765 if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "* ";
766 }
767 for (unsigned int i=(TOWER_number+1)%9; i<9; i++) f_out << left << setw(7) << "";
768 f_out << right << setw(3)<<"*"<<"\n";
769 f_out << left << setw(55) <<"* Tower sizes in phi, for eta>0 [degree]:"<<"" << right << setw(15)<<"*"<<"\n";
770 f_out << "* ";
771 for (unsigned int i=0; i<TOWER_number; i++) {
772 f_out << left << setw(7) << TOWER_dphi[i];
773 if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "* ";
774 }
775 for (unsigned int i=(TOWER_number)%9; i<9; i++) f_out << left << setw(7) << "";
776 f_out << right << setw(3)<<"*"<<"\n";
777 f_out<<"* *"<<"\n";
778 f_out<<"#******************* *"<<"\n";
779 f_out<<"# Minimum pT's [GeV] *"<<"\n";
780 f_out<<"#******************* *"<<"\n";
781 f_out<<"* *"<<"\n";
782 f_out << left << setw(40) <<"* Minimum pT for electrons: "<<""
783 << left << setw(20) <<PTCUT_elec <<""<< right << setw(10)<<"*"<<"\n";
784 f_out << left << setw(40) <<"* Minimum pT for muons: "<<""
785 << left << setw(20) <<PTCUT_muon <<""<< right << setw(10)<<"*"<<"\n";
786 f_out << left << setw(40) <<"* Minimum pT for jets: "<<""
787 << left << setw(20) <<PTCUT_jet <<""<< right << setw(10)<<"*"<<"\n";
788 f_out << left << setw(40) <<"* Minimum pT for Tau-jets: "<<""
789 << left << setw(20) <<PTCUT_taujet <<""<< right << setw(10)<<"*"<<"\n";
790 f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
791 << left << setw(20) <<PTCUT_gamma <<""<< right << setw(10)<<"*"<<"\n";
792 f_out<<"* *"<<"\n";
793 f_out<<"#******************* *"<<"\n";
794 f_out<<"# Isolation criteria *"<<"\n";
795 f_out<<"#******************* *"<<"\n";
796 f_out<<"* *"<<"\n";
797 f_out << left << setw(40) <<"* Minimum pT for tracks [GeV]: "<<""
798 << left << setw(20) <<ISOL_PT <<""<< right << setw(10)<<"*"<<"\n";
799 f_out << left << setw(40) <<"* Cone for isolation criteria: "<<""
800 << left << setw(20) <<ISOL_Cone <<""<< right << setw(10)<<"*"<<"\n";
801
802 if(ISOL_Calo_ET > 1E98) f_out<<"# No Calorimetric isolation applied *"<<"\n";
803 else {
804 f_out << left << setw(40) <<"* Minimum ET for towers [GeV]: "<<""
805 << left << setw(20) <<ISOL_Calo_ET <<""<< right << setw(10)<<"*"<<"\n";
806 f_out << left << setw(40) <<"* Cone for calorimetric isolation: "<<""
807 << left << setw(20) <<ISOL_Calo_Cone <<""<< right << setw(10)<<"*"<<"\n";
808 f_out << left << setw(40) <<"* Grid size (NxN) for calorimetric isolation: "<<""
809 << left << setw(20) <<ISOL_Calo_Grid <<""<< right << setw(10)<<"*"<<"\n";
810 }
811
812
813 f_out<<"* *"<<"\n";
814 f_out<<"#*************** *"<<"\n";
815 f_out<<"# Jet definition *"<<"\n";
816 f_out<<"#*************** *"<<"\n";
817 f_out<<"* *"<<"\n";
818 f_out<<"* Six algorithms are currently available: *"<<"\n";
819 f_out<<"* - 1) CDF cone algorithm, *"<<"\n";
820 f_out<<"* - 2) CDF MidPoint algorithm, *"<<"\n";
821 f_out<<"* - 3) SIScone algorithm, *"<<"\n";
822 f_out<<"* - 4) kt algorithm, *"<<"\n";
823 f_out<<"* - 5) Cambrigde/Aachen algorithm, *"<<"\n";
824 f_out<<"* - 6) Anti-kt algorithm. *"<<"\n";
825 f_out<<"* *"<<"\n";
826 f_out<<"* You have chosen *"<<"\n";
827 switch(JET_jetalgo) {
828 default:
829 case 1: {
830 f_out<<"* CDF JetClu jet algorithm with parameters: *"<<"\n";
831 f_out << left << setw(40) <<"* - Seed threshold: "<<""
832 << left << setw(10) <<JET_seed <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
833 f_out << left << setw(40) <<"* - Cone radius: "<<""
834 << left << setw(10) <<JET_coneradius <<""<< right << setw(20)<<"*"<<"\n";
835 f_out << left << setw(40) <<"* - Adjacency cut: "<<""
836 << left << setw(10) <<JET_C_adjacencycut <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
837 f_out << left << setw(40) <<"* - Max iterations: "<<""
838 << left << setw(10) <<JET_C_maxiterations <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
839 f_out << left << setw(40) <<"* - Iratch: "<<""
840 << left << setw(10) <<JET_C_iratch <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
841 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
842 << left << setw(10) <<JET_overlap <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
843 }
844 break;
845 case 2: {
846 f_out<<"* CDF midpoint jet algorithm with parameters: *"<<"\n";
847 f_out << left << setw(40) <<"* - Seed threshold: "<<""
848 << left << setw(20) <<JET_seed <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
849 f_out << left << setw(40) <<"* - Cone radius: "<<""
850 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
851 f_out << left << setw(40) <<"* - Cone area fraction:"<<""
852 << left << setw(20) <<JET_M_coneareafraction <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
853 f_out << left << setw(40) <<"* - Maximum pair size: "<<""
854 << left << setw(20) <<JET_M_maxpairsize <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
855 f_out << left << setw(40) <<"* - Max iterations: "<<""
856 << left << setw(20) <<JET_M_maxiterations <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
857 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
858 << left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
859 }
860 break;
861 case 3: {
862 f_out <<"* SISCone jet algorithm with parameters: *"<<"\n";
863 f_out << left << setw(40) <<"* - Cone radius: "<<""
864 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
865 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
866 << left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
867 f_out << left << setw(40) <<"* - Number pass max: "<<""
868 << left << setw(20) <<JET_S_npass <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
869 f_out << left << setw(40) <<"* - Minimum pT for protojet: "<<""
870 << left << setw(20) <<JET_S_protojet_ptmin <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
871 }
872 break;
873 case 4: {
874 f_out <<"* KT jet algorithm with parameters: *"<<"\n";
875 f_out << left << setw(40) <<"* - Cone radius: "<<""
876 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
877 }
878 break;
879 case 5: {
880 f_out <<"* Cambridge/Aachen jet algorithm with parameters: *"<<"\n";
881 f_out << left << setw(40) <<"* - Cone radius: "<<""
882 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
883 }
884 break;
885 case 6: {
886 f_out <<"* Anti-kt jet algorithm with parameters: *"<<"\n";
887 f_out << left << setw(40) <<"* - Cone radius: "<<""
888 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
889 }
890 break;
891 }
892 f_out<<"* *"<<"\n";
893 f_out<<"#****************************** *"<<"\n";
894 f_out<<"# Tau-jet definition parameters *"<<"\n";
895 f_out<<"#****************************** *"<<"\n";
896 f_out<<"* *"<<"\n";
897 f_out << left << setw(45) <<"* Cone radius for calorimeter tagging: "<<""
898 << left << setw(5) <<TAU_energy_scone <<""<< right << setw(20)<<"*"<<"\n";
899 f_out << left << setw(45) <<"* Fraction of energy in the small cone: "<<""
900 << left << setw(5) <<TAU_energy_frac*100 <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
901 f_out << left << setw(45) <<"* Cone radius for tracking tagging: "<<""
902 << left << setw(5) <<TAU_track_scone <<""<< right << setw(20)<<"*"<<"\n";
903 f_out << left << setw(45) <<"* Minimum track pT [GeV]: "<<""
904 << left << setw(5) <<TAU_track_pt <<""<< right << setw(20)<<"*"<<"\n";
905 f_out<<"* *"<<"\n";
906 f_out<<"#*************************** *"<<"\n";
907 f_out<<"# B-tagging efficiencies [%] *"<<"\n";
908 f_out<<"#*************************** *"<<"\n";
909 f_out<<"* *"<<"\n";
910 f_out << left << setw(50) <<"* Efficiency to tag a \"b\" as a b-jet: "<<""
911 << left << setw(10) <<BTAG_b <<""<< right << setw(10)<<"*"<<"\n";
912 f_out << left << setw(50) <<"* Efficiency to mistag a c-jet as a b-jet: "<<""
913 << left << setw(10) <<BTAG_mistag_c <<""<< right << setw(10)<<"*"<<"\n";
914 f_out << left << setw(50) <<"* Efficiency to mistag a light jet as a b-jet: "<<""
915 << left << setw(10) <<BTAG_mistag_l <<""<< right << setw(10)<<"*"<<"\n";
916 f_out<<"* *"<<"\n";
917 f_out<<"* *"<<"\n";
918 f_out<<"#....................................................................*"<<"\n";
919 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
920
921}
922
923// **********Provides the smeared TLorentzVector for the electrons********
924// Smears the electron energy, and changes the 4-momentum accordingly
925// different smearing if the electron is central (eta < 2.5) or forward
926void RESOLution::SmearElectron(TLorentzVector &electron) {
927 // the 'electron' variable will be changed by the function
928 float energy = electron.E(); // before smearing
929 float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
930
931 if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
932 energyS = gRandom->Gaus(energy, sqrt(
933 pow(ELG_Ncen,2) +
934 pow(ELG_Ccen*energy,2) +
935 pow(ELG_Scen*sqrt(energy),2) ));
936 }
937 if(fabs(electron.Eta()) > CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_fwd){
938 energyS = gRandom->Gaus(energy, sqrt(
939 pow(ELG_Nfwd,2) +
940 pow(ELG_Cfwd*energy,2) +
941 pow(ELG_Sfwd*sqrt(energy),2) ) );
942 }
943 electron.SetPtEtaPhiE(energyS/cosh(electron.Eta()), electron.Eta(), electron.Phi(), energyS);
944 if(electron.E() < 0)electron.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
945}
946
947
948// **********Provides the smeared TLorentzVector for the muons********
949// Smears the muon pT and changes the 4-momentum accordingly
950void RESOLution::SmearMu(TLorentzVector &muon) {
951 // the 'muon' variable will be changed by the function
952 float pt = muon.Pt(); // before smearing
953 float ptS=pt;
954
955 if(fabs(muon.Eta()) < CEN_max_mu )
956 {
957 ptS = gRandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
958 }
959 muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
960
961 if(muon.E() < 0)muon.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
962}
963
964
965// **********Provides the smeared TLorentzVector for the hadrons********
966// Smears the hadron 4-momentum
967void RESOLution::SmearHadron(TLorentzVector &hadron, const float frac)
968 // the 'hadron' variable will be changed by the function
969 // the 'frac' variable describes the long-living particles. Should be 0.7 for K0S and Lambda, 1. otherwise
970{
971 float energy = hadron.E(); // before smearing
972 float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
973 float energy_ecal = (1.0 - frac)*energy; // electromagnetic calorimeter
974 float energy_hcal = frac*energy; // hadronic calorimeter
975 // frac takes into account the decay of long-living particles, that decay in the calorimeters
976 // some of the particles decay mostly in the ecal, some mostly in the hcal
977
978 float energyS1,energyS2;
979 if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
980 energyS1 = gRandom->Gaus(energy_hcal, sqrt(
981 pow(HAD_Nhcal,2) +
982 pow(HAD_Chcal*energy_hcal,2) +
983 pow(HAD_Shcal*sqrt(energy_hcal),2) )) ;
984
985
986 energyS2 = gRandom->Gaus(energy_ecal, sqrt(
987 pow(ELG_Ncen,2) +
988 pow(ELG_Ccen*energy_ecal,2) +
989 pow(ELG_Scen*sqrt(energy_ecal),2) ) );
990
991 energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
992 }
993 if(fabs(hadron.Eta()) > CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_fwd){
994 energyS = gRandom->Gaus(energy, sqrt(
995 pow(HAD_Nhf,2) +
996 pow(HAD_Chf*energy,2) +
997 pow(HAD_Shf*sqrt(energy),2) ));
998}
999
1000
1001
1002 hadron.SetPtEtaPhiE(energyS/cosh(hadron.Eta()),hadron.Eta(), hadron.Phi(), energyS);
1003
1004 if(hadron.E() < 0)hadron.SetPxPyPzE(0,0,0,0);
1005}
1006
1007//******************************************************************************************
1008
1009//void RESOLution::SortedVector(vector<ParticleUtil> &vect)
1010void RESOLution::SortedVector(vector<D_Particle> &vect)
1011{
1012 int i,j = 0;
1013 TLorentzVector tmp;
1014 bool en_desordre = true;
1015 int entries=vect.size();
1016 for(i = 0 ; (i < entries) && en_desordre; i++)
1017 {
1018 en_desordre = false;
1019 for(j = 1 ; j < entries - i ; j++)
1020 {
1021 if ( vect[j].Pt() > vect[j-1].Pt() )
1022 {
1023 //ParticleUtil tmp = vect[j-1];
1024 D_Particle tmp = vect[j-1];
1025 vect[j-1] = vect[j];
1026 vect[j] = tmp;
1027 en_desordre = true;
1028 }
1029 }
1030 }
1031}
1032
1033// **********Provides the energy in the cone of radius TAU_CONE_ENERGY for the tau identification********
1034// to be taken into account, a calo tower should
1035// 1) have a transverse energy \f$ E_T = \sqrt{E_X^2 + E_Y^2} \f$ above a given threshold
1036// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
1037double RESOLution::EnergySmallCone(const vector<PhysicsTower> &towers, const float eta, const float phi) {
1038 double Energie=0;
1039 for(unsigned int i=0; i < towers.size(); i++) {
1040 if(towers[i].fourVector.pt() < JET_seed) continue;
1041 if((DeltaR(phi,eta,towers[i].fourVector.phi(),towers[i].fourVector.eta()) < TAU_energy_scone)) {
1042 Energie += towers[i].fourVector.E;
1043 }
1044 }
1045 return Energie;
1046}
1047
1048
1049// **********Provides the number of tracks in the cone of radius TAU_CONE_TRACKS for the tau identification********
1050// to be taken into account, a track should
1051// 1) avec a transverse momentum \$f p_T \$ above a given threshold
1052// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
1053// IMPORTANT REMARK !!!!!
1054// NEW : "charge" will contain the sum of all charged tracks in the cone TAU_track_scone
1055unsigned int RESOLution::NumTracks(float& charge, const vector<TRootTracks> &tracks, const float pt_track, const float eta, const float phi) {
1056 unsigned int numbtrack=0; // number of track in the tau-jet cone, which is smaller than R;
1057 charge=0;
1058 for(unsigned int i=0; i < tracks.size(); i++) {
1059 if(tracks[i].PT < pt_track ) continue;
1060 //float dr = DeltaR(phi,eta,tracks[i].PhiOuter,tracks[i].EtaOuter);
1061 float dr = DeltaR(phi,eta,tracks[i].Phi,tracks[i].Eta);
1062 if (dr > TAU_track_scone) continue;
1063 numbtrack++;
1064 charge += tracks[i].Charge; // total charge in the cone for Tau-jet
1065 }
1066 return numbtrack;
1067}
1068
1069//*** Returns the PID of the particle with the highest energy, in a cone with a radius CONERADIUS and an axis (eta,phi) *********
1070//used by Btaggedjet
1071///// Attention : bug removed => CONERADIUS/2 -> CONERADIUS !!
1072int RESOLution::Bjets(const TSimpleArray<GenParticle> &subarray, const float& eta, const float& phi) {
1073 float emax=0;
1074 int Ppid=0;
1075 if(subarray.GetEntries()>0) {
1076 for(int i=0; i < subarray.GetEntries();i++) { // should have pt>PT_JETMIN and a small cone radius (r<CONE_JET)
1077 float genDeltaR = DeltaR(subarray[i]->Phi,subarray[i]->Eta,phi,eta);
1078 if(genDeltaR < JET_coneradius && subarray[i]->E > emax) {
1079 emax=subarray[i]->E;
1080 Ppid=abs(subarray[i]->PID);
1081 }
1082 }
1083 }
1084 return Ppid;
1085}
1086
1087
1088//******************** Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************
1089bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<GenParticle> &subarray) {
1090 if( rand()%100 < (BTAG_b+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
1091 else if( rand()%100 < (BTAG_mistag_c+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
1092 else if( rand()%100 < (BTAG_mistag_l+1) && Bjets(subarray,JET.Eta(),JET.Phi())!=0) return true; // b-tag of light jets is 1%
1093 return false;
1094}
1095
1096//***********************Isolation criteria***********************
1097//****************************************************************
1098bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, float& ptiso )
1099{
1100 bool isolated = false;
1101 ptiso = 0; // sum of all track pt in isolation cone
1102 float deltar=1E99; // Initial value; should be high; no further repercussion
1103
1104 // loop on all tracks, with p_t above threshold, close enough from the charged lepton
1105 for(unsigned int i=0; i < tracks.size(); i++) {
1106 if(tracks[i].PT < pt_second_track) continue; // ptcut on tracks
1107 float genDeltaR = DeltaR(part.Phi(),part.Eta(),tracks[i].Phi,tracks[i].Eta);
1108 if(
1109 (genDeltaR > deltar) ||
1110 (genDeltaR==0) // rejets the track of the particle itself
1111 ) continue ;
1112 deltar=genDeltaR; // finds the closest track
1113
1114 // as long as (genDeltaR==0) is put above, the particle itself is not taken into account
1115 if( genDeltaR < ISOL_Cone) ptiso += tracks[i].PT; // dR cut on tracks
1116 }
1117 if(deltar > isolCone) isolated = true;
1118 return isolated;
1119}
1120
1121// ******* Calorimetric isolation
1122float RESOLution::CaloIsolation(const D_Particle& part, const D_CaloTowerList & towers) {
1123 // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy
1124 // in a 3×3 grid surrounding the muon to the pT of the muon. For well-isolated muons, both ptiso and etrat will be small.
1125 if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
1126 float etrat=0;
1127 // available parameters: ISOL_Calo_ET , ISOL_Calo_Cone ,
1128/* for(unsigned int i=0; i < towers.size(); i++) {
1129 if(towers[i].E > ISOL_Calo_ET) {
1130 float genDeltaR = DeltaR(part.Phi(),part.Eta(),towers[i].getPhi(),towers[i].getEta());
1131 if(genDeltaR < ISOL_Calo_Cone) {
1132 ptiso += towers[i].getET();
1133 }
1134 }
1135 } // loop on towers
1136 ptiso -=
1137*/
1138 etrat = 100*etrat/part.Pt();
1139 if(etrat<0) cout << "Error: negative etrat in CaloIsolation (" << etrat <<")\n";
1140 else if(etrat>99) cout << "Error: etrat shoud be in [0;99] in CaloIsolation (" << etrat <<")\n";
1141 return etrat;
1142}
1143
1144
1145 //********** returns a segmented value for eta and phi, for calo towers *****
1146void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
1147 iEta = UNDEFINED;
1148 int index= iUNDEFINED;
1149 for (unsigned int i=1; i< TOWER_number+1; i++) {
1150 if(fabs(eta)>TOWER_eta_edges[i-1] && fabs(eta)<TOWER_eta_edges[i]) {
1151 iEta = (eta>0) ? TOWER_eta_edges[i-1] : -TOWER_eta_edges[i];
1152 index = i-1;
1153 break;
1154 }
1155 }
1156 if(index==UNDEFINED) return;
1157 iPhi = UNDEFINED;
1158 float dphi = TOWER_dphi[index]*pi/180.;
1159 for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
1160 float low = -pi+(i-1)*dphi;
1161 float high= low+dphi;
1162 if(phi > low && phi < high ){
1163 iPhi = low;
1164 break;
1165 }
1166 }
1167 if (phi > pi-dphi) iPhi = pi-dphi;
1168}
1169
1170
1171
1172//**************************** Returns the delta Phi ****************************
1173float DeltaPhi(const float phi1, const float phi2) {
1174 float deltaphi=phi1-phi2; // in here, -pi < phi < pi
1175 if(fabs(deltaphi) > pi) {
1176 deltaphi=2.*pi -fabs(deltaphi);// put deltaphi between 0 and pi
1177 }
1178 else deltaphi=fabs(deltaphi);
1179
1180 return deltaphi;
1181}
1182
1183//**************************** Returns the delta R****************************
1184float DeltaR(const float phi1, const float eta1, const float phi2, const float eta2) {
1185 return sqrt(pow(DeltaPhi(phi1,phi2),2) + pow(eta1-eta2,2));
1186}
1187
1188int sign(const int myint) {
1189 if (myint >0) return 1;
1190 else if (myint <0) return -1;
1191 else return 0;
1192}
1193
1194int sign(const float myfloat) {
1195 if (myfloat >0) return 1;
1196 else if (myfloat <0) return -1;
1197 else return 0;
1198}
1199
1200int ChargeVal(const int pid)
1201{
1202 int charge;
1203 if(
1204 (pid == pGAMMA) ||
1205 (pid == pPI0) ||
1206 (pid == pK0L) ||
1207 (pid == pN) ||
1208 (pid == pSIGMA0) ||
1209 (pid == pDELTA0) ||
1210 (pid == pK0S) // not charged particles : invisible by tracker
1211 )
1212 charge = 0;
1213 else charge = (sign(pid));
1214 return charge;
1215
1216}
Note: See TracBrowser for help on using the repository browser.