Fork me on GitHub

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

Last change on this file since 665 was 665, checked in by cp3-support, 12 years ago

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