Fork me on GitHub

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

Last change on this file since 259 was 257, checked in by severine ovyn, 16 years ago

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