Fork me on GitHub

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

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

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