Fork me on GitHub

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

Last change on this file since 557 was 557, checked in by severine ovyn, 15 years ago

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