Fork me on GitHub

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

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

update
File size: 96.0 KB
RevLine 
[260]1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
[443]6** \------------- arXiv:0903.2225v1 ------------/ **
[260]7** **
8** **
9** This package uses: **
10** ------------------ **
[443]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] **
[260]14** FROG: [hep-ex/0901.2718v1] **
[443]15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
[260]16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
[443]22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
[260]24** **
[443]25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
[260]29** Copyright (C) 2008-2009, **
[443]30** All rights reserved. **
[260]31** **
32***********************************************************************/
[2]33
34/// \file SmearUtil.cc
35/// \brief RESOLution class, and some generic definitions
36
37
[219]38#include "SmearUtil.h"
[399]39#include "TStopwatch.h"
[2]40
41#include <iostream>
[219]42#include <fstream>
[2]43#include <sstream>
[44]44#include <iomanip>
[380]45#include <map>
[454]46#include <vector>
[526]47#include <cmath>
[537]48#include <cstdlib> // for exit()
[219]49using namespace std;
[44]50
[2]51//------------------------------------------------------------------------------
52
53RESOLution::RESOLution() {
54
[94]55 // Detector characteristics
56 CEN_max_tracker = 2.5; // Maximum tracker coverage
[494]57 CEN_max_calo_cen = 1.7; // central calorimeter coverage
58 CEN_max_calo_ec = 3.0; // calorimeter endcap coverage
[94]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
[494]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
[257]70 ELG_Sfwd = 2.084; // S term for FCAL
[494]71 ELG_Nfwd = 0.0; // N term
72 ELG_Cfwd = 0.107; // C term
[374]73 ELG_Szdc = 0.70; // S term for ZDC
[494]74 ELG_Nzdc = 0.0; // N term
75 ELG_Czdc = 0.08; // C term
[2]76
[494]77 // Energy resolution for hadrons in ecal/hcal/fwd
[94]78 // \sigma/E = C + N/E + S/\sqrt{E}
[494]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
[374]88 HAD_Szdc = 1.38; // S term for ZDC
[494]89 HAD_Nzdc = 0.; // N term
90 HAD_Czdc = 0.13; // C term
[2]91
[94]92 // Muon smearing
93 MU_SmearPt = 0.01;
[2]94
[374]95 // time resolution
96 ZDC_T_resolution = 0; // resolution for time measurement [s]
97 RP220_T_resolution = 0;
98 RP420_T_resolution = 0;
99
[94]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
[2]103
[94]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];
[2]118
119
[374]120 // Thresholds for reconstructed objetcs (GeV)
[94]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;
[33]126
[374]127 ZDC_gamma_E = 20; // GeV
128 ZDC_n_E = 50; // GeV
129
[321]130 // Isolation
[305]131 ISOL_PT = 2.0; //minimal pt of tracks for isolation criteria
132 ISOL_Cone = 0.5; //Cone for isolation criteria
[321]133 ISOL_Calo_ET = 1E99; //minimal tower energy for isolation criteria. Default off = 1E99
[392]134 ISOL_Calo_Grid = 3; //Grid size (N x N) for calorimetric isolation -- should be odd
[305]135
[94]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
[383]140 JET_Eflow = 1; // 1 for Energy flow in jets reco ; 0 if not
[33]141
[94]142 // Tagging definition
[526]143 BTAG_b = 40.;
144 BTAG_mistag_c = 10.;
145 BTAG_mistag_l = 1.;
[2]146
[94]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
[307]150 FLAG_RP = 1; //1 to run the zero degree calorimeter else 0
[94]151 FLAG_trigger = 1; //1 to run the trigger selection else 0
152 FLAG_frog = 1; //1 to run the FROG event display
[307]153 FLAG_lhco = 1;
[2]154
[94]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
[2]161
[94]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]
[2]167
[94]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
[257]172 RP_IP_name = "IP5";
[252]173 RP_beam1Card = "data/LHCB1IR5_v6.500.tfs";
174 RP_beam2Card = "data/LHCB1IR5_v6.500.tfs";
[2]175
[94]176 // In case FROG event display allowed
177 NEvents_Frog = 10;
[2]178
[422]179 // Number of events to be processed
180 NEvents = -1;
181
[94]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
[252]206
207 //for very forward detectors
[399]208 RP_offsetEl_s = 120; // distance of beam separation point, from IP
[404]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
[399]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
[380]215
[399]216
[380]217 PdgTableFilename = "data/particle.tbl";
[454]218 inputfilelist = "";
219 detectorcard = "";
220 triggercard = "";
[526]221
222 grandom = new TRandom3(0); // a new seed is set everytime Delphes is run
[2]223}
224
[219]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;
[494]230 CEN_max_calo_ec = DET.CEN_max_calo_ec;
[219]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;
[494]238 ELG_Sec = DET.ELG_Sec;
239 ELG_Nec = DET.ELG_Nec;
240 ELG_Cec = DET.ELG_Cec;
[219]241 ELG_Cfwd = DET.ELG_Cfwd;
242 ELG_Sfwd = DET.ELG_Sfwd;
243 ELG_Nfwd = DET.ELG_Nfwd;
[374]244 ELG_Czdc = DET.ELG_Czdc;
245 ELG_Szdc = DET.ELG_Szdc;
246 ELG_Nzdc = DET.ELG_Nzdc;
[219]247
[494]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;
[374]258 HAD_Szdc = DET.HAD_Szdc;
259 HAD_Nzdc = DET.HAD_Nzdc;
260 HAD_Czdc = DET.HAD_Czdc;
[219]261
[374]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
[219]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
[374]289 ZDC_gamma_E = DET.ZDC_gamma_E;
290 ZDC_n_E = DET.ZDC_n_E;
291
[321]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;
[305]297
298
[219]299 // General jet variable
300 JET_coneradius = DET.JET_coneradius;
301 JET_jetalgo = DET.JET_jetalgo;
302 JET_seed = DET.JET_seed;
[383]303 JET_Eflow = DET.JET_Eflow;
[219]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;
[306]313 FLAG_RP = DET.FLAG_RP;
[219]314 FLAG_trigger = DET.FLAG_trigger;
315 FLAG_frog = DET.FLAG_frog;
[307]316 FLAG_lhco = DET.FLAG_lhco;
[219]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;
[252]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;
[404]339 RP_offsetEl_y = DET.RP_offsetEl_y;
[254]340 RP_cross_x = DET.RP_cross_x;
341 RP_cross_y = DET.RP_cross_y;
[399]342 RP_cross_ang_x = DET.RP_cross_ang_x;
343 RP_cross_ang_y = DET.RP_cross_ang_y;
[257]344 RP_IP_name = DET.RP_IP_name;
[219]345
346 // In case FROG event display allowed
347 NEvents_Frog = DET.NEvents_Frog;
348
[422]349 // Number of events to be processed
350 NEvents = DET.NEvents;
351
[219]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;
[380]372 PdgTableFilename = DET.PdgTableFilename;
373 PDGtable = DET.PDGtable;
[454]374 inputfilelist = DET.inputfilelist;
375 detectorcard = DET.detectorcard;
376 triggercard = DET.triggercard;
[526]377
378 grandom = new TRandom3(*(DET.grandom));
[219]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;
[494]386 CEN_max_calo_ec = DET.CEN_max_calo_ec;
[219]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;
[494]394 ELG_Sec = DET.ELG_Sec;
395 ELG_Nec = DET.ELG_Nec;
396 ELG_Cec = DET.ELG_Cec;
[219]397 ELG_Cfwd = DET.ELG_Cfwd;
398 ELG_Sfwd = DET.ELG_Sfwd;
399 ELG_Nfwd = DET.ELG_Nfwd;
[374]400 ELG_Czdc = DET.ELG_Czdc;
401 ELG_Szdc = DET.ELG_Szdc;
402 ELG_Nzdc = DET.ELG_Nzdc;
[219]403
[494]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;
[374]414 HAD_Szdc = DET.HAD_Szdc;
415 HAD_Nzdc = DET.HAD_Nzdc;
416 HAD_Czdc = DET.HAD_Czdc;
[219]417
[374]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
[219]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
[374]445 ZDC_gamma_E = DET.ZDC_gamma_E;
446 ZDC_n_E = DET.ZDC_n_E;
447
[321]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;
[305]453
[219]454 // General jet variable
455 JET_coneradius = DET.JET_coneradius;
456 JET_jetalgo = DET.JET_jetalgo;
457 JET_seed = DET.JET_seed;
[383]458 JET_Eflow = DET.JET_Eflow;
[219]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;
[306]468 FLAG_RP = DET.FLAG_RP;
[219]469 FLAG_trigger = DET.FLAG_trigger;
470 FLAG_frog = DET.FLAG_frog;
[307]471 FLAG_lhco = DET.FLAG_lhco;
[219]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;
[252]490 RP_offsetEl_s = DET.RP_offsetEl_s;
491 RP_offsetEl_x = DET.RP_offsetEl_x;
[404]492 RP_offsetEl_y = DET.RP_offsetEl_y;
[252]493 RP_beam1Card = DET.RP_beam1Card;
494 RP_beam2Card = DET.RP_beam2Card;
[254]495 RP_cross_x = DET.RP_cross_x;
496 RP_cross_y = DET.RP_cross_y;
[399]497 RP_cross_ang_x = DET.RP_cross_ang_x;
498 RP_cross_ang_y = DET.RP_cross_ang_y;
[257]499 RP_IP_name = DET.RP_IP_name;
[219]500
[252]501
[219]502 // In case FROG event display allowed
503 NEvents_Frog = DET.NEvents_Frog;
504
[422]505 // Number of events to be processed
506 NEvents = DET.NEvents;
507
[219]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;
[380]528
529 PdgTableFilename = DET.PdgTableFilename;
530 PDGtable = DET.PDGtable;
[454]531
532 inputfilelist = DET.inputfilelist;
533 detectorcard = DET.detectorcard;
534 triggercard = DET.triggercard;
535
[526]536 grandom = new TRandom3(*(DET.grandom));
537
[219]538 return *this;
539}
540
[454]541void RESOLution::setNames(const string& list, const string& det, const string& trig) {
542 inputfilelist = list;
543 detectorcard = det;
544 triggercard = trig;
545}
[219]546
[2]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()) {
[249]555 cout <<"** WARNING: Datadard not found, use default values **" << endl;
[94]556 return;
[2]557 }
[494]558 bool CEN_max_calo_ec_flag = false;
[94]559
[2]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;
[252]564 float value; int ivalue; string svalue;
[494]565
[2]566 if(strstr(temp_string.c_str(),"#")) { }
[94]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;}
[494]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;}
[94]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;}
[494]572
[94]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;}
[257]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;}
[399]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;}
[404]588 else if(strstr(temp_string.c_str(),"RP_offsetEl_y")) {curstring >> varname >> value; RP_offsetEl_y = value;}
[399]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;}
[94]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;}
[494]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;}
[94]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;}
[374]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
[494]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;}
[374]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;}
[494]628
[374]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;}
[94]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;}
[531]641 else if(strstr(temp_string.c_str(),"TRACK_eff")) {curstring >> varname >> value; TRACK_eff = value;}
[33]642
[94]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;} }
[2]648
[94]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;}
[374]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;}
[43]656
[321]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;}
[305]661
[94]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;}
[384]665 else if(strstr(temp_string.c_str(),"JET_Eflow")) {curstring >> varname >> ivalue; JET_Eflow = ivalue;}
[94]666
[526]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;}
[2]670
[94]671 else if(strstr(temp_string.c_str(),"FLAG_vfd")) {curstring >> varname >> ivalue; FLAG_vfd = ivalue;}
[526]672 else if(strstr(temp_string.c_str(),"FLAG_RP")) {curstring >> varname >> ivalue; FLAG_RP = ivalue;}
[94]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;}
[307]675 else if(strstr(temp_string.c_str(),"FLAG_lhco")) {curstring >> varname >> ivalue; FLAG_lhco = ivalue;}
[94]676 else if(strstr(temp_string.c_str(),"NEvents_Frog")) {curstring >> varname >> ivalue; NEvents_Frog = ivalue;}
[422]677 else if(strstr(temp_string.c_str(),"NEvents")) {curstring >> varname >> ivalue; NEvents = ivalue;}
[380]678
679 else if(strstr(temp_string.c_str(),"PdgTableFilename")) {curstring >> varname >> svalue; PdgTableFilename = svalue;}
[94]680 }
[494]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 }
[392]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
[94]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
[2]726}
727
[219]728void RESOLution::Logfile(const string& LogName) {
[454]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
[44]745 ofstream f_out(LogName.c_str());
[260]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;
[528]765 f_out <<"** --- Version 1.8 of Delphes --- **"<< endl;
766 f_out <<"** Last date of change: 16 August 2009 **"<< endl;
[260]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;
[528]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;
[260]795 f_out <<"** **"<< endl;
796 f_out <<"**********************************************************************"<< endl;
797 f_out <<"** **"<< endl;
[380]798 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
[454]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
[44]823 f_out<<"* *"<<"\n";
[380]824 f_out<<"* *"<<"\n";
[44]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: "<<""
[94]830 << left << setw(10) <<CEN_max_tracker <<""<< right << setw(15)<<"*"<<"\n";
[44]831 f_out << left << setw(30) <<"* Maximum central calorimeter: "<<""
[94]832 << left << setw(10) <<CEN_max_calo_cen <<""<< right << setw(15)<<"*"<<"\n";
[494]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: "<<""
[94]836 << left << setw(10) <<CEN_max_calo_fwd <<""<< right << setw(15)<<"*"<<"\n";
[44]837 f_out << left << setw(30) <<"* Muon chambers coverage: "<<""
[94]838 << left << setw(10) <<CEN_max_mu <<""<< right << setw(15)<<"*"<<"\n";
[44]839 f_out<<"* *"<<"\n";
[306]840 if(FLAG_RP==1){
841 f_out<<"#************************************ *"<<"\n";
842 f_out<<"# Very forward Roman Pots switched on *"<<"\n";
843 f_out<<"#************************************ *"<<"\n";
[94]844 f_out<<"* *"<<"\n";
[306]845 f_out << left << setw(55) <<"* Distance of the 220 RP to the IP in meters:"<<""
[94]846 << left << setw(5) <<RP_220_s <<""<< right << setw(10)<<"*"<<"\n";
[306]847 f_out << left << setw(55) <<"* Distance of the 220 RP to the beam in meters:"<<""
[94]848 << left << setw(5) <<RP_220_x <<""<< right << setw(10)<<"*"<<"\n";
[306]849 f_out << left << setw(55) <<"* Distance of the 420 RP to the IP in meters:"<<""
[94]850 << left << setw(5) <<RP_420_s <<""<< right << setw(10)<<"*"<<"\n";
[306]851 f_out << left << setw(55) <<"* Distance of the 420 RP to the beam in meters:"<<""
[94]852 << left << setw(5) <<RP_420_x <<""<< right << setw(10)<<"*"<<"\n";
[257]853 f_out << left << setw(55) <<"* Interaction point at the LHC named: "<<""
854 << left << setw(5) <<RP_IP_name <<""<< right << setw(10)<<"*"<<"\n";
[252]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";
[494]859 f_out << left << setw(54) <<"* Beam separation, in meters(hor):"<<""
[399]860 << left << setw(6) << RP_offsetEl_x <<""<< right << setw(10)<<"*"<<"\n";
[494]861 f_out << left << setw(54) <<"* Beam separation, in meters(ver):"<<""
[404]862 << left << setw(6) << RP_offsetEl_y <<""<< right << setw(10)<<"*"<<"\n";
[494]863 f_out << left << setw(54) <<"* Distance from IP for Beam separation (m):"<<""
[399]864 << left << setw(6) <<RP_offsetEl_s <<""<< right << setw(10)<<"*"<<"\n";
[494]865 f_out << left << setw(54) <<"* X offset of beam crossing in micrometers:"<<""
[399]866 << left << setw(6) <<RP_cross_x <<""<< right << setw(10)<<"*"<<"\n";
[494]867 f_out << left << setw(54) <<"* Y offset of beam crossing in micrometers:"<<""
[399]868 << left << setw(6) <<RP_cross_y <<""<< right << setw(10)<<"*"<<"\n";
[494]869 f_out << left << setw(54) <<"* X Angle of beam crossing:"<<""
[399]870 << left << setw(6) <<RP_cross_ang_x <<""<< right << setw(10)<<"*"<<"\n";
[494]871 f_out << left << setw(54) <<"* Y Angle of beam crossing:"<<""
[399]872 << left << setw(6) <<RP_cross_ang_y <<""<< right << setw(10)<<"*"<<"\n";
[94]873 f_out<<"* *"<<"\n";
874 }
875 else {
[306]876 f_out<<"#************************************* *"<<"\n";
877 f_out<<"# Very forward Roman Pots switched off *"<<"\n";
878 f_out<<"#************************************* *"<<"\n";
[94]879 f_out<<"* *"<<"\n";
880 }
[306]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
[44]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";
[494]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
[257]922 f_out << left << setw(30) <<"* S term for FCAL: "<<""
[44]923 << left << setw(30) <<ELG_Sfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]924 f_out << left << setw(30) <<"* N term for FCAL: "<<""
[44]925 << left << setw(30) <<ELG_Nfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]926 f_out << left << setw(30) <<"* C term for FCAL: "<<""
[44]927 << left << setw(30) <<ELG_Cfwd <<""<< right << setw(10)<<"*"<<"\n";
[374]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
[44]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: "<<""
[494]941 << left << setw(30) <<HAD_Scen <<""<< right << setw(10)<<"*"<<"\n";
[44]942 f_out << left << setw(30) <<"* N term for central HCAL: "<<""
[494]943 << left << setw(30) <<HAD_Ncen <<""<< right << setw(10)<<"*"<<"\n";
[44]944 f_out << left << setw(30) <<"* C term for central HCAL: "<<""
[494]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";
[257]952 f_out << left << setw(30) <<"* S term for FCAL: "<<""
[494]953 << left << setw(30) <<HAD_Sfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]954 f_out << left << setw(30) <<"* N term for FCAL: "<<""
[494]955 << left << setw(30) <<HAD_Nfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]956 f_out << left << setw(30) <<"* C term for FCAL: "<<""
[494]957 << left << setw(30) <<HAD_Cfwd <<""<< right << setw(10)<<"*"<<"\n";
[374]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
[44]965 f_out<<"* *"<<"\n";
966 f_out<<"#************************* *"<<"\n";
[374]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";
[44]980 f_out<<"# Muon smearing parameters *"<<"\n";
981 f_out<<"#************************* *"<<"\n";
982 f_out<<"* *"<<"\n";
[94]983 f_out << left << setw(55) <<"* PT resolution for muons : "<<""
984 << left << setw(5) <<MU_SmearPt <<""<< right << setw(10)<<"*"<<"\n";
[44]985 f_out<<"* *"<<"\n";
[94]986 if(FLAG_bfield==1){
987 f_out<<"#*************************** *"<<"\n";
[264]988 f_out<<"# Magnetic field switched on *"<<"\n";
[94]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";
[264]1009 f_out<<"# Magnetic field switched off *"<<"\n";
[94]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";
[44]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: "<<""
[94]1044 << left << setw(20) <<PTCUT_elec <<""<< right << setw(10)<<"*"<<"\n";
[44]1045 f_out << left << setw(40) <<"* Minimum pT for muons: "<<""
[94]1046 << left << setw(20) <<PTCUT_muon <<""<< right << setw(10)<<"*"<<"\n";
[44]1047 f_out << left << setw(40) <<"* Minimum pT for jets: "<<""
[94]1048 << left << setw(20) <<PTCUT_jet <<""<< right << setw(10)<<"*"<<"\n";
[44]1049 f_out << left << setw(40) <<"* Minimum pT for Tau-jets: "<<""
[94]1050 << left << setw(20) <<PTCUT_taujet <<""<< right << setw(10)<<"*"<<"\n";
[74]1051 f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
[94]1052 << left << setw(20) <<PTCUT_gamma <<""<< right << setw(10)<<"*"<<"\n";
[374]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
[44]1058 f_out<<"* *"<<"\n";
[305]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";
[321]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: "<<""
[399]1073 << left << setw(20) <<ISOL_Calo_Grid <<""<< right << setw(4)<<"*"<<"\n";
[321]1074 }
1075
1076
[305]1077 f_out<<"* *"<<"\n";
[44]1078 f_out<<"#*************** *"<<"\n";
1079 f_out<<"# Jet definition *"<<"\n";
1080 f_out<<"#*************** *"<<"\n";
[383]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 }
[44]1091 f_out<<"* *"<<"\n";
[49]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";
[94]1101 switch(JET_jetalgo) {
[44]1102 default:
1103 case 1: {
[94]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";
[44]1117 }
1118 break;
1119 case 2: {
[94]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";
[44]1133 }
1134 break;
1135 case 3: {
[94]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";
[44]1145 }
1146 break;
1147 case 4: {
[94]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";
[44]1151 }
1152 break;
[49]1153 case 5: {
[94]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";
[44]1157 }
[49]1158 break;
1159 case 6: {
[94]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";
[49]1163 }
1164 break;
1165 }
[44]1166 f_out<<"* *"<<"\n";
[94]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";
[44]1192 f_out<<"#....................................................................*"<<"\n";
1193 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
[399]1194
1195 f_out.close();
[44]1196}
1197
[2]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}
[71]1205
[94]1206 if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
[526]1207 energyS = grandom->Gaus(energy, sqrt(
[2]1208 pow(ELG_Ncen,2) +
1209 pow(ELG_Ccen*energy,2) +
[22]1210 pow(ELG_Scen*sqrt(energy),2) ));
[55]1211 }
[494]1212 else if(fabs(electron.Eta()) >= CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_ec){
[526]1213 energyS = grandom->Gaus(energy, sqrt(
[494]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){
[526]1219 energyS = grandom->Gaus(energy, sqrt(
[2]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
[61]1234 float ptS=pt;
1235
[94]1236 if(fabs(muon.Eta()) < CEN_max_mu )
[61]1237 {
[526]1238 ptS = grandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
[61]1239 }
1240 muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
[2]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
[31]1259 float energyS1,energyS2;
[94]1260 if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
[526]1261 energyS1 = grandom->Gaus(energy_hcal, sqrt(
[494]1262 pow(HAD_Ncen,2) +
1263 pow(HAD_Ccen*energy_hcal,2) +
1264 pow(HAD_Scen*sqrt(energy_hcal),2) )) ;
[9]1265
[526]1266 energyS2 = grandom->Gaus(energy_ecal, sqrt(
[32]1267 pow(ELG_Ncen,2) +
1268 pow(ELG_Ccen*energy_ecal,2) +
1269 pow(ELG_Scen*sqrt(energy_ecal),2) ) );
[9]1270
[10]1271 energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
[55]1272 }
[494]1273 else if(fabs(hadron.Eta()) >= CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_ec) {
[526]1274 energyS1 = grandom->Gaus(energy_hcal, sqrt(
[494]1275 pow(HAD_Nec,2) +
1276 pow(HAD_Cec*energy_hcal,2) +
1277 pow(HAD_Sec*sqrt(energy_hcal),2) )) ;
1278
[526]1279 energyS2 = grandom->Gaus(energy_ecal, sqrt(
[494]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){
[526]1287 energyS = grandom->Gaus(energy, sqrt(
[494]1288 pow(HAD_Nfwd,2) +
1289 pow(HAD_Cfwd*energy,2) +
1290 pow(HAD_Sfwd*sqrt(energy),2) ));
[55]1291}
1292
[10]1293
1294
[2]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
[74]1300//******************************************************************************************
1301
[264]1302//void RESOLution::SortedVector(vector<ParticleUtil> &vect)
1303void RESOLution::SortedVector(vector<D_Particle> &vect)
[74]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 {
[264]1316 //ParticleUtil tmp = vect[j-1];
1317 D_Particle tmp = vect[j-1];
[74]1318 vect[j-1] = vect[j];
1319 vect[j] = tmp;
1320 en_desordre = true;
1321 }
1322 }
1323 }
1324}
1325
[2]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++) {
[94]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)) {
[2]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 !!!!!
[287]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;
[2]1351 for(unsigned int i=0; i < tracks.size(); i++) {
[287]1352 if(tracks[i].PT < pt_track ) continue;
[319]1353 //float dr = DeltaR(phi,eta,tracks[i].PhiOuter,tracks[i].EtaOuter);
[287]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
[2]1358 }
[287]1359 return numbtrack;
[2]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 !!
[350]1365int RESOLution::Bjets(const TSimpleArray<TRootC::GenParticle> &subarray, const float& eta, const float& phi) {
[2]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);
[94]1371 if(genDeltaR < JET_coneradius && subarray[i]->E > emax) {
[2]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****************
[526]1382
[350]1383bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootC::GenParticle> &subarray) {
[526]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;
[2]1413}
1414
[526]1415
[31]1416//***********************Isolation criteria***********************
1417//****************************************************************
[321]1418bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, float& ptiso )
[31]1419{
1420 bool isolated = false;
[321]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);
[31]1428 if(
1429 (genDeltaR > deltar) ||
[321]1430 (genDeltaR==0) // rejets the track of the particle itself
[31]1431 ) continue ;
[321]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
[31]1436 }
[305]1437 if(deltar > isolCone) isolated = true;
[31]1438 return isolated;
1439}
1440
[556]1441float RESOLution::IsolationSumEt(const float iPhi, const float iEta, D_CaloTowerList & towers)
1442{
1443 float etiso = 0; // sum of all track pt in isolation cone
1444cout<<"*************************"<<endl;
1445 for (unsigned int i=0; i< towers.size(); i++) {
1446// cout<<" eta part "<<iEta<<" eta tour "<<towers[i].getEta()<<endl;
1447// cout<<" phi part "<<iPhi<<" phi tour "<<towers[i].getPhi()<<endl;
1448
1449 if(DeltaR(iPhi,iEta,towers[i].getPhi(),towers[i].getEta())<ISOL_Cone) {
1450 cout<<"delta R pour lequel on add "<<DeltaR(iPhi,iEta,towers[i].getPhi(),towers[i].getEta())<<endl;
1451 etiso += towers[i].getET(); // dR cut on tracks
1452 }
1453 }
1454
1455 return etiso;
1456}
1457
1458
1459
[321]1460// ******* Calorimetric isolation
[550]1461float RESOLution::CaloIsolation(const D_Particle& part, const D_CaloTowerList & towers, const float iPhi, const float iEta, const int iTower) {
1462 // iPhi and iEta are the coordinates of the center of the tower crossed by the particle
1463 // iPhi is in [-pi ; pi] and iEta is in [-etamax ; etamax]
1464 // iTower is the index of the tower, in [0, n_tower]. iTower points only towers in positive range
1465
[321]1466 if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
[392]1467 float et_sum=0;
1468
1469 unsigned int N = ISOL_Calo_Grid;
[550]1470 if(N%2 ==0 || N==0) { cout << "Warning: ISOL_Calo_Grid must be odd. ISOL_Calo_Grid = 3 will be used\n"; N=3;}
1471 int index= iTower; // index of the central tower of the grid in TOWER_eta_edges[.];
1472 // !! TOWER_eta_edges is only with eta>0
1473
[551]1474 //cout << "iEta= " << iEta << "\tiPhi= " << iPhi << endl;
[392]1475 if(index != iUNDEFINED) {
[550]1476 int range = (int) (N-1)/2; // the (int) conversion is needed, as N is "unsigned int"
1477 for(int i= -range; i<= range; i++) { // shift of index in eta: e.g.: i = -1, 0, 1 if N=3;
1478 // goal : to identify the center of each cell in the NxN grid.
1479 // the eta/phi coordinates of each center will be in (eta_ith_tower,phi_ith_tower)
[551]1480 unsigned int i_eta = (index+i>=0) ? index + i : (int) -1*(index + i +1); // i_eta is the index in [0; index_eta_max]
1481 if(i_eta>=TOWER_number) continue; // avoid going too far: TOWER_dphi[TOWER_number], TOWER_eta_edges[TOWER_number+1]
1482 //cout << "i_eta" << i_eta << "\t TOWER_number = " << TOWER_number << endl;
[392]1483
[550]1484 float dphi = TOWER_dphi[i_eta]*pi/180.; // in rad // size in phi of the cells for this eta
1485 float eta_ith_tower = (TOWER_eta_edges[i_eta]+TOWER_eta_edges[i_eta+1])/2.0;
1486 if(iEta<0) eta_ith_tower *= -1; // needed if the central tower is in negative region
[551]1487 if(index+i<0) eta_ith_tower *= -1; // needed if the "eta=0"-axis is crossed by the grid
[550]1488 //cout << "pour eta_ith_tower=" << eta_ith_tower << ", on va dphi = " << dphi << endl;
1489 // !!! at this point, eta_ith_tower is the value in eta of the center of the current calo cell
1490
1491 for(int j= -range; j<= range; j++) { // shift of index in phi. e.g.: j = -1, 0, 1 if N=3;
1492 float phi_ith_tower = iPhi + j*dphi; // in radian
1493 //cout << "eta_ith_tower= " << eta_ith_tower << "\tphi_ith_tower= " << phi_ith_tower << "\tj= " << j << "\tdphi=" << dphi << endl;
1494 if(phi_ith_tower > pi) phi_ith_tower -= 2.*pi;
1495 else if (phi_ith_tower < -pi) phi_ith_tower += 2.*pi;
1496
1497 D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
[551]1498 //cout << "eta_ith_tower= " << eta_ith_tower << "\tphi_ith_tower= " << phi_ith_tower <<"\t"
1499 // << "calMuon.getEta= " << calMuon.getEta() << "\tcalMuon.getPhi()= " << calMuon.getPhi() <<"\t";
[550]1500
[551]1501 if(calMuon.getEta() != UNDEFINED && calMuon.getET() > ISOL_Calo_ET) {
1502 et_sum += calMuon.getET();
1503 //cout << calMuon.getET() << " GeV";
[550]1504 }
[551]1505 //else cout << " not active";
1506 } // j-loop (phi indices)
1507 } // i-loop (eta indices)
[550]1508 } // undefined index
1509 else {
1510 cout << "out of range" << endl;
1511 if (CEN_max_mu < CEN_max_calo_fwd)
[392]1512 cout << "** ERROR in RESOLution::CaloIsolation: 'muon'-tower not found! **" << endl;
[551]1513 } // should never happen ! this would be a bug
[392]1514
[551]1515 return (part.Pt()==0)? UNDEFINED: et_sum/part.Pt();
[321]1516}
[31]1517
[321]1518
[71]1519 //********** returns a segmented value for eta and phi, for calo towers *****
[550]1520int RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
[423]1521 iEta = UNDEFINED;
1522 int index= iUNDEFINED;
1523 for (unsigned int i=1; i< TOWER_number+1; i++) {
1524 if(fabs(eta)>TOWER_eta_edges[i-1] && fabs(eta)<=TOWER_eta_edges[i]) {
1525 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);
1526 index = i-1;
1527 break;
1528 }
1529 }
1530 iPhi = UNDEFINED;
[550]1531 if(index==iUNDEFINED) return index;
1532
[423]1533 float dphi = TOWER_dphi[index]*pi/180.;
1534 for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
1535 float low = -pi+(i-1)*dphi;
1536 float high= low+dphi;
1537 if(phi > low && phi <= high ){
1538 iPhi = (low+high)/2.0;
1539 break;
1540 }
1541 }
1542 if (phi > pi-dphi) iPhi = pi-dphi;
[550]1543 return index;
[423]1544}
[71]1545
[2]1546//**************************** Returns the delta Phi ****************************
1547float DeltaPhi(const float phi1, const float phi2) {
[244]1548 float deltaphi=phi1-phi2; // in here, -pi < phi < pi
1549 if(fabs(deltaphi) > pi) {
1550 deltaphi=2.*pi -fabs(deltaphi);// put deltaphi between 0 and pi
[219]1551 }
[2]1552 else deltaphi=fabs(deltaphi);
1553
1554 return deltaphi;
1555}
1556
1557//**************************** Returns the delta R****************************
1558float DeltaR(const float phi1, const float eta1, const float phi2, const float eta2) {
1559 return sqrt(pow(DeltaPhi(phi1,phi2),2) + pow(eta1-eta2,2));
1560}
1561
1562int sign(const int myint) {
1563 if (myint >0) return 1;
1564 else if (myint <0) return -1;
1565 else return 0;
1566}
1567
1568int sign(const float myfloat) {
1569 if (myfloat >0) return 1;
1570 else if (myfloat <0) return -1;
1571 else return 0;
1572}
1573
[270]1574int ChargeVal(const int pid)
[55]1575{
[380]1576 cout << "ChargeVal :: deprecated function, do not use it anymore" << endl;
[55]1577 int charge;
1578 if(
1579 (pid == pGAMMA) ||
1580 (pid == pPI0) ||
1581 (pid == pK0L) ||
1582 (pid == pN) ||
1583 (pid == pSIGMA0) ||
1584 (pid == pDELTA0) ||
1585 (pid == pK0S) // not charged particles : invisible by tracker
1586 )
1587 charge = 0;
[376]1588 else charge = sign(pid);
[55]1589 return charge;
1590
[2]1591}
[380]1592
1593//------------------------------------------------------------------------------
1594void RESOLution::ReadParticleDataGroupTable() {
1595
1596 string temp_string;
1597 istringstream curstring;
1598
1599 ifstream fichier_a_lire(PdgTableFilename.c_str());
1600 if(!fichier_a_lire.good()) {
1601 cout <<"** ERROR: PDG Table ("<< PdgTableFilename
1602 << ") not found! exit. **" << endl;
1603 exit(1);
1604 return;
1605 }
1606 // first three lines of the file are useless
1607 getline(fichier_a_lire,temp_string);
1608 getline(fichier_a_lire,temp_string);
1609 getline(fichier_a_lire,temp_string);
1610
1611
1612 while (getline(fichier_a_lire,temp_string)) {
1613 curstring.clear(); // needed when using several times istringstream::str(string)
1614 curstring.str(temp_string);
[469]1615 long int ID; std::string name; int charge; float mass; float width; float lifetime;
[380]1616 // ID name chg mass total width lifetime
1617 // 1 d -1 0.33000 0.00000 0.00000E+00
[404]1618 // in the table, the charge is in units of e+/3
1619 // the total width is in GeV
1620 // the lifetime is ctau in mm
[380]1621 curstring >> ID >> name >> charge >> mass >> width >> lifetime;
[404]1622 PdgParticle particle(ID,name,mass,charge/3.,width,lifetime/1000.);
[380]1623 PDGtable.insert(ID,particle);
1624 //PdgTable.insert(pair<int,PdgParticle>(ID,particle));
1625 //cout << PDGtable[ID].name() << "\t" << PDGtable[ID].mass() << "\t" << PDGtable[ID].charge() << endl;
1626 }
1627
1628} // ReadParticleDataGroupTable
[399]1629
1630
1631// to be improved in order to avoid code repetition
1632// sorry, no time to do it right now (XR, 19/05/2009)
1633void 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) {
1634
1635TStopwatch globalwatch(global), loopwatch(loop), triggerwatch(trigger), frogwatch(frog), lhcowatch(lhco);
1636
1637 cout <<"** **"<< endl;
1638 cout <<"** ################## Time report ################# **"<< endl;
1639 cout << left << setw(32) <<"** Time report for "<<""
1640 << left << setw(15) << allEntries <<""
1641 << right << setw(22) <<"events **"<<endl;
1642 cout <<"** **"<< endl;
1643 cout << left << setw(10) <<"**"<<""
1644 << left << setw(15) <<"Time (s):"<<""
1645 << right << setw(15) <<"CPU"<<""
1646 << right << setw(15) <<"Real"<<""
1647 << right << setw(14) <<"**"<<endl;
1648 cout << left << setw(10) <<"**"<<""
1649 << left << setw(15) <<" + Global:"<<""
1650 << right << setw(15) <<globalwatch.CpuTime()<<""
1651 << right << setw(15) <<globalwatch.RealTime()<<""
1652 << right << setw(14) <<"**"<<endl;
1653 cout << left << setw(10) <<"**"<<""
1654 << left << setw(15) <<" + Events:"<<""
1655 << right << setw(15) <<loopwatch.CpuTime()<<""
1656 << right << setw(15) <<loopwatch.RealTime()<<""
1657 << right << setw(14) <<"**"<<endl;
1658 if(flag_trigger == 1)
1659 {
1660 cout << left << setw(10) <<"**"<<""
1661 << left << setw(15) <<" + Trigger:"<<""
1662 << right << setw(15) <<triggerwatch.CpuTime()<<""
1663 << right << setw(15) <<triggerwatch.RealTime()<<""
1664 << right << setw(14) <<"**"<<endl;
1665 }
1666 if(flag_frog == 1)
1667 {
1668 cout << left << setw(10) <<"**"<<""
1669 << left << setw(15) <<" + Frog:"<<""
1670 << right << setw(15) <<frogwatch.CpuTime()<<""
1671 << right << setw(15) <<frogwatch.RealTime()<<""
1672 << right << setw(14) <<"**"<<endl;
1673 }
1674 if(flag_lhco == 1)
1675 {
1676 cout << left << setw(10) <<"**"<<""
1677 << left << setw(15) <<" + LHCO:"<<""
1678 << right << setw(15) <<lhcowatch.CpuTime()<<""
1679 << right << setw(15) <<lhcowatch.RealTime()<<""
1680 << right << setw(14) <<"**"<<endl;
1681 }
1682
1683
1684 ofstream f_out(LogName.c_str(),ios_base::app);
1685
1686 f_out <<"** *"<< endl;
1687 f_out <<"** ################## Time report ################# *"<< endl;
1688 f_out << left << setw(32) <<"** Time report for "<<""
1689 << left << setw(15) << allEntries <<""
[472]1690 << right << setw(23) <<"events *"<<endl;
[399]1691 f_out <<"** *"<< endl;
1692 f_out << left << setw(10) <<"**"<<""
1693 << left << setw(15) <<"Time (s):"<<""
1694 << right << setw(15) <<"CPU"<<""
1695 << right << setw(15) <<"Real"<<""
1696 << right << setw(15) <<" *"<<endl;
1697 f_out << left << setw(10) <<"**"<<""
1698 << left << setw(15) <<" + Global:"<<""
1699 << right << setw(15) <<globalwatch.CpuTime()<<""
1700 << right << setw(15) <<globalwatch.RealTime()<<""
1701 << right << setw(15) <<" *"<<endl;
1702 f_out << left << setw(10) <<"**"<<""
1703 << left << setw(15) <<" + Events:"<<""
1704 << right << setw(15) <<loopwatch.CpuTime()<<""
1705 << right << setw(15) <<loopwatch.RealTime()<<""
1706 << right << setw(15) <<" *"<<endl;
1707 if(flag_trigger == 1)
1708 {
1709 f_out << left << setw(10) <<"**"<<""
1710 << left << setw(15) <<" + Trigger:"<<""
1711 << right << setw(15) <<triggerwatch.CpuTime()<<""
1712 << right << setw(15) <<triggerwatch.RealTime()<<""
1713 << right << setw(15) <<" *"<<endl;
1714 }
1715 if(flag_frog == 1)
1716 {
1717 f_out << left << setw(10) <<"**"<<""
1718 << left << setw(15) <<" + Frog:"<<""
1719 << right << setw(15) <<frogwatch.CpuTime()<<""
1720 << right << setw(15) <<frogwatch.RealTime()<<""
1721 << right << setw(15) <<" *"<<endl;
1722 }
1723 if(flag_lhco == 1)
1724 {
1725 f_out << left << setw(10) <<"**"<<""
1726 << left << setw(15) <<" + LHCO:"<<""
1727 << right << setw(15) <<lhcowatch.CpuTime()<<""
1728 << right << setw(15) <<lhcowatch.RealTime()<<""
1729 << right << setw(15) <<" *"<<endl;
1730 }
[472]1731 f_out << left << setw(16) << "** " << ""
1732 << left << setw(52) << get_time_date() << "**" << endl;
1733
1734
[399]1735 f_out<<"* *"<<"\n";
1736 f_out<<"* *"<<"\n";
1737 f_out<<"#....................................................................*"<<"\n";
1738 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
1739
1740 f_out.close();
1741
1742}
[404]1743
1744void print_header() {
1745 cout << endl << endl;
1746
1747 cout <<"*********************************************************************"<< endl;
1748 cout <<"*********************************************************************"<< endl;
1749 cout <<"** **"<< endl;
1750 cout <<"** Welcome to **"<< endl;
1751 cout <<"** **"<< endl;
1752 cout <<"** **"<< endl;
1753 cout <<"** .ddddddd- lL hH **"<< endl;
1754 cout <<"** -Dd` `dD: Ll hH` **"<< endl;
1755 cout <<"** dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss **"<< endl;
1756 cout <<"** -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs **"<< endl;
1757 cout <<"** dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. **"<< endl;
1758 cout <<"** .Dd :dd eE. LlL PpppPP Hh Hh eE sSS **"<< endl;
1759 cout <<"** dddddd:. eee+: lL. pp. hh. hh eee+ sssssS **"<< endl;
1760 cout <<"** Pp **"<< endl;
1761 cout <<"** **"<< endl;
1762 cout <<"** Delphes, a framework for the fast simulation **"<< endl;
1763 cout <<"** of a generic collider experiment **"<< endl;
1764 cout <<"** arXiv:0903.2225v1 [hep-ph] **"<< endl;
1765 cout <<"** **"<< endl;
[528]1766 cout <<"** --- Version 1.8 of Delphes --- **"<< endl;
1767 cout <<"** Last date of change: 16 August 2009 **"<< endl;
[404]1768 cout <<"** **"<< endl;
1769 cout <<"** **"<< endl;
1770 cout <<"** This package uses: **"<< endl;
1771 cout <<"** ------------------ **"<< endl;
1772 cout <<"** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
1773 cout <<"** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **"<< endl;
1774 cout <<"** FROG: [hep-ex/0901.2718v1] **"<< endl;
1775 cout <<"** **"<< endl;
1776 cout <<"**-----------------------------------------------------------------**"<< endl;
1777 cout <<"** **"<< endl;
1778 cout <<"** Main authors: **"<< endl;
1779 cout <<"** ------------- **"<< endl;
1780 cout <<"** **"<< endl;
1781 cout <<"** Séverine Ovyn Xavier Rouby **"<< endl;
1782 cout <<"** severine.ovyn@uclouvain.be xavier.rouby@cern **"<< endl;
1783 cout <<"** Center for Particle Physics and Phenomenology (CP3) **"<< endl;
1784 cout <<"** Universite Catholique de Louvain (UCL) **"<< endl;
1785 cout <<"** Louvain-la-Neuve, Belgium **"<< endl;
1786 cout <<"** **"<< endl;
1787 cout <<"**-----------------------------------------------------------------**"<< endl;
1788 cout <<"** **"<< endl;
1789 cout <<"** Former Delphes versions and documentation can be found on : **"<< endl;
1790 cout <<"** http://www.fynu.ucl.ac.be/delphes.html **"<< endl;
1791 cout <<"** **"<< endl;
1792 cout <<"** **"<< endl;
[528]1793 cout <<"** Disclaimer: this program comes without guarantees. **"<< endl;
1794 cout <<"** Beware of errors and please give us **"<< endl;
1795 cout <<"** your feedbacks about potential bugs. **"<< endl;
[404]1796 cout <<"** **"<< endl;
1797 cout <<"*********************************************************************"<< endl;
1798 cout <<"*********************************************************************"<< endl;
1799
1800}
[465]1801
1802string get_time_date() {
1803 time_t rawtime;
1804 struct tm * timeinfo;
1805
1806 time ( &rawtime );
1807 timeinfo = localtime ( &rawtime );
1808
1809 char temp[100];
[553]1810 //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);
1811 sprintf(temp,"%i",timeinfo->tm_mday);
1812 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);
1813 sprintf(temp,"%s/%d %d",temp,timeinfo->tm_year+1900,timeinfo->tm_hour);
1814 if(timeinfo->tm_min<10) sprintf(temp,"%s:0%d",temp,timeinfo->tm_min); else sprintf(temp,"%s:%d",temp,timeinfo->tm_min);
1815 if(timeinfo->tm_sec<10) sprintf(temp,"%s:0%d",temp,timeinfo->tm_sec); else sprintf(temp,"%s:%d",temp,timeinfo->tm_sec);
[465]1816 string tempstring(temp);
1817 return tempstring;
1818}
1819
[494]1820void warning(const string oldname, const string newname) {
1821 string text = "** Warning in datacard: " + oldname + " deprecated. Use " + newname + " instead.";
1822 cout << left << setw(67) << text
1823 << right << setw(2) <<"**" << endl;
1824}
Note: See TracBrowser for help on using the repository browser.