Fork me on GitHub

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

Last change on this file since 550 was 550, checked in by severine ovyn, 14 years ago

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