Fork me on GitHub

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

Last change on this file since 377 was 376, checked in by Xavier Rouby, 16 years ago

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