Fork me on GitHub

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

Last change on this file since 483 was 472, checked in by Xavier Rouby, 15 years ago

minor change: add the date/time in logfile
File size: 87.5 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \------------- arXiv:0903.2225v1 ------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
12** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
13** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
14** FROG: [hep-ex/0901.2718v1] **
15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
24** **
25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
29** Copyright (C) 2008-2009, **
30** All rights reserved. **
31** **
32***********************************************************************/
33
34/// \file SmearUtil.cc
35/// \brief RESOLution class, and some generic definitions
36
37
38#include "SmearUtil.h"
39#include "TRandom.h"
40#include "TStopwatch.h"
41
42#include <iostream>
43#include <fstream>
44#include <sstream>
45#include <iomanip>
46#include <map>
47#include <vector>
48using namespace std;
49
50//------------------------------------------------------------------------------
51
52RESOLution::RESOLution() {
53
54 // Detector characteristics
55 CEN_max_tracker = 2.5; // Maximum tracker coverage
56 CEN_max_calo_cen = 3.0; // central calorimeter coverage
57 CEN_max_calo_fwd = 5.0; // forward calorimeter pseudorapidity coverage
58 CEN_max_mu = 2.4; // muon chambers pseudorapidity coverage
59
60 // Energy resolution for electron/photon
61 // \sigma/E = C + N/E + S/\sqrt{E}
62 ELG_Scen = 0.05; // S term for central ECAL
63 ELG_Ncen = 0.25; // N term for central ECAL
64 ELG_Ccen = 0.005; // C term for central ECAL
65 ELG_Sfwd = 2.084; // S term for FCAL
66 ELG_Nfwd = 0.0; // N term for FCAL
67 ELG_Cfwd = 0.107; // C term for FCAL
68 ELG_Szdc = 0.70; // S term for ZDC
69 ELG_Nzdc = 0.0; // N term for ZDC
70 ELG_Czdc = 0.08; // C term for ZDC
71
72 // Energy resolution for hadrons in ecal/hcal/hf
73 // \sigma/E = C + N/E + S/\sqrt{E}
74 HAD_Shcal = 1.5; // S term for central HCAL
75 HAD_Nhcal = 0.; // N term for central HCAL
76 HAD_Chcal = 0.05; // C term for central HCAL
77 HAD_Shf = 2.7; // S term for FCAL
78 HAD_Nhf = 0.; // N term for FCAL
79 HAD_Chf = 0.13; // C term for FCAL
80 HAD_Szdc = 1.38; // S term for ZDC
81 HAD_Nzdc = 0.; // N term for ZDC
82 HAD_Czdc = 0.13; // C term for ZDC
83
84 // Muon smearing
85 MU_SmearPt = 0.01;
86
87 // time resolution
88 ZDC_T_resolution = 0; // resolution for time measurement [s]
89 RP220_T_resolution = 0;
90 RP420_T_resolution = 0;
91
92 // Tracking efficiencies
93 TRACK_ptmin = 0.9; // minimal pt needed to reach the calorimeter in GeV
94 TRACK_eff = 100; // efficiency associated to the tracking
95
96 // Calorimetric towers
97 TOWER_number = 40;
98 const float tower_eta_edges[41] = {
99 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,
100 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,
101 4.350, 4.525, 4.700, 5.000}; // temporary object
102 TOWER_eta_edges = new float[TOWER_number+1];
103 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = tower_eta_edges[i];
104
105 const float tower_dphi[40] = {
106 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 10,
107 10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; // temporary object
108 TOWER_dphi = new float[TOWER_number];
109 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = tower_dphi[i];
110
111
112 // Thresholds for reconstructed objetcs (GeV)
113 PTCUT_elec = 10.0;
114 PTCUT_muon = 10.0;
115 PTCUT_jet = 20.0;
116 PTCUT_gamma = 10.0;
117 PTCUT_taujet = 10.0;
118
119 ZDC_gamma_E = 20; // GeV
120 ZDC_n_E = 50; // GeV
121
122 // Isolation
123 ISOL_PT = 2.0; //minimal pt of tracks for isolation criteria
124 ISOL_Cone = 0.5; //Cone for isolation criteria
125 ISOL_Calo_ET = 1E99; //minimal tower energy for isolation criteria. Default off = 1E99
126 ISOL_Calo_Grid = 3; //Grid size (N x N) for calorimetric isolation -- should be odd
127
128 // General jet variable
129 JET_coneradius = 0.7; // generic jet radius ; not for tau's !!!
130 JET_jetalgo = 1; // 1 for Cone algorithm, 2 for MidPoint algorithm, 3 for SIScone algorithm, 4 for kt algorithm
131 JET_seed = 1.0; // minimum seed to start jet reconstruction
132 JET_Eflow = 1; // 1 for Energy flow in jets reco ; 0 if not
133
134 // Tagging definition
135 BTAG_b = 40;
136 BTAG_mistag_c = 10;
137 BTAG_mistag_l = 1;
138
139 // FLAGS
140 FLAG_bfield = 1; //1 to run the bfield propagation else 0
141 FLAG_vfd = 1; //1 to run the very forward detectors else 0
142 FLAG_RP = 1; //1 to run the zero degree calorimeter else 0
143 FLAG_trigger = 1; //1 to run the trigger selection else 0
144 FLAG_frog = 1; //1 to run the FROG event display
145 FLAG_lhco = 1;
146
147 // In case BField propagation allowed
148 TRACK_radius = 129; //radius of the BField coverage
149 TRACK_length = 300; //length of the BField coverage
150 TRACK_bfield_x = 0; //X composant of the BField
151 TRACK_bfield_y = 0; //Y composant of the BField
152 TRACK_bfield_z = 3.8; //Z composant of the BField
153
154 // In case Very forward detectors allowed
155 VFD_min_calo_vfd = 5.2; // very forward calorimeter (if any) like CASTOR
156 VFD_max_calo_vfd = 6.6;
157 VFD_min_zdc = 8.3;
158 VFD_s_zdc = 140; // distance of the Zero Degree Calorimeter, from the Interaction poin, in [m]
159
160 RP_220_s = 220; // distance of the RP to the IP, in meters
161 RP_220_x = 0.002; // distance of the RP to the beam, in meters
162 RP_420_s = 420; // distance of the RP to the IP, in meters
163 RP_420_x = 0.004; // distance of the RP to the beam, in meters
164 RP_IP_name = "IP5";
165 RP_beam1Card = "data/LHCB1IR5_v6.500.tfs";
166 RP_beam2Card = "data/LHCB1IR5_v6.500.tfs";
167
168 // In case FROG event display allowed
169 NEvents_Frog = 10;
170
171 // Number of events to be processed
172 NEvents = -1;
173
174 //********************************************
175 //jet stuffs not defined in the input datacard
176 //********************************************
177
178 JET_overlap = 0.75;
179 // MidPoint algorithm definition
180 JET_M_coneareafraction = 0.25;
181 JET_M_maxpairsize = 2;
182 JET_M_maxiterations = 100;
183 // Define Cone algorithm.
184 JET_C_adjacencycut = 2;
185 JET_C_maxiterations = 100;
186 JET_C_iratch = 1;
187 //Define SISCone algorithm.
188 JET_S_npass = 0;
189 JET_S_protojet_ptmin= 0.0;
190
191 //For Tau-jet definition
192 TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
193 TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
194 TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
195 TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
196
197 PT_QUARKS_MIN = 2.0 ; // minimal pt needed by quarks to do b-tag
198
199 //for very forward detectors
200 RP_offsetEl_s = 120; // distance of beam separation point, from IP
201 RP_offsetEl_x = -0.097; // half distance of separation of beams in horizontal plan, in m
202 RP_offsetEl_y = 0; // half distance of separation of beams in vertical plan, in m
203 RP_cross_x = -500; // IP offset in horizontal plane, in micrometers
204 RP_cross_y = 0.0; // IP offset in vertical plane, in micrometers
205 RP_cross_ang_x = 142.5; // half-crossing angle in horizontal plane, in microrad
206 RP_cross_ang_y = 0.0; // half-crossing angle in vertical plane, in microrad
207
208
209 PdgTableFilename = "data/particle.tbl";
210 inputfilelist = "";
211 detectorcard = "";
212 triggercard = "";
213}
214
215
216RESOLution::RESOLution(const RESOLution & DET) {
217 // Detector characteristics
218 CEN_max_tracker = DET.CEN_max_tracker;
219 CEN_max_calo_cen = DET.CEN_max_calo_cen;
220 CEN_max_calo_fwd = DET.CEN_max_calo_fwd;
221 CEN_max_mu = DET.CEN_max_mu;
222
223 // Energy resolution for electron/photon
224 ELG_Scen = DET.ELG_Scen;
225 ELG_Ncen = DET.ELG_Ncen;
226 ELG_Ccen = DET.ELG_Ccen;
227 ELG_Cfwd = DET.ELG_Cfwd;
228 ELG_Sfwd = DET.ELG_Sfwd;
229 ELG_Nfwd = DET.ELG_Nfwd;
230 ELG_Czdc = DET.ELG_Czdc;
231 ELG_Szdc = DET.ELG_Szdc;
232 ELG_Nzdc = DET.ELG_Nzdc;
233
234 // Energy resolution for hadrons in ecal/hcal/hf/zdc
235 HAD_Shcal = DET.HAD_Shcal;
236 HAD_Nhcal = DET.HAD_Nhcal;
237 HAD_Chcal = DET.HAD_Chcal;
238 HAD_Shf = DET.HAD_Shf;
239 HAD_Nhf = DET.HAD_Nhf;
240 HAD_Chf = DET.HAD_Chf;
241 HAD_Szdc = DET.HAD_Szdc;
242 HAD_Nzdc = DET.HAD_Nzdc;
243 HAD_Czdc = DET.HAD_Czdc;
244
245 // time resolution
246 ZDC_T_resolution = DET.ZDC_T_resolution; // resolution for time measurement [s]
247 RP220_T_resolution = DET.RP220_T_resolution;
248 RP420_T_resolution = DET.RP420_T_resolution;
249
250 // Muon smearing
251 MU_SmearPt = DET.MU_SmearPt;
252
253 // Tracking efficiencies
254 TRACK_ptmin = DET.TRACK_ptmin;
255 TRACK_eff = DET.TRACK_eff;
256
257 // Calorimetric towers
258 TOWER_number = DET.TOWER_number;
259 TOWER_eta_edges = new float[TOWER_number+1];
260 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];
261
262 TOWER_dphi = new float[TOWER_number];
263 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];
264
265 // Thresholds for reconstructed objetcs
266 PTCUT_elec = DET.PTCUT_elec;
267 PTCUT_muon = DET.PTCUT_muon;
268 PTCUT_jet = DET.PTCUT_jet;
269 PTCUT_gamma = DET.PTCUT_gamma;
270 PTCUT_taujet = DET.PTCUT_taujet;
271
272 ZDC_gamma_E = DET.ZDC_gamma_E;
273 ZDC_n_E = DET.ZDC_n_E;
274
275 // Isolation
276 ISOL_PT = DET.ISOL_PT; // tracking isolation
277 ISOL_Cone = DET.ISOL_Cone;
278 ISOL_Calo_ET = DET.ISOL_Calo_ET; // calorimeter isolation, defaut off
279 ISOL_Calo_Grid = DET.ISOL_Calo_Grid;
280
281
282 // General jet variable
283 JET_coneradius = DET.JET_coneradius;
284 JET_jetalgo = DET.JET_jetalgo;
285 JET_seed = DET.JET_seed;
286 JET_Eflow = DET.JET_Eflow;
287
288 // Tagging definition
289 BTAG_b = DET.BTAG_b;
290 BTAG_mistag_c = DET.BTAG_mistag_c;
291 BTAG_mistag_l = DET.BTAG_mistag_l;
292
293 // FLAGS
294 FLAG_bfield = DET.FLAG_bfield;
295 FLAG_vfd = DET.FLAG_vfd;
296 FLAG_RP = DET.FLAG_RP;
297 FLAG_trigger = DET.FLAG_trigger;
298 FLAG_frog = DET.FLAG_frog;
299 FLAG_lhco = DET.FLAG_lhco;
300
301 // In case BField propagation allowed
302 TRACK_radius = DET.TRACK_radius;
303 TRACK_length = DET.TRACK_length;
304 TRACK_bfield_x = DET.TRACK_bfield_x;
305 TRACK_bfield_y = DET.TRACK_bfield_y;
306 TRACK_bfield_z = DET.TRACK_bfield_z;
307
308 // In case Very forward detectors allowed
309 VFD_min_calo_vfd = DET.VFD_min_calo_vfd;
310 VFD_max_calo_vfd = DET.VFD_max_calo_vfd;
311 VFD_min_zdc = DET.VFD_min_zdc;
312 VFD_s_zdc = DET.VFD_s_zdc;
313
314 RP_220_s = DET.RP_220_s;
315 RP_220_x = DET.RP_220_x;
316 RP_420_s = DET.RP_420_s;
317 RP_420_x = DET.RP_420_x;
318 RP_beam1Card = DET.RP_beam1Card;
319 RP_beam2Card = DET.RP_beam2Card;
320 RP_offsetEl_s = DET.RP_offsetEl_s;
321 RP_offsetEl_x = DET.RP_offsetEl_x;
322 RP_offsetEl_y = DET.RP_offsetEl_y;
323 RP_cross_x = DET.RP_cross_x;
324 RP_cross_y = DET.RP_cross_y;
325 RP_cross_ang_x = DET.RP_cross_ang_x;
326 RP_cross_ang_y = DET.RP_cross_ang_y;
327 RP_IP_name = DET.RP_IP_name;
328
329 // In case FROG event display allowed
330 NEvents_Frog = DET.NEvents_Frog;
331
332 // Number of events to be processed
333 NEvents = DET.NEvents;
334
335 JET_overlap = DET.JET_overlap;
336 // MidPoint algorithm definition
337 JET_M_coneareafraction = DET.JET_M_coneareafraction;
338 JET_M_maxpairsize = DET.JET_M_maxpairsize;
339 JET_M_maxiterations = DET.JET_M_maxiterations;
340 // Define Cone algorithm.
341 JET_C_adjacencycut = DET.JET_C_adjacencycut;
342 JET_C_maxiterations = DET.JET_C_maxiterations;
343 JET_C_iratch = DET.JET_C_iratch;
344 //Define SISCone algorithm.
345 JET_S_npass = DET.JET_S_npass;
346 JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin;
347
348 //For Tau-jet definition
349 TAU_energy_scone = DET.TAU_energy_scone;
350 TAU_track_scone = DET.TAU_track_scone;
351 TAU_track_pt = DET.TAU_track_pt;
352 TAU_energy_frac = DET.TAU_energy_frac;
353
354 PT_QUARKS_MIN = DET.PT_QUARKS_MIN;
355 PdgTableFilename = DET.PdgTableFilename;
356 PDGtable = DET.PDGtable;
357 inputfilelist = DET.inputfilelist;
358 detectorcard = DET.detectorcard;
359 triggercard = DET.triggercard;
360}
361
362RESOLution& RESOLution::operator=(const RESOLution& DET) {
363 if(this==&DET) return *this;
364 // Detector characteristics
365 CEN_max_tracker = DET.CEN_max_tracker;
366 CEN_max_calo_cen = DET.CEN_max_calo_cen;
367 CEN_max_calo_fwd = DET.CEN_max_calo_fwd;
368 CEN_max_mu = DET.CEN_max_mu;
369
370 // Energy resolution for electron/photon
371 ELG_Scen = DET.ELG_Scen;
372 ELG_Ncen = DET.ELG_Ncen;
373 ELG_Ccen = DET.ELG_Ccen;
374 ELG_Cfwd = DET.ELG_Cfwd;
375 ELG_Sfwd = DET.ELG_Sfwd;
376 ELG_Nfwd = DET.ELG_Nfwd;
377 ELG_Czdc = DET.ELG_Czdc;
378 ELG_Szdc = DET.ELG_Szdc;
379 ELG_Nzdc = DET.ELG_Nzdc;
380
381 // Energy resolution for hadrons in ecal/hcal/hf
382 HAD_Shcal = DET.HAD_Shcal;
383 HAD_Nhcal = DET.HAD_Nhcal;
384 HAD_Chcal = DET.HAD_Chcal;
385 HAD_Shf = DET.HAD_Shf;
386 HAD_Nhf = DET.HAD_Nhf;
387 HAD_Chf = DET.HAD_Chf;
388 HAD_Szdc = DET.HAD_Szdc;
389 HAD_Nzdc = DET.HAD_Nzdc;
390 HAD_Czdc = DET.HAD_Czdc;
391
392 // time resolution
393 ZDC_T_resolution = DET.ZDC_T_resolution; // resolution for time measurement [s]
394 RP220_T_resolution = DET.RP220_T_resolution;
395 RP420_T_resolution = DET.RP420_T_resolution;
396
397 // Muon smearing
398 MU_SmearPt = DET.MU_SmearPt;
399
400 // Tracking efficiencies
401 TRACK_ptmin = DET.TRACK_ptmin;
402 TRACK_eff = DET.TRACK_eff;
403
404 // Calorimetric towers
405 TOWER_number = DET.TOWER_number;
406 TOWER_eta_edges = new float[TOWER_number+1];
407 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];
408
409 TOWER_dphi = new float[TOWER_number];
410 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];
411
412 // Thresholds for reconstructed objetcs
413 PTCUT_elec = DET.PTCUT_elec;
414 PTCUT_muon = DET.PTCUT_muon;
415 PTCUT_jet = DET.PTCUT_jet;
416 PTCUT_gamma = DET.PTCUT_gamma;
417 PTCUT_taujet = DET.PTCUT_taujet;
418
419 ZDC_gamma_E = DET.ZDC_gamma_E;
420 ZDC_n_E = DET.ZDC_n_E;
421
422 // Isolation
423 ISOL_PT = DET.ISOL_PT; // tracking isolation
424 ISOL_Cone = DET.ISOL_Cone;
425 ISOL_Calo_ET = DET.ISOL_Calo_ET; // calorimeter isolation, defaut off
426 ISOL_Calo_Grid = DET.ISOL_Calo_Grid;
427
428 // General jet variable
429 JET_coneradius = DET.JET_coneradius;
430 JET_jetalgo = DET.JET_jetalgo;
431 JET_seed = DET.JET_seed;
432 JET_Eflow = DET.JET_Eflow;
433
434 // Tagging definition
435 BTAG_b = DET.BTAG_b;
436 BTAG_mistag_c = DET.BTAG_mistag_c;
437 BTAG_mistag_l = DET.BTAG_mistag_l;
438
439 // FLAGS
440 FLAG_bfield = DET.FLAG_bfield;
441 FLAG_vfd = DET.FLAG_vfd;
442 FLAG_RP = DET.FLAG_RP;
443 FLAG_trigger = DET.FLAG_trigger;
444 FLAG_frog = DET.FLAG_frog;
445 FLAG_lhco = DET.FLAG_lhco;
446
447 // In case BField propagation allowed
448 TRACK_radius = DET.TRACK_radius;
449 TRACK_length = DET.TRACK_length;
450 TRACK_bfield_x = DET.TRACK_bfield_x;
451 TRACK_bfield_y = DET.TRACK_bfield_y;
452 TRACK_bfield_z = DET.TRACK_bfield_z;
453
454 // In case Very forward detectors allowed
455 VFD_min_calo_vfd = DET.VFD_min_calo_vfd;
456 VFD_max_calo_vfd = DET.VFD_max_calo_vfd;
457 VFD_min_zdc = DET.VFD_min_zdc;
458 VFD_s_zdc = DET.VFD_s_zdc;
459
460 RP_220_s = DET.RP_220_s;
461 RP_220_x = DET.RP_220_x;
462 RP_420_s = DET.RP_420_s;
463 RP_420_x = DET.RP_420_x;
464 RP_offsetEl_s = DET.RP_offsetEl_s;
465 RP_offsetEl_x = DET.RP_offsetEl_x;
466 RP_offsetEl_y = DET.RP_offsetEl_y;
467 RP_beam1Card = DET.RP_beam1Card;
468 RP_beam2Card = DET.RP_beam2Card;
469 RP_cross_x = DET.RP_cross_x;
470 RP_cross_y = DET.RP_cross_y;
471 RP_cross_ang_x = DET.RP_cross_ang_x;
472 RP_cross_ang_y = DET.RP_cross_ang_y;
473 RP_IP_name = DET.RP_IP_name;
474
475
476 // In case FROG event display allowed
477 NEvents_Frog = DET.NEvents_Frog;
478
479 // Number of events to be processed
480 NEvents = DET.NEvents;
481
482 JET_overlap = DET.JET_overlap;
483 // MidPoint algorithm definition
484 JET_M_coneareafraction = DET.JET_M_coneareafraction;
485 JET_M_maxpairsize = DET.JET_M_maxpairsize;
486 JET_M_maxiterations = DET.JET_M_maxiterations;
487 // Define Cone algorithm.
488 JET_C_adjacencycut = DET.JET_C_adjacencycut;
489 JET_C_maxiterations = DET.JET_C_maxiterations;
490 JET_C_iratch = DET.JET_C_iratch;
491 //Define SISCone algorithm.
492 JET_S_npass = DET.JET_S_npass;
493 JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin;
494
495 //For Tau-jet definition
496 TAU_energy_scone = DET.TAU_energy_scone;
497 TAU_track_scone = DET.TAU_track_scone;
498 TAU_track_pt = DET.TAU_track_pt;
499 TAU_energy_frac = DET.TAU_energy_frac;
500
501 PT_QUARKS_MIN = DET.PT_QUARKS_MIN;
502
503 PdgTableFilename = DET.PdgTableFilename;
504 PDGtable = DET.PDGtable;
505
506 inputfilelist = DET.inputfilelist;
507 detectorcard = DET.detectorcard;
508 triggercard = DET.triggercard;
509
510 return *this;
511}
512
513void RESOLution::setNames(const string& list, const string& det, const string& trig) {
514 inputfilelist = list;
515 detectorcard = det;
516 triggercard = trig;
517}
518
519//------------------------------------------------------------------------------
520void RESOLution::ReadDataCard(const string datacard) {
521
522 string temp_string;
523 istringstream curstring;
524
525 ifstream fichier_a_lire(datacard.c_str());
526 if(!fichier_a_lire.good()) {
527 cout <<"** WARNING: Datadard not found, use default values **" << endl;
528 return;
529 }
530
531 while (getline(fichier_a_lire,temp_string)) {
532 curstring.clear(); // needed when using several times istringstream::str(string)
533 curstring.str(temp_string);
534 string varname;
535 float value; int ivalue; string svalue;
536
537 if(strstr(temp_string.c_str(),"#")) { }
538 else if(strstr(temp_string.c_str(),"CEN_max_tracker")) {curstring >> varname >> value; CEN_max_tracker = value;}
539 else if(strstr(temp_string.c_str(),"CEN_max_calo_cen")) {curstring >> varname >> value; CEN_max_calo_cen = value;}
540 else if(strstr(temp_string.c_str(),"CEN_max_calo_fwd")) {curstring >> varname >> value; CEN_max_calo_fwd = value;}
541 else if(strstr(temp_string.c_str(),"CEN_max_mu")) {curstring >> varname >> value; CEN_max_mu = value;}
542
543 else if(strstr(temp_string.c_str(),"VFD_min_calo_vfd")) {curstring >> varname >> value; VFD_min_calo_vfd = value;}
544 else if(strstr(temp_string.c_str(),"VFD_max_calo_vfd")) {curstring >> varname >> value; VFD_max_calo_vfd = value;}
545 else if(strstr(temp_string.c_str(),"VFD_min_zdc")) {curstring >> varname >> value; VFD_min_zdc = value;}
546 else if(strstr(temp_string.c_str(),"VFD_s_zdc")) {curstring >> varname >> value; VFD_s_zdc = value;}
547
548 else if(strstr(temp_string.c_str(),"RP_220_s")) {curstring >> varname >> value; RP_220_s = value;}
549 else if(strstr(temp_string.c_str(),"RP_220_x")) {curstring >> varname >> value; RP_220_x = value;}
550 else if(strstr(temp_string.c_str(),"RP_420_s")) {curstring >> varname >> value; RP_420_s = value;}
551 else if(strstr(temp_string.c_str(),"RP_420_x")) {curstring >> varname >> value; RP_420_x = value;}
552 else if(strstr(temp_string.c_str(),"RP_beam1Card")) {curstring >> varname >> svalue;RP_beam1Card = svalue;}
553 else if(strstr(temp_string.c_str(),"RP_beam2Card")) {curstring >> varname >> svalue;RP_beam2Card = svalue;}
554 else if(strstr(temp_string.c_str(),"RP_IP_name")) {curstring >> varname >> svalue;RP_IP_name = svalue;}
555
556 else if(strstr(temp_string.c_str(),"RP_offsetEl_s")) {curstring >> varname >> value; RP_offsetEl_s = value;}
557 else if(strstr(temp_string.c_str(),"RP_offsetEl_x")) {curstring >> varname >> value; RP_offsetEl_x = value;}
558 else if(strstr(temp_string.c_str(),"RP_offsetEl_y")) {curstring >> varname >> value; RP_offsetEl_y = value;}
559 else if(strstr(temp_string.c_str(),"RP_cross_x")) {curstring >> varname >> value; RP_cross_x = value;}
560 else if(strstr(temp_string.c_str(),"RP_cross_y")) {curstring >> varname >> value; RP_cross_y = value;}
561 else if(strstr(temp_string.c_str(),"RP_cross_ang_x")) {curstring >> varname >> value; RP_cross_ang_x = value;}
562 else if(strstr(temp_string.c_str(),"RP_cross_ang_y")) {curstring >> varname >> value; RP_cross_ang_y = value;}
563
564 else if(strstr(temp_string.c_str(),"ELG_Scen")) {curstring >> varname >> value; ELG_Scen = value;}
565 else if(strstr(temp_string.c_str(),"ELG_Ncen")) {curstring >> varname >> value; ELG_Ncen = value;}
566 else if(strstr(temp_string.c_str(),"ELG_Ccen")) {curstring >> varname >> value; ELG_Ccen = value;}
567 else if(strstr(temp_string.c_str(),"ELG_Sfwd")) {curstring >> varname >> value; ELG_Sfwd = value;}
568 else if(strstr(temp_string.c_str(),"ELG_Cfwd")) {curstring >> varname >> value; ELG_Cfwd = value;}
569 else if(strstr(temp_string.c_str(),"ELG_Nfwd")) {curstring >> varname >> value; ELG_Nfwd = value;}
570 else if(strstr(temp_string.c_str(),"ELG_Szdc")) {curstring >> varname >> value; ELG_Szdc = value;}
571 else if(strstr(temp_string.c_str(),"ELG_Czdc")) {curstring >> varname >> value; ELG_Czdc = value;}
572 else if(strstr(temp_string.c_str(),"ELG_Nzdc")) {curstring >> varname >> value; ELG_Nzdc = value;}
573
574 else if(strstr(temp_string.c_str(),"HAD_Shcal")) {curstring >> varname >> value; HAD_Shcal = value;}
575 else if(strstr(temp_string.c_str(),"HAD_Nhcal")) {curstring >> varname >> value; HAD_Nhcal = value;}
576 else if(strstr(temp_string.c_str(),"HAD_Chcal")) {curstring >> varname >> value; HAD_Chcal = value;}
577 else if(strstr(temp_string.c_str(),"HAD_Shf")) {curstring >> varname >> value; HAD_Shf = value;}
578 else if(strstr(temp_string.c_str(),"HAD_Nhf")) {curstring >> varname >> value; HAD_Nhf = value;}
579 else if(strstr(temp_string.c_str(),"HAD_Chf")) {curstring >> varname >> value; HAD_Chf = value;}
580 else if(strstr(temp_string.c_str(),"HAD_Szdc")) {curstring >> varname >> value; HAD_Szdc = value;}
581 else if(strstr(temp_string.c_str(),"HAD_Nzdc")) {curstring >> varname >> value; HAD_Nzdc = value;}
582 else if(strstr(temp_string.c_str(),"HAD_Czdc")) {curstring >> varname >> value; HAD_Czdc = value;}
583 else if(strstr(temp_string.c_str(),"ZDC_T_resolution")) {curstring >> varname >> value; ZDC_T_resolution = value;}
584 else if(strstr(temp_string.c_str(),"RP220_T_resolution")) {curstring >> varname >> value; RP220_T_resolution = value;}
585 else if(strstr(temp_string.c_str(),"RP420_T_resolution")) {curstring >> varname >> value; RP420_T_resolution = value;}
586 else if(strstr(temp_string.c_str(),"MU_SmearPt")) {curstring >> varname >> value; MU_SmearPt = value;}
587
588 else if(strstr(temp_string.c_str(),"TRACK_radius")) {curstring >> varname >> ivalue;TRACK_radius = ivalue;}
589 else if(strstr(temp_string.c_str(),"TRACK_length")) {curstring >> varname >> ivalue;TRACK_length = ivalue;}
590 else if(strstr(temp_string.c_str(),"TRACK_bfield_x")) {curstring >> varname >> value; TRACK_bfield_x = value;}
591 else if(strstr(temp_string.c_str(),"TRACK_bfield_y")) {curstring >> varname >> value; TRACK_bfield_y = value;}
592 else if(strstr(temp_string.c_str(),"TRACK_bfield_z")) {curstring >> varname >> value; TRACK_bfield_z = value;}
593 else if(strstr(temp_string.c_str(),"FLAG_bfield")) {curstring >> varname >> ivalue; FLAG_bfield = ivalue;}
594 else if(strstr(temp_string.c_str(),"TRACK_ptmin")) {curstring >> varname >> value; TRACK_ptmin = value;}
595 else if(strstr(temp_string.c_str(),"TRACK_eff")) {curstring >> varname >> ivalue;TRACK_eff = ivalue;}
596
597 else if(strstr(temp_string.c_str(),"TOWER_number")) {curstring >> varname >> ivalue;TOWER_number = ivalue;}
598 else if(strstr(temp_string.c_str(),"TOWER_eta_edges")){
599 curstring >> varname; for(unsigned int i=0; i<TOWER_number+1; i++) {curstring >> value; TOWER_eta_edges[i] = value;} }
600 else if(strstr(temp_string.c_str(),"TOWER_dphi")){
601 curstring >> varname; for(unsigned int i=0; i<TOWER_number; i++) {curstring >> value; TOWER_dphi[i] = value;} }
602
603 else if(strstr(temp_string.c_str(),"PTCUT_elec")) {curstring >> varname >> value; PTCUT_elec = value;}
604 else if(strstr(temp_string.c_str(),"PTCUT_muon")) {curstring >> varname >> value; PTCUT_muon = value;}
605 else if(strstr(temp_string.c_str(),"PTCUT_jet")) {curstring >> varname >> value; PTCUT_jet = value;}
606 else if(strstr(temp_string.c_str(),"PTCUT_gamma")) {curstring >> varname >> value; PTCUT_gamma = value;}
607 else if(strstr(temp_string.c_str(),"PTCUT_taujet")) {curstring >> varname >> value; PTCUT_taujet = value;}
608 else if(strstr(temp_string.c_str(),"ZDC_gamma_E")) {curstring >> varname >> value; ZDC_gamma_E = value;}
609 else if(strstr(temp_string.c_str(),"ZDC_n_E")) {curstring >> varname >> value; ZDC_n_E = value;}
610
611 else if(strstr(temp_string.c_str(),"ISOL_PT")) {curstring >> varname >> value; ISOL_PT = value;}
612 else if(strstr(temp_string.c_str(),"ISOL_Cone")) {curstring >> varname >> value; ISOL_Cone = value;}
613 else if(strstr(temp_string.c_str(),"ISOL_Calo_ET")) {curstring >> varname >> value; ISOL_Calo_ET = value;}
614 else if(strstr(temp_string.c_str(),"ISOL_Calo_Grid")) {curstring >> varname >> ivalue; ISOL_Calo_Grid = ivalue;}
615
616 else if(strstr(temp_string.c_str(),"JET_coneradius")) {curstring >> varname >> value; JET_coneradius = value;}
617 else if(strstr(temp_string.c_str(),"JET_jetalgo")) {curstring >> varname >> ivalue;JET_jetalgo = ivalue;}
618 else if(strstr(temp_string.c_str(),"JET_seed")) {curstring >> varname >> value; JET_seed = value;}
619 else if(strstr(temp_string.c_str(),"JET_Eflow")) {curstring >> varname >> ivalue; JET_Eflow = ivalue;}
620
621 else if(strstr(temp_string.c_str(),"BTAG_b")) {curstring >> varname >> ivalue;BTAG_b = ivalue;}
622 else if(strstr(temp_string.c_str(),"BTAG_mistag_c")) {curstring >> varname >> ivalue;BTAG_mistag_c = ivalue;}
623 else if(strstr(temp_string.c_str(),"BTAG_mistag_l")) {curstring >> varname >> ivalue;BTAG_mistag_l = ivalue;}
624
625 else if(strstr(temp_string.c_str(),"FLAG_vfd")) {curstring >> varname >> ivalue; FLAG_vfd = ivalue;}
626 else if(strstr(temp_string.c_str(),"FLAG_RP")) {curstring >> varname >> ivalue; FLAG_RP = ivalue;}
627 else if(strstr(temp_string.c_str(),"FLAG_trigger")) {curstring >> varname >> ivalue; FLAG_trigger = ivalue;}
628 else if(strstr(temp_string.c_str(),"FLAG_frog")) {curstring >> varname >> ivalue; FLAG_frog = ivalue;}
629 else if(strstr(temp_string.c_str(),"FLAG_lhco")) {curstring >> varname >> ivalue; FLAG_lhco = ivalue;}
630 else if(strstr(temp_string.c_str(),"NEvents_Frog")) {curstring >> varname >> ivalue; NEvents_Frog = ivalue;}
631 else if(strstr(temp_string.c_str(),"NEvents")) {curstring >> varname >> ivalue; NEvents = ivalue;}
632
633 else if(strstr(temp_string.c_str(),"PdgTableFilename")) {curstring >> varname >> svalue; PdgTableFilename = svalue;}
634 }
635
636 if(ISOL_Calo_Grid%2 ==0) {
637 ISOL_Calo_Grid++;
638 cout <<"** WARNING: ISOL_Calo_Grid is not odd. Set it to "<< ISOL_Calo_Grid << " **" << endl;
639 }
640
641 //jet stuffs not defined in the input datacard
642 JET_overlap = 0.75;
643 // MidPoint algorithm definition
644 JET_M_coneareafraction = 0.25;
645 JET_M_maxpairsize = 2;
646 JET_M_maxiterations = 100;
647 // Define Cone algorithm.
648 JET_C_adjacencycut = 2;
649 JET_C_maxiterations = 100;
650 JET_C_iratch = 1;
651 //Define SISCone algorithm.
652 JET_S_npass = 0;
653 JET_S_protojet_ptmin= 0.0;
654
655 //For Tau-jet definition
656 TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
657 TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
658 TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
659 TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
660
661}
662
663void RESOLution::Logfile(const string& LogName) {
664
665 // creates the list of good input files
666 // this list is vector<string> inputfiles.
667 ifstream infile(inputfilelist.c_str());
668 vector<string> inputfiles;
669 string filename;
670 while(1) {
671 infile >> filename; // reads the first line of the list
672 if(!infile.good()) break; // quits when at the end of the list
673 ifstream checking_the_file(filename.c_str()); // try to open the file
674 if(!checking_the_file.good()) continue; // skips bad/unknown files
675 else checking_the_file.close(); // close file if found
676 inputfiles.push_back(filename); // append the name to the vector
677 }
678 infile.close();
679
680 ofstream f_out(LogName.c_str());
681
682 f_out <<"**********************************************************************"<< endl;
683 f_out <<"**********************************************************************"<< endl;
684 f_out <<"** **"<< endl;
685 f_out <<"** Welcome to **"<< endl;
686 f_out <<"** **"<< endl;
687 f_out <<"** **"<< endl;
688 f_out <<"** .ddddddd- lL hH **"<< endl;
689 f_out <<"** -Dd` `dD: Ll hH` **"<< endl;
690 f_out <<"** dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss **"<< endl;
691 f_out <<"** -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs **"<< endl;
692 f_out <<"** dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. **"<< endl;
693 f_out <<"** .Dd :dd eE. LlL PpppPP Hh Hh eE sSS **"<< endl;
694 f_out <<"** dddddd:. eee+: lL. pp. hh. hh eee+ sssssS **"<< endl;
695 f_out <<"** Pp **"<< endl;
696 f_out <<"** **"<< endl;
697 f_out <<"** Delphes, a framework for the fast simulation **"<< endl;
698 f_out <<"** of a generic collider experiment **"<< endl;
699 f_out <<"** **"<< endl;
700 f_out <<"** --- Version 1.7 of Delphes --- **"<< endl;
701 f_out <<"** Last date of change: 7 May 2009 **"<< endl;
702 f_out <<"** **"<< endl;
703 f_out <<"** **"<< endl;
704 f_out <<"** This package uses: **"<< endl;
705 f_out <<"** ------------------ **"<< endl;
706 f_out <<"** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
707 f_out <<"** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **"<< endl;
708 f_out <<"** FROG: L. Quertenmont, V. Roberfroid [hep-ex/0901.2718v1] **"<< endl;
709 f_out <<"** **"<< endl;
710 f_out <<"** ---------------------------------------------------------------- **"<< endl;
711 f_out <<"** **"<< endl;
712 f_out <<"** Main authors: **"<< endl;
713 f_out <<"** ------------- **"<< endl;
714 f_out <<"** **"<< endl;
715 f_out <<"** Séverine Ovyn Xavier Rouby **"<< endl;
716 f_out <<"** severine.ovyn@uclouvain.be xavier.rouby@cern **"<< endl;
717 f_out <<"** Center for Particle Physics and Phenomenology (CP3) **"<< endl;
718 f_out <<"** Universite Catholique de Louvain (UCL) **"<< endl;
719 f_out <<"** Louvain-la-Neuve, Belgium **"<< endl;
720 f_out <<"** **"<< endl;
721 f_out <<"** ---------------------------------------------------------------- **"<< endl;
722 f_out <<"** **"<< endl;
723 f_out <<"** Former Delphes versions and documentation can be found on : **"<< endl;
724 f_out <<"** http://www.fynu.ucl.ac.be/delphes.html **"<< endl;
725 f_out <<"** **"<< endl;
726 f_out <<"** **"<< endl;
727 f_out <<"** Disclaimer: this program is a beta version of Delphes and **"<< endl;
728 f_out <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;
729 f_out <<"** give us your feedbacks about potential bugs **"<< endl;
730 f_out <<"** **"<< endl;
731 f_out <<"**********************************************************************"<< endl;
732 f_out <<"** **"<< endl;
733 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
734 f_out <<"* *"<< endl;
735 f_out <<"#******************************** *"<<"\n";
736 f_out <<"# Input files *"<<"\n";
737 f_out <<"#******************************** *"<<"\n";
738 f_out << left << setw(22) <<"* Input list "<<""
739 << left << setw(39) << inputfilelist << "" << right << setw(9) << "*"<<"\n";
740 for (unsigned int i =0; i<inputfiles.size(); i++) {
741 f_out << left << setw(22) <<"* - file "<<""
742 << left << setw(43) << inputfiles[i] << "" << right << setw(5) << "*"<<"\n";
743 }
744 if(detectorcard != "")
745 f_out << left << setw(22) <<"* Detector card "<<""
746 << left << setw(39) << detectorcard << "" << right << setw(9) << "*"<<"\n";
747 if(triggercard != "")
748 f_out << left << setw(22) <<"* Trigger card "<<""
749 << left << setw(39) << triggercard << "" << right << setw(9) << "*"<<"\n";
750 f_out<<"* Beam optics : *"<<"\n";
751 f_out << left << setw(22) <<"* - beam 1 "<<""
752 << left << setw(33) << RP_beam1Card << "" << right << setw(15) << "*"<<"\n";
753 f_out << left << setw(22) <<"* - beam 2 "<<""
754 << left << setw(33) << RP_beam2Card << "" << right << setw(15) << "*"<<"\n";
755 f_out << left << setw(22) <<"* Input PDG table " << ""
756 << left << setw(39) << PdgTableFilename << "" << right << setw(9) << "*"<<"\n";
757
758 f_out<<"* *"<<"\n";
759 f_out<<"* *"<<"\n";
760 f_out<<"#******************************** *"<<"\n";
761 f_out<<"# Central detector caracteristics *"<<"\n";
762 f_out<<"#******************************** *"<<"\n";
763 f_out<<"* *"<<"\n";
764 f_out << left << setw(30) <<"* Maximum tracking system: "<<""
765 << left << setw(10) <<CEN_max_tracker <<""<< right << setw(15)<<"*"<<"\n";
766 f_out << left << setw(30) <<"* Maximum central calorimeter: "<<""
767 << left << setw(10) <<CEN_max_calo_cen <<""<< right << setw(15)<<"*"<<"\n";
768 f_out << left << setw(30) <<"* Maximum forward calorimeter: "<<""
769 << left << setw(10) <<CEN_max_calo_fwd <<""<< right << setw(15)<<"*"<<"\n";
770 f_out << left << setw(30) <<"* Muon chambers coverage: "<<""
771 << left << setw(10) <<CEN_max_mu <<""<< right << setw(15)<<"*"<<"\n";
772 f_out<<"* *"<<"\n";
773 if(FLAG_RP==1){
774 f_out<<"#************************************ *"<<"\n";
775 f_out<<"# Very forward Roman Pots switched on *"<<"\n";
776 f_out<<"#************************************ *"<<"\n";
777 f_out<<"* *"<<"\n";
778 f_out << left << setw(55) <<"* Distance of the 220 RP to the IP in meters:"<<""
779 << left << setw(5) <<RP_220_s <<""<< right << setw(10)<<"*"<<"\n";
780 f_out << left << setw(55) <<"* Distance of the 220 RP to the beam in meters:"<<""
781 << left << setw(5) <<RP_220_x <<""<< right << setw(10)<<"*"<<"\n";
782 f_out << left << setw(55) <<"* Distance of the 420 RP to the IP in meters:"<<""
783 << left << setw(5) <<RP_420_s <<""<< right << setw(10)<<"*"<<"\n";
784 f_out << left << setw(55) <<"* Distance of the 420 RP to the beam in meters:"<<""
785 << left << setw(5) <<RP_420_x <<""<< right << setw(10)<<"*"<<"\n";
786 f_out << left << setw(55) <<"* Interaction point at the LHC named: "<<""
787 << left << setw(5) <<RP_IP_name <<""<< right << setw(10)<<"*"<<"\n";
788 f_out << left << setw(35) <<"* Datacard for beam 1: "<<""
789 << left << setw(25) <<RP_beam1Card <<""<< right << setw(10)<<"*"<<"\n";
790 f_out << left << setw(35) <<"* Datacard for beam 2: "<<""
791 << left << setw(25) <<RP_beam2Card <<""<< right << setw(10)<<"*"<<"\n";
792 f_out << left << setw(44) <<"* Beam separation, in meters(hor):"<<""
793 << left << setw(6) << RP_offsetEl_x <<""<< right << setw(10)<<"*"<<"\n";
794 f_out << left << setw(44) <<"* Beam separation, in meters(ver):"<<""
795 << left << setw(6) << RP_offsetEl_y <<""<< right << setw(10)<<"*"<<"\n";
796 f_out << left << setw(44) <<"* Distance from IP for Beam separation (m):"<<""
797 << left << setw(6) <<RP_offsetEl_s <<""<< right << setw(10)<<"*"<<"\n";
798 f_out << left << setw(44) <<"* X offset of beam crossing in micrometers:"<<""
799 << left << setw(6) <<RP_cross_x <<""<< right << setw(10)<<"*"<<"\n";
800 f_out << left << setw(44) <<"* Y offset of beam crossing in micrometers:"<<""
801 << left << setw(6) <<RP_cross_y <<""<< right << setw(10)<<"*"<<"\n";
802 f_out << left << setw(44) <<"* X Angle of beam crossing:"<<""
803 << left << setw(6) <<RP_cross_ang_x <<""<< right << setw(10)<<"*"<<"\n";
804 f_out << left << setw(44) <<"* Y Angle of beam crossing:"<<""
805 << left << setw(6) <<RP_cross_ang_y <<""<< right << setw(10)<<"*"<<"\n";
806 f_out<<"* *"<<"\n";
807 }
808 else {
809 f_out<<"#************************************* *"<<"\n";
810 f_out<<"# Very forward Roman Pots switched off *"<<"\n";
811 f_out<<"#************************************* *"<<"\n";
812 f_out<<"* *"<<"\n";
813 }
814 if(FLAG_vfd==1){
815 f_out<<"#************************************** *"<<"\n";
816 f_out<<"# Very forward calorimeters switched on *"<<"\n";
817 f_out<<"#************************************** *"<<"\n";
818 f_out<<"* *"<<"\n";
819 f_out << left << setw(55) <<"* Minimum very forward calorimeter: "<<""
820 << left << setw(5) <<VFD_min_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
821 f_out << left << setw(55) <<"* Maximum very forward calorimeter: "<<""
822 << left << setw(5) <<VFD_max_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
823 f_out << left << setw(55) <<"* Minimum coverage zero_degree calorimeter "<<""
824 << left << setw(5) <<VFD_min_zdc <<""<< right << setw(10)<<"*"<<"\n";
825 f_out << left << setw(55) <<"* Distance of the ZDC to the IP, in meters: "<<""
826 << left << setw(5) <<VFD_s_zdc <<""<< right << setw(10)<<"*"<<"\n";
827 f_out<<"* *"<<"\n";
828 }
829 else {
830 f_out<<"#*************************************** *"<<"\n";
831 f_out<<"# Very forward calorimeters switched off *"<<"\n";
832 f_out<<"#*************************************** *"<<"\n";
833 f_out<<"* *"<<"\n";
834 }
835
836 f_out<<"#************************************ *"<<"\n";
837 f_out<<"# Electromagnetic smearing parameters *"<<"\n";
838 f_out<<"#************************************ *"<<"\n";
839 f_out<<"* *"<<"\n";
840 //# \sigma/E = C + N/E + S/\sqrt{E}
841 f_out << left << setw(30) <<"* S term for central ECAL: "<<""
842 << left << setw(30) <<ELG_Scen <<""<< right << setw(10)<<"*"<<"\n";
843 f_out << left << setw(30) <<"* N term for central ECAL: "<<""
844 << left << setw(30) <<ELG_Ncen <<""<< right << setw(10)<<"*"<<"\n";
845 f_out << left << setw(30) <<"* C term for central ECAL: "<<""
846 << left << setw(30) <<ELG_Ccen <<""<< right << setw(10)<<"*"<<"\n";
847 f_out << left << setw(30) <<"* S term for FCAL: "<<""
848 << left << setw(30) <<ELG_Sfwd <<""<< right << setw(10)<<"*"<<"\n";
849 f_out << left << setw(30) <<"* N term for FCAL: "<<""
850 << left << setw(30) <<ELG_Nfwd <<""<< right << setw(10)<<"*"<<"\n";
851 f_out << left << setw(30) <<"* C term for FCAL: "<<""
852 << left << setw(30) <<ELG_Cfwd <<""<< right << setw(10)<<"*"<<"\n";
853 f_out << left << setw(30) <<"* S term for ZDC: "<<""
854 << left << setw(30) <<ELG_Szdc <<""<< right << setw(10)<<"*"<<"\n";
855 f_out << left << setw(30) <<"* N term for ZDC: "<<""
856 << left << setw(30) <<ELG_Nzdc <<""<< right << setw(10)<<"*"<<"\n";
857 f_out << left << setw(30) <<"* C term for ZDC: "<<""
858 << left << setw(30) <<ELG_Czdc <<""<< right << setw(10)<<"*"<<"\n";
859
860 f_out<<"* *"<<"\n";
861 f_out<<"#***************************** *"<<"\n";
862 f_out<<"# Hadronic smearing parameters *"<<"\n";
863 f_out<<"#***************************** *"<<"\n";
864 f_out<<"* *"<<"\n";
865 f_out << left << setw(30) <<"* S term for central HCAL: "<<""
866 << left << setw(30) <<HAD_Shcal <<""<< right << setw(10)<<"*"<<"\n";
867 f_out << left << setw(30) <<"* N term for central HCAL: "<<""
868 << left << setw(30) <<HAD_Nhcal <<""<< right << setw(10)<<"*"<<"\n";
869 f_out << left << setw(30) <<"* C term for central HCAL: "<<""
870 << left << setw(30) <<HAD_Chcal <<""<< right << setw(10)<<"*"<<"\n";
871 f_out << left << setw(30) <<"* S term for FCAL: "<<""
872 << left << setw(30) <<HAD_Shf <<""<< right << setw(10)<<"*"<<"\n";
873 f_out << left << setw(30) <<"* N term for FCAL: "<<""
874 << left << setw(30) <<HAD_Nhf <<""<< right << setw(10)<<"*"<<"\n";
875 f_out << left << setw(30) <<"* C term for FCAL: "<<""
876 << left << setw(30) <<HAD_Chf <<""<< right << setw(10)<<"*"<<"\n";
877 f_out << left << setw(30) <<"* S term for ZDC: "<<""
878 << left << setw(30) <<HAD_Szdc <<""<< right << setw(10)<<"*"<<"\n";
879 f_out << left << setw(30) <<"* N term for ZDC: "<<""
880 << left << setw(30) <<HAD_Nzdc <<""<< right << setw(10)<<"*"<<"\n";
881 f_out << left << setw(30) <<"* C term for ZDC: "<<""
882 << left << setw(30) <<HAD_Czdc <<""<< right << setw(10)<<"*"<<"\n";
883
884 f_out<<"* *"<<"\n";
885 f_out<<"#************************* *"<<"\n";
886 f_out<<"# Time smearing parameters *"<<"\n";
887 f_out<<"#************************* *"<<"\n";
888 f_out<<"* *"<<"\n";
889 f_out << left << setw(55) <<"* Time resolution for ZDC : "<<""
890 << left << setw(5) <<ZDC_T_resolution <<""<< right << setw(10)<<"*"<<"\n";
891 f_out << left << setw(55) <<"* Time resolution for RP220 : "<<""
892 << left << setw(5) <<RP220_T_resolution <<""<< right << setw(10)<<"*"<<"\n";
893 f_out << left << setw(55) <<"* Time resolution for RP420 : "<<""
894 << left << setw(5) <<RP420_T_resolution <<""<< right << setw(10)<<"*"<<"\n";
895 f_out<<"* *"<<"\n";
896
897 f_out<<"* *"<<"\n";
898 f_out<<"#************************* *"<<"\n";
899 f_out<<"# Muon smearing parameters *"<<"\n";
900 f_out<<"#************************* *"<<"\n";
901 f_out<<"* *"<<"\n";
902 f_out << left << setw(55) <<"* PT resolution for muons : "<<""
903 << left << setw(5) <<MU_SmearPt <<""<< right << setw(10)<<"*"<<"\n";
904 f_out<<"* *"<<"\n";
905 if(FLAG_bfield==1){
906 f_out<<"#*************************** *"<<"\n";
907 f_out<<"# Magnetic field switched on *"<<"\n";
908 f_out<<"#*************************** *"<<"\n";
909 f_out<<"* *"<<"\n";
910 f_out << left << setw(55) <<"* Radius of the BField coverage: "<<""
911 << left << setw(5) <<TRACK_radius <<""<< right << setw(10)<<"*"<<"\n";
912 f_out << left << setw(55) <<"* Length of the BField coverage: "<<""
913 << left << setw(5) <<TRACK_length <<""<< right << setw(10)<<"*"<<"\n";
914 f_out << left << setw(55) <<"* BField X component: "<<""
915 << left << setw(5) <<TRACK_bfield_x <<""<< right << setw(10)<<"*"<<"\n";
916 f_out << left << setw(55) <<"* BField Y component: "<<""
917 << left << setw(5) <<TRACK_bfield_y <<""<< right << setw(10)<<"*"<<"\n";
918 f_out << left << setw(55) <<"* BField Z component: "<<""
919 << left << setw(5) <<TRACK_bfield_z <<""<< right << setw(10)<<"*"<<"\n";
920 f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
921 << left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
922 f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
923 << left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
924 f_out<<"* *"<<"\n";
925 }
926 else {
927 f_out<<"#**************************** *"<<"\n";
928 f_out<<"# Magnetic field switched off *"<<"\n";
929 f_out<<"#**************************** *"<<"\n";
930 f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
931 << left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
932 f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
933 << left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
934 f_out<<"* *"<<"\n";
935 }
936 f_out<<"#******************** *"<<"\n";
937 f_out<<"# Calorimetric Towers *"<<"\n";
938 f_out<<"#******************** *"<<"\n";
939 f_out << left << setw(55) <<"* Number of calorimetric towers in eta, for eta>0: "<<""
940 << left << setw(5) << TOWER_number <<""<< right << setw(10)<<"*"<<"\n";
941 f_out << left << setw(55) <<"* Tower edges in eta, for eta>0: "<<"" << right << setw(15)<<"*"<<"\n";
942 f_out << "* ";
943 for (unsigned int i=0; i<TOWER_number+1; i++) {
944 f_out << left << setw(7) << TOWER_eta_edges[i];
945 if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "* ";
946 }
947 for (unsigned int i=(TOWER_number+1)%9; i<9; i++) f_out << left << setw(7) << "";
948 f_out << right << setw(3)<<"*"<<"\n";
949 f_out << left << setw(55) <<"* Tower sizes in phi, for eta>0 [degree]:"<<"" << right << setw(15)<<"*"<<"\n";
950 f_out << "* ";
951 for (unsigned int i=0; i<TOWER_number; i++) {
952 f_out << left << setw(7) << TOWER_dphi[i];
953 if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "* ";
954 }
955 for (unsigned int i=(TOWER_number)%9; i<9; i++) f_out << left << setw(7) << "";
956 f_out << right << setw(3)<<"*"<<"\n";
957 f_out<<"* *"<<"\n";
958 f_out<<"#******************* *"<<"\n";
959 f_out<<"# Minimum pT's [GeV] *"<<"\n";
960 f_out<<"#******************* *"<<"\n";
961 f_out<<"* *"<<"\n";
962 f_out << left << setw(40) <<"* Minimum pT for electrons: "<<""
963 << left << setw(20) <<PTCUT_elec <<""<< right << setw(10)<<"*"<<"\n";
964 f_out << left << setw(40) <<"* Minimum pT for muons: "<<""
965 << left << setw(20) <<PTCUT_muon <<""<< right << setw(10)<<"*"<<"\n";
966 f_out << left << setw(40) <<"* Minimum pT for jets: "<<""
967 << left << setw(20) <<PTCUT_jet <<""<< right << setw(10)<<"*"<<"\n";
968 f_out << left << setw(40) <<"* Minimum pT for Tau-jets: "<<""
969 << left << setw(20) <<PTCUT_taujet <<""<< right << setw(10)<<"*"<<"\n";
970 f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
971 << left << setw(20) <<PTCUT_gamma <<""<< right << setw(10)<<"*"<<"\n";
972 f_out << left << setw(40) <<"* Minimum E for photons in ZDC: "<<""
973 << left << setw(20) <<ZDC_gamma_E <<""<< right << setw(10)<<"*"<<"\n";
974 f_out << left << setw(40) <<"* Minimum E for neutrons in ZDC: "<<""
975 << left << setw(20) <<ZDC_n_E <<""<< right << setw(10)<<"*"<<"\n";
976
977 f_out<<"* *"<<"\n";
978 f_out<<"#******************* *"<<"\n";
979 f_out<<"# Isolation criteria *"<<"\n";
980 f_out<<"#******************* *"<<"\n";
981 f_out<<"* *"<<"\n";
982 f_out << left << setw(40) <<"* Minimum pT for tracks [GeV]: "<<""
983 << left << setw(20) <<ISOL_PT <<""<< right << setw(10)<<"*"<<"\n";
984 f_out << left << setw(40) <<"* Cone for isolation criteria: "<<""
985 << left << setw(20) <<ISOL_Cone <<""<< right << setw(10)<<"*"<<"\n";
986
987 if(ISOL_Calo_ET > 1E98) f_out<<"# No Calorimetric isolation applied *"<<"\n";
988 else {
989 f_out << left << setw(40) <<"* Minimum ET for towers [GeV]: "<<""
990 << left << setw(20) <<ISOL_Calo_ET <<""<< right << setw(10)<<"*"<<"\n";
991 f_out << left << setw(40) <<"* Grid size (NxN) for calorimetric isolation: "<<""
992 << left << setw(20) <<ISOL_Calo_Grid <<""<< right << setw(4)<<"*"<<"\n";
993 }
994
995
996 f_out<<"* *"<<"\n";
997 f_out<<"#*************** *"<<"\n";
998 f_out<<"# Jet definition *"<<"\n";
999 f_out<<"#*************** *"<<"\n";
1000 if(JET_Eflow)
1001 {
1002 f_out<<"#*************** *"<<"\n";
1003 f_out<<"#* Running considering perfect energy flow on the tracker coverage *"<<"\n";
1004 }
1005 else
1006 {
1007 f_out<<"#* Running considering no energy flow on the tracker coverage *"<<"\n";
1008 f_out<<"#* --> jet algo applied on the calorimetric towers *"<<"\n";
1009 }
1010 f_out<<"* *"<<"\n";
1011 f_out<<"* Six algorithms are currently available: *"<<"\n";
1012 f_out<<"* - 1) CDF cone algorithm, *"<<"\n";
1013 f_out<<"* - 2) CDF MidPoint algorithm, *"<<"\n";
1014 f_out<<"* - 3) SIScone algorithm, *"<<"\n";
1015 f_out<<"* - 4) kt algorithm, *"<<"\n";
1016 f_out<<"* - 5) Cambrigde/Aachen algorithm, *"<<"\n";
1017 f_out<<"* - 6) Anti-kt algorithm. *"<<"\n";
1018 f_out<<"* *"<<"\n";
1019 f_out<<"* You have chosen *"<<"\n";
1020 switch(JET_jetalgo) {
1021 default:
1022 case 1: {
1023 f_out<<"* CDF JetClu jet algorithm with parameters: *"<<"\n";
1024 f_out << left << setw(40) <<"* - Seed threshold: "<<""
1025 << left << setw(10) <<JET_seed <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
1026 f_out << left << setw(40) <<"* - Cone radius: "<<""
1027 << left << setw(10) <<JET_coneradius <<""<< right << setw(20)<<"*"<<"\n";
1028 f_out << left << setw(40) <<"* - Adjacency cut: "<<""
1029 << left << setw(10) <<JET_C_adjacencycut <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
1030 f_out << left << setw(40) <<"* - Max iterations: "<<""
1031 << left << setw(10) <<JET_C_maxiterations <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
1032 f_out << left << setw(40) <<"* - Iratch: "<<""
1033 << left << setw(10) <<JET_C_iratch <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
1034 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
1035 << left << setw(10) <<JET_overlap <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
1036 }
1037 break;
1038 case 2: {
1039 f_out<<"* CDF midpoint jet algorithm with parameters: *"<<"\n";
1040 f_out << left << setw(40) <<"* - Seed threshold: "<<""
1041 << left << setw(20) <<JET_seed <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1042 f_out << left << setw(40) <<"* - Cone radius: "<<""
1043 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
1044 f_out << left << setw(40) <<"* - Cone area fraction:"<<""
1045 << left << setw(20) <<JET_M_coneareafraction <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1046 f_out << left << setw(40) <<"* - Maximum pair size: "<<""
1047 << left << setw(20) <<JET_M_maxpairsize <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1048 f_out << left << setw(40) <<"* - Max iterations: "<<""
1049 << left << setw(20) <<JET_M_maxiterations <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1050 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
1051 << left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1052 }
1053 break;
1054 case 3: {
1055 f_out <<"* SISCone jet algorithm with parameters: *"<<"\n";
1056 f_out << left << setw(40) <<"* - Cone radius: "<<""
1057 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
1058 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
1059 << left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1060 f_out << left << setw(40) <<"* - Number pass max: "<<""
1061 << left << setw(20) <<JET_S_npass <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1062 f_out << left << setw(40) <<"* - Minimum pT for protojet: "<<""
1063 << left << setw(20) <<JET_S_protojet_ptmin <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
1064 }
1065 break;
1066 case 4: {
1067 f_out <<"* KT jet algorithm with parameters: *"<<"\n";
1068 f_out << left << setw(40) <<"* - Cone radius: "<<""
1069 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
1070 }
1071 break;
1072 case 5: {
1073 f_out <<"* Cambridge/Aachen jet algorithm with parameters: *"<<"\n";
1074 f_out << left << setw(40) <<"* - Cone radius: "<<""
1075 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
1076 }
1077 break;
1078 case 6: {
1079 f_out <<"* Anti-kt jet algorithm with parameters: *"<<"\n";
1080 f_out << left << setw(40) <<"* - Cone radius: "<<""
1081 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
1082 }
1083 break;
1084 }
1085 f_out<<"* *"<<"\n";
1086 f_out<<"#****************************** *"<<"\n";
1087 f_out<<"# Tau-jet definition parameters *"<<"\n";
1088 f_out<<"#****************************** *"<<"\n";
1089 f_out<<"* *"<<"\n";
1090 f_out << left << setw(45) <<"* Cone radius for calorimeter tagging: "<<""
1091 << left << setw(5) <<TAU_energy_scone <<""<< right << setw(20)<<"*"<<"\n";
1092 f_out << left << setw(45) <<"* Fraction of energy in the small cone: "<<""
1093 << left << setw(5) <<TAU_energy_frac*100 <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
1094 f_out << left << setw(45) <<"* Cone radius for tracking tagging: "<<""
1095 << left << setw(5) <<TAU_track_scone <<""<< right << setw(20)<<"*"<<"\n";
1096 f_out << left << setw(45) <<"* Minimum track pT [GeV]: "<<""
1097 << left << setw(5) <<TAU_track_pt <<""<< right << setw(20)<<"*"<<"\n";
1098 f_out<<"* *"<<"\n";
1099 f_out<<"#*************************** *"<<"\n";
1100 f_out<<"# B-tagging efficiencies [%] *"<<"\n";
1101 f_out<<"#*************************** *"<<"\n";
1102 f_out<<"* *"<<"\n";
1103 f_out << left << setw(50) <<"* Efficiency to tag a \"b\" as a b-jet: "<<""
1104 << left << setw(10) <<BTAG_b <<""<< right << setw(10)<<"*"<<"\n";
1105 f_out << left << setw(50) <<"* Efficiency to mistag a c-jet as a b-jet: "<<""
1106 << left << setw(10) <<BTAG_mistag_c <<""<< right << setw(10)<<"*"<<"\n";
1107 f_out << left << setw(50) <<"* Efficiency to mistag a light jet as a b-jet: "<<""
1108 << left << setw(10) <<BTAG_mistag_l <<""<< right << setw(10)<<"*"<<"\n";
1109 f_out<<"* *"<<"\n";
1110 f_out<<"* *"<<"\n";
1111 f_out<<"#....................................................................*"<<"\n";
1112 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
1113
1114 f_out.close();
1115}
1116
1117// **********Provides the smeared TLorentzVector for the electrons********
1118// Smears the electron energy, and changes the 4-momentum accordingly
1119// different smearing if the electron is central (eta < 2.5) or forward
1120void RESOLution::SmearElectron(TLorentzVector &electron) {
1121 // the 'electron' variable will be changed by the function
1122 float energy = electron.E(); // before smearing
1123 float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
1124
1125 if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
1126 energyS = gRandom->Gaus(energy, sqrt(
1127 pow(ELG_Ncen,2) +
1128 pow(ELG_Ccen*energy,2) +
1129 pow(ELG_Scen*sqrt(energy),2) ));
1130 }
1131 if(fabs(electron.Eta()) > CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_fwd){
1132 energyS = gRandom->Gaus(energy, sqrt(
1133 pow(ELG_Nfwd,2) +
1134 pow(ELG_Cfwd*energy,2) +
1135 pow(ELG_Sfwd*sqrt(energy),2) ) );
1136 }
1137 electron.SetPtEtaPhiE(energyS/cosh(electron.Eta()), electron.Eta(), electron.Phi(), energyS);
1138 if(electron.E() < 0)electron.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
1139}
1140
1141
1142// **********Provides the smeared TLorentzVector for the muons********
1143// Smears the muon pT and changes the 4-momentum accordingly
1144void RESOLution::SmearMu(TLorentzVector &muon) {
1145 // the 'muon' variable will be changed by the function
1146 float pt = muon.Pt(); // before smearing
1147 float ptS=pt;
1148
1149 if(fabs(muon.Eta()) < CEN_max_mu )
1150 {
1151 ptS = gRandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
1152 }
1153 muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
1154
1155 if(muon.E() < 0)muon.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
1156}
1157
1158
1159// **********Provides the smeared TLorentzVector for the hadrons********
1160// Smears the hadron 4-momentum
1161void RESOLution::SmearHadron(TLorentzVector &hadron, const float frac)
1162 // the 'hadron' variable will be changed by the function
1163 // the 'frac' variable describes the long-living particles. Should be 0.7 for K0S and Lambda, 1. otherwise
1164{
1165 float energy = hadron.E(); // before smearing
1166 float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
1167 float energy_ecal = (1.0 - frac)*energy; // electromagnetic calorimeter
1168 float energy_hcal = frac*energy; // hadronic calorimeter
1169 // frac takes into account the decay of long-living particles, that decay in the calorimeters
1170 // some of the particles decay mostly in the ecal, some mostly in the hcal
1171
1172 float energyS1,energyS2;
1173 if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
1174 energyS1 = gRandom->Gaus(energy_hcal, sqrt(
1175 pow(HAD_Nhcal,2) +
1176 pow(HAD_Chcal*energy_hcal,2) +
1177 pow(HAD_Shcal*sqrt(energy_hcal),2) )) ;
1178
1179
1180 energyS2 = gRandom->Gaus(energy_ecal, sqrt(
1181 pow(ELG_Ncen,2) +
1182 pow(ELG_Ccen*energy_ecal,2) +
1183 pow(ELG_Scen*sqrt(energy_ecal),2) ) );
1184
1185 energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
1186 }
1187 if(fabs(hadron.Eta()) > CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_fwd){
1188 energyS = gRandom->Gaus(energy, sqrt(
1189 pow(HAD_Nhf,2) +
1190 pow(HAD_Chf*energy,2) +
1191 pow(HAD_Shf*sqrt(energy),2) ));
1192}
1193
1194
1195
1196 hadron.SetPtEtaPhiE(energyS/cosh(hadron.Eta()),hadron.Eta(), hadron.Phi(), energyS);
1197
1198 if(hadron.E() < 0)hadron.SetPxPyPzE(0,0,0,0);
1199}
1200
1201//******************************************************************************************
1202
1203//void RESOLution::SortedVector(vector<ParticleUtil> &vect)
1204void RESOLution::SortedVector(vector<D_Particle> &vect)
1205{
1206 int i,j = 0;
1207 TLorentzVector tmp;
1208 bool en_desordre = true;
1209 int entries=vect.size();
1210 for(i = 0 ; (i < entries) && en_desordre; i++)
1211 {
1212 en_desordre = false;
1213 for(j = 1 ; j < entries - i ; j++)
1214 {
1215 if ( vect[j].Pt() > vect[j-1].Pt() )
1216 {
1217 //ParticleUtil tmp = vect[j-1];
1218 D_Particle tmp = vect[j-1];
1219 vect[j-1] = vect[j];
1220 vect[j] = tmp;
1221 en_desordre = true;
1222 }
1223 }
1224 }
1225}
1226
1227// **********Provides the energy in the cone of radius TAU_CONE_ENERGY for the tau identification********
1228// to be taken into account, a calo tower should
1229// 1) have a transverse energy \f$ E_T = \sqrt{E_X^2 + E_Y^2} \f$ above a given threshold
1230// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
1231double RESOLution::EnergySmallCone(const vector<PhysicsTower> &towers, const float eta, const float phi) {
1232 double Energie=0;
1233 for(unsigned int i=0; i < towers.size(); i++) {
1234 if(towers[i].fourVector.pt() < JET_seed) continue;
1235 if((DeltaR(phi,eta,towers[i].fourVector.phi(),towers[i].fourVector.eta()) < TAU_energy_scone)) {
1236 Energie += towers[i].fourVector.E;
1237 }
1238 }
1239 return Energie;
1240}
1241
1242
1243// **********Provides the number of tracks in the cone of radius TAU_CONE_TRACKS for the tau identification********
1244// to be taken into account, a track should
1245// 1) avec a transverse momentum \$f p_T \$ above a given threshold
1246// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
1247// IMPORTANT REMARK !!!!!
1248// NEW : "charge" will contain the sum of all charged tracks in the cone TAU_track_scone
1249unsigned int RESOLution::NumTracks(float& charge, const vector<TRootTracks> &tracks, const float pt_track, const float eta, const float phi) {
1250 unsigned int numbtrack=0; // number of track in the tau-jet cone, which is smaller than R;
1251 charge=0;
1252 for(unsigned int i=0; i < tracks.size(); i++) {
1253 if(tracks[i].PT < pt_track ) continue;
1254 //float dr = DeltaR(phi,eta,tracks[i].PhiOuter,tracks[i].EtaOuter);
1255 float dr = DeltaR(phi,eta,tracks[i].Phi,tracks[i].Eta);
1256 if (dr > TAU_track_scone) continue;
1257 numbtrack++;
1258 charge += tracks[i].Charge; // total charge in the cone for Tau-jet
1259 }
1260 return numbtrack;
1261}
1262
1263//*** Returns the PID of the particle with the highest energy, in a cone with a radius CONERADIUS and an axis (eta,phi) *********
1264//used by Btaggedjet
1265///// Attention : bug removed => CONERADIUS/2 -> CONERADIUS !!
1266int RESOLution::Bjets(const TSimpleArray<TRootC::GenParticle> &subarray, const float& eta, const float& phi) {
1267 float emax=0;
1268 int Ppid=0;
1269 if(subarray.GetEntries()>0) {
1270 for(int i=0; i < subarray.GetEntries();i++) { // should have pt>PT_JETMIN and a small cone radius (r<CONE_JET)
1271 float genDeltaR = DeltaR(subarray[i]->Phi,subarray[i]->Eta,phi,eta);
1272 if(genDeltaR < JET_coneradius && subarray[i]->E > emax) {
1273 emax=subarray[i]->E;
1274 Ppid=abs(subarray[i]->PID);
1275 }
1276 }
1277 }
1278 return Ppid;
1279}
1280
1281
1282//******************** Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************
1283bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootC::GenParticle> &subarray) {
1284 if( rand()%100 < (BTAG_b+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
1285 else if( rand()%100 < (BTAG_mistag_c+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
1286 else if( rand()%100 < (BTAG_mistag_l+1) && Bjets(subarray,JET.Eta(),JET.Phi())!=0) return true; // b-tag of light jets is 1%
1287 return false;
1288}
1289
1290//***********************Isolation criteria***********************
1291//****************************************************************
1292bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, float& ptiso )
1293{
1294 bool isolated = false;
1295 ptiso = 0; // sum of all track pt in isolation cone
1296 float deltar=1E99; // Initial value; should be high; no further repercussion
1297
1298 // loop on all tracks, with p_t above threshold, close enough from the charged lepton
1299 for(unsigned int i=0; i < tracks.size(); i++) {
1300 if(tracks[i].PT < pt_second_track) continue; // ptcut on tracks
1301 float genDeltaR = DeltaR(part.Phi(),part.Eta(),tracks[i].Phi,tracks[i].Eta);
1302 if(
1303 (genDeltaR > deltar) ||
1304 (genDeltaR==0) // rejets the track of the particle itself
1305 ) continue ;
1306 deltar=genDeltaR; // finds the closest track
1307
1308 // as long as (genDeltaR==0) is put above, the particle itself is not taken into account
1309 if( genDeltaR < ISOL_Cone) ptiso += tracks[i].PT; // dR cut on tracks
1310 }
1311 if(deltar > isolCone) isolated = true;
1312 return isolated;
1313}
1314
1315// ******* Calorimetric isolation
1316float RESOLution::CaloIsolation(const D_Particle& part, const D_CaloTowerList & towers, const float iPhi, const float iEta) {
1317 // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy
1318 // in a 3×3 grid surrounding the muon to the pT of the muon. For well-isolated muons, both ptiso and etrat will be small.
1319 if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
1320 float et_sum=0;
1321 // available parameters: ISOL_Calo_ET , ISOL_Calo_Grid
1322 // Get the EtaCalo/PhiCalo of the muon ;
1323 // transform it into iEta/iPhi to get the towers, and their neighbourh (i-1, i-2, etc)
1324
1325 unsigned int N = ISOL_Calo_Grid;
1326 int index= iUNDEFINED; // index of the central tower of the grid in TOWER_eta_edges[.];
1327 // !! TOWER_eta_edges is only with eta>0
1328 // finds the index of the central tower of the NxN grid
1329 for (unsigned int i=1; i< TOWER_number+1; i++) {
1330 if(fabs(iEta) >= TOWER_eta_edges[i-1] && fabs(iEta) < TOWER_eta_edges[i]) {
1331 index = i-1;
1332 break;
1333 }
1334 }
1335 if(index != iUNDEFINED) {
1336 // finds the size in phi of the cells for this eta
1337 float dphi = TOWER_dphi[index]*pi/180.; // in rad
1338
1339 //cout << "Grid " << " ----------\n";
1340 for (unsigned int i_eta=0; i_eta<N; i_eta++) {
1341 unsigned int real_index = (iEta>0) ? index+i_eta-(N-1)/2 : index+1+i_eta-(N-1)/2 ;
1342 float eta_ith_tower = TOWER_eta_edges[real_index];
1343 if(iEta<0) eta_ith_tower *= -1;
1344
1345 for (unsigned int i_phi=0; i_phi<N; i_phi++) {
1346 float phi_ith_tower = iPhi + (float)(i_phi - (N-1.)/2.)*dphi;
1347 D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
1348 if(calMuon.getEta() != UNDEFINED && calMuon.getE() > ISOL_Calo_ET) {
1349 et_sum += calMuon.getE();
1350 //cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\t" << calMuon.getE() << " GeV" << endl;
1351 }
1352 //else cout << "eta/phi = " << eta_ith_tower << "\t" << phi_ith_tower << "\tnot active\n";
1353 }
1354
1355 } // NxN grid
1356 //cout << "-----------" << etrat << endl;
1357 }
1358 else if (CEN_max_mu < CEN_max_calo_fwd)
1359 cout << "** ERROR in RESOLution::CaloIsolation: 'muon'-tower not found! **" << endl;
1360 // should never happen ! this would be a bug
1361 //cout << "etrat = " << et_sum << "\t Pt=" << part.Pt() << endl;
1362
1363 // should return a number between 0 and 99 (due to LHCO definitions)
1364 // which is Pt(muon) / sum(ET)
1365 float etrat = 0.;
1366 if(et_sum==0) etrat = 99.;
1367 else if(et_sum>0) etrat = 100*part.Pt()/et_sum;
1368 if(etrat<0) cout << "Error: negative etrat in CaloIsolation (" << etrat <<")\n";
1369 //else if(etrat>99) cout << "Error: etrat should be in [0;99] in CaloIsolation (" << etrat <<")\n";
1370 if(etrat>99) etrat = 99;
1371 return etrat;
1372}
1373
1374
1375 //********** returns a segmented value for eta and phi, for calo towers *****
1376void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
1377 iEta = UNDEFINED;
1378 int index= iUNDEFINED;
1379 for (unsigned int i=1; i< TOWER_number+1; i++) {
1380 if(fabs(eta)>TOWER_eta_edges[i-1] && fabs(eta)<=TOWER_eta_edges[i]) {
1381 iEta = (eta>0) ? ((TOWER_eta_edges[i-1]+TOWER_eta_edges[i])/2.0) : -((TOWER_eta_edges[i-1]+TOWER_eta_edges[i])/2.0);
1382 index = i-1;
1383 break;
1384 }
1385 }
1386 if(index==UNDEFINED) return;
1387 iPhi = UNDEFINED;
1388 float dphi = TOWER_dphi[index]*pi/180.;
1389 for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
1390 float low = -pi+(i-1)*dphi;
1391 float high= low+dphi;
1392 if(phi > low && phi <= high ){
1393 iPhi = (low+high)/2.0;
1394 break;
1395 }
1396 }
1397 if (phi > pi-dphi) iPhi = pi-dphi;
1398}
1399
1400//**************************** Returns the delta Phi ****************************
1401float DeltaPhi(const float phi1, const float phi2) {
1402 float deltaphi=phi1-phi2; // in here, -pi < phi < pi
1403 if(fabs(deltaphi) > pi) {
1404 deltaphi=2.*pi -fabs(deltaphi);// put deltaphi between 0 and pi
1405 }
1406 else deltaphi=fabs(deltaphi);
1407
1408 return deltaphi;
1409}
1410
1411//**************************** Returns the delta R****************************
1412float DeltaR(const float phi1, const float eta1, const float phi2, const float eta2) {
1413 return sqrt(pow(DeltaPhi(phi1,phi2),2) + pow(eta1-eta2,2));
1414}
1415
1416int sign(const int myint) {
1417 if (myint >0) return 1;
1418 else if (myint <0) return -1;
1419 else return 0;
1420}
1421
1422int sign(const float myfloat) {
1423 if (myfloat >0) return 1;
1424 else if (myfloat <0) return -1;
1425 else return 0;
1426}
1427
1428int ChargeVal(const int pid)
1429{
1430 cout << "ChargeVal :: deprecated function, do not use it anymore" << endl;
1431 int charge;
1432 if(
1433 (pid == pGAMMA) ||
1434 (pid == pPI0) ||
1435 (pid == pK0L) ||
1436 (pid == pN) ||
1437 (pid == pSIGMA0) ||
1438 (pid == pDELTA0) ||
1439 (pid == pK0S) // not charged particles : invisible by tracker
1440 )
1441 charge = 0;
1442 else charge = sign(pid);
1443 return charge;
1444
1445}
1446
1447//------------------------------------------------------------------------------
1448void RESOLution::ReadParticleDataGroupTable() {
1449
1450 string temp_string;
1451 istringstream curstring;
1452
1453 ifstream fichier_a_lire(PdgTableFilename.c_str());
1454 if(!fichier_a_lire.good()) {
1455 cout <<"** ERROR: PDG Table ("<< PdgTableFilename
1456 << ") not found! exit. **" << endl;
1457 exit(1);
1458 return;
1459 }
1460 // first three lines of the file are useless
1461 getline(fichier_a_lire,temp_string);
1462 getline(fichier_a_lire,temp_string);
1463 getline(fichier_a_lire,temp_string);
1464
1465
1466 while (getline(fichier_a_lire,temp_string)) {
1467 curstring.clear(); // needed when using several times istringstream::str(string)
1468 curstring.str(temp_string);
1469 long int ID; std::string name; int charge; float mass; float width; float lifetime;
1470 // ID name chg mass total width lifetime
1471 // 1 d -1 0.33000 0.00000 0.00000E+00
1472 // in the table, the charge is in units of e+/3
1473 // the total width is in GeV
1474 // the lifetime is ctau in mm
1475 curstring >> ID >> name >> charge >> mass >> width >> lifetime;
1476 PdgParticle particle(ID,name,mass,charge/3.,width,lifetime/1000.);
1477 PDGtable.insert(ID,particle);
1478 //PdgTable.insert(pair<int,PdgParticle>(ID,particle));
1479 //cout << PDGtable[ID].name() << "\t" << PDGtable[ID].mass() << "\t" << PDGtable[ID].charge() << endl;
1480 }
1481
1482} // ReadParticleDataGroupTable
1483
1484
1485// to be improved in order to avoid code repetition
1486// sorry, no time to do it right now (XR, 19/05/2009)
1487void time_report(const TStopwatch& global,const TStopwatch& loop,const TStopwatch& trigger,const TStopwatch& frog,const TStopwatch& lhco, const int flag_frog, const int flag_trigger, const int flag_lhco, const string& LogName, const Long64_t allEntries) {
1488
1489TStopwatch globalwatch(global), loopwatch(loop), triggerwatch(trigger), frogwatch(frog), lhcowatch(lhco);
1490
1491 cout <<"** **"<< endl;
1492 cout <<"** ################## Time report ################# **"<< endl;
1493 cout << left << setw(32) <<"** Time report for "<<""
1494 << left << setw(15) << allEntries <<""
1495 << right << setw(22) <<"events **"<<endl;
1496 cout <<"** **"<< endl;
1497 cout << left << setw(10) <<"**"<<""
1498 << left << setw(15) <<"Time (s):"<<""
1499 << right << setw(15) <<"CPU"<<""
1500 << right << setw(15) <<"Real"<<""
1501 << right << setw(14) <<"**"<<endl;
1502 cout << left << setw(10) <<"**"<<""
1503 << left << setw(15) <<" + Global:"<<""
1504 << right << setw(15) <<globalwatch.CpuTime()<<""
1505 << right << setw(15) <<globalwatch.RealTime()<<""
1506 << right << setw(14) <<"**"<<endl;
1507 cout << left << setw(10) <<"**"<<""
1508 << left << setw(15) <<" + Events:"<<""
1509 << right << setw(15) <<loopwatch.CpuTime()<<""
1510 << right << setw(15) <<loopwatch.RealTime()<<""
1511 << right << setw(14) <<"**"<<endl;
1512 if(flag_trigger == 1)
1513 {
1514 cout << left << setw(10) <<"**"<<""
1515 << left << setw(15) <<" + Trigger:"<<""
1516 << right << setw(15) <<triggerwatch.CpuTime()<<""
1517 << right << setw(15) <<triggerwatch.RealTime()<<""
1518 << right << setw(14) <<"**"<<endl;
1519 }
1520 if(flag_frog == 1)
1521 {
1522 cout << left << setw(10) <<"**"<<""
1523 << left << setw(15) <<" + Frog:"<<""
1524 << right << setw(15) <<frogwatch.CpuTime()<<""
1525 << right << setw(15) <<frogwatch.RealTime()<<""
1526 << right << setw(14) <<"**"<<endl;
1527 }
1528 if(flag_lhco == 1)
1529 {
1530 cout << left << setw(10) <<"**"<<""
1531 << left << setw(15) <<" + LHCO:"<<""
1532 << right << setw(15) <<lhcowatch.CpuTime()<<""
1533 << right << setw(15) <<lhcowatch.RealTime()<<""
1534 << right << setw(14) <<"**"<<endl;
1535 }
1536
1537
1538 ofstream f_out(LogName.c_str(),ios_base::app);
1539
1540 f_out <<"** *"<< endl;
1541 f_out <<"** ################## Time report ################# *"<< endl;
1542 f_out << left << setw(32) <<"** Time report for "<<""
1543 << left << setw(15) << allEntries <<""
1544 << right << setw(23) <<"events *"<<endl;
1545 f_out <<"** *"<< endl;
1546 f_out << left << setw(10) <<"**"<<""
1547 << left << setw(15) <<"Time (s):"<<""
1548 << right << setw(15) <<"CPU"<<""
1549 << right << setw(15) <<"Real"<<""
1550 << right << setw(15) <<" *"<<endl;
1551 f_out << left << setw(10) <<"**"<<""
1552 << left << setw(15) <<" + Global:"<<""
1553 << right << setw(15) <<globalwatch.CpuTime()<<""
1554 << right << setw(15) <<globalwatch.RealTime()<<""
1555 << right << setw(15) <<" *"<<endl;
1556 f_out << left << setw(10) <<"**"<<""
1557 << left << setw(15) <<" + Events:"<<""
1558 << right << setw(15) <<loopwatch.CpuTime()<<""
1559 << right << setw(15) <<loopwatch.RealTime()<<""
1560 << right << setw(15) <<" *"<<endl;
1561 if(flag_trigger == 1)
1562 {
1563 f_out << left << setw(10) <<"**"<<""
1564 << left << setw(15) <<" + Trigger:"<<""
1565 << right << setw(15) <<triggerwatch.CpuTime()<<""
1566 << right << setw(15) <<triggerwatch.RealTime()<<""
1567 << right << setw(15) <<" *"<<endl;
1568 }
1569 if(flag_frog == 1)
1570 {
1571 f_out << left << setw(10) <<"**"<<""
1572 << left << setw(15) <<" + Frog:"<<""
1573 << right << setw(15) <<frogwatch.CpuTime()<<""
1574 << right << setw(15) <<frogwatch.RealTime()<<""
1575 << right << setw(15) <<" *"<<endl;
1576 }
1577 if(flag_lhco == 1)
1578 {
1579 f_out << left << setw(10) <<"**"<<""
1580 << left << setw(15) <<" + LHCO:"<<""
1581 << right << setw(15) <<lhcowatch.CpuTime()<<""
1582 << right << setw(15) <<lhcowatch.RealTime()<<""
1583 << right << setw(15) <<" *"<<endl;
1584 }
1585 f_out << left << setw(16) << "** " << ""
1586 << left << setw(52) << get_time_date() << "**" << endl;
1587
1588
1589 f_out<<"* *"<<"\n";
1590 f_out<<"* *"<<"\n";
1591 f_out<<"#....................................................................*"<<"\n";
1592 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
1593
1594 f_out.close();
1595
1596}
1597
1598void print_header() {
1599 cout << endl << endl;
1600
1601 cout <<"*********************************************************************"<< endl;
1602 cout <<"*********************************************************************"<< endl;
1603 cout <<"** **"<< endl;
1604 cout <<"** Welcome to **"<< endl;
1605 cout <<"** **"<< endl;
1606 cout <<"** **"<< endl;
1607 cout <<"** .ddddddd- lL hH **"<< endl;
1608 cout <<"** -Dd` `dD: Ll hH` **"<< endl;
1609 cout <<"** dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss **"<< endl;
1610 cout <<"** -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs **"<< endl;
1611 cout <<"** dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. **"<< endl;
1612 cout <<"** .Dd :dd eE. LlL PpppPP Hh Hh eE sSS **"<< endl;
1613 cout <<"** dddddd:. eee+: lL. pp. hh. hh eee+ sssssS **"<< endl;
1614 cout <<"** Pp **"<< endl;
1615 cout <<"** **"<< endl;
1616 cout <<"** Delphes, a framework for the fast simulation **"<< endl;
1617 cout <<"** of a generic collider experiment **"<< endl;
1618 cout <<"** arXiv:0903.2225v1 [hep-ph] **"<< endl;
1619 cout <<"** **"<< endl;
1620 cout <<"** --- Version 1.7 of Delphes --- **"<< endl;
1621 cout <<"** Last date of change: 7 May 2009 **"<< endl;
1622 cout <<"** **"<< endl;
1623 cout <<"** **"<< endl;
1624 cout <<"** This package uses: **"<< endl;
1625 cout <<"** ------------------ **"<< endl;
1626 cout <<"** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
1627 cout <<"** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **"<< endl;
1628 cout <<"** FROG: [hep-ex/0901.2718v1] **"<< endl;
1629 cout <<"** **"<< endl;
1630 cout <<"**-----------------------------------------------------------------**"<< endl;
1631 cout <<"** **"<< endl;
1632 cout <<"** Main authors: **"<< endl;
1633 cout <<"** ------------- **"<< endl;
1634 cout <<"** **"<< endl;
1635 cout <<"** Séverine Ovyn Xavier Rouby **"<< endl;
1636 cout <<"** severine.ovyn@uclouvain.be xavier.rouby@cern **"<< endl;
1637 cout <<"** Center for Particle Physics and Phenomenology (CP3) **"<< endl;
1638 cout <<"** Universite Catholique de Louvain (UCL) **"<< endl;
1639 cout <<"** Louvain-la-Neuve, Belgium **"<< endl;
1640 cout <<"** **"<< endl;
1641 cout <<"**-----------------------------------------------------------------**"<< endl;
1642 cout <<"** **"<< endl;
1643 cout <<"** Former Delphes versions and documentation can be found on : **"<< endl;
1644 cout <<"** http://www.fynu.ucl.ac.be/delphes.html **"<< endl;
1645 cout <<"** **"<< endl;
1646 cout <<"** **"<< endl;
1647 cout <<"** Disclaimer: this program is a beta version of Delphes and **"<< endl;
1648 cout <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;
1649 cout <<"** give us your feedbacks about potential bugs **"<< endl;
1650 cout <<"** **"<< endl;
1651 cout <<"*********************************************************************"<< endl;
1652 cout <<"*********************************************************************"<< endl;
1653
1654}
1655
1656string get_time_date() {
1657 time_t rawtime;
1658 struct tm * timeinfo;
1659
1660 time ( &rawtime );
1661 timeinfo = localtime ( &rawtime );
1662
1663 char temp[100];
1664 sprintf(temp,"%i/%i/%i %i:%i:%i",timeinfo->tm_mday,timeinfo->tm_mon+1,timeinfo->tm_year+1900,timeinfo->tm_hour,timeinfo->tm_min,timeinfo->tm_sec);
1665 string tempstring(temp);
1666 return tempstring;
1667}
1668
Note: See TracBrowser for help on using the repository browser.