Fork me on GitHub

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

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