Fork me on GitHub

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

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

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