Fork me on GitHub

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

Last change on this file since 172 was 94, checked in by severine ovyn, 16 years ago

Add frog plus cleaning + bugs remove
File size: 44.6 KB
Line 
1/*
2 ---- Delphes ----
3 A Fast Simulator for general purpose LHC detector
4 S. Ovyn ~~~~ severine.ovyn@uclouvain.be
5
6 Center for Particle Physics and Phenomenology (CP3)
7 Universite Catholique de Louvain (UCL)
8 Louvain-la-Neuve, Belgium
9*/
10
11/// \file SmearUtil.cc
12/// \brief RESOLution class, and some generic definitions
13
14
15#include "interface/SmearUtil.h"
16#include "TRandom.h"
17
18#include <iostream>
19#include <sstream>
20#include <fstream>
21#include <iomanip>
22
23
24
25using namespace std;
26
27//------------------------------------------------------------------------------
28
29RESOLution::RESOLution() {
30
31 // Detector characteristics
32 CEN_max_tracker = 2.5; // Maximum tracker coverage
33 CEN_max_calo_cen = 3.0; // central calorimeter coverage
34 CEN_max_calo_fwd = 5.0; // forward calorimeter pseudorapidity coverage
35 CEN_max_mu = 2.4; // muon chambers pseudorapidity coverage
36
37 // Energy resolution for electron/photon
38 // \sigma/E = C + N/E + S/\sqrt{E}
39 ELG_Scen = 0.05; // S term for central ECAL
40 ELG_Ncen = 0.25; // N term for central ECAL
41 ELG_Ccen = 0.005; // C term for central ECAL
42 ELG_Cfwd = 0.107; // S term for forward ECAL
43 ELG_Sfwd = 2.084; // C term for forward ECAL
44 ELG_Nfwd = 0.0; // N term for central ECAL
45
46 // Energy resolution for hadrons in ecal/hcal/hf
47 // \sigma/E = C + N/E + S/\sqrt{E}
48 HAD_Shcal = 1.5; // S term for central HCAL // hadronic calorimeter
49 HAD_Nhcal = 0.; // N term for central HCAL
50 HAD_Chcal = 0.05; // C term for central HCAL
51 HAD_Shf = 2.7; // S term for HF // forward calorimeter
52 HAD_Nhf = 0.; // N term for HF
53 HAD_Chf = 0.13; // C term for HF
54
55 // Muon smearing
56 MU_SmearPt = 0.01;
57
58 // Tracking efficiencies
59 TRACK_ptmin = 0.9; // minimal pt needed to reach the calorimeter in GeV
60 TRACK_eff = 100; // efficiency associated to the tracking
61
62 // Calorimetric towers
63 TOWER_number = 40;
64 const float tower_eta_edges[41] = {
65 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,
66 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,
67 4.350, 4.525, 4.700, 5.000}; // temporary object
68 TOWER_eta_edges = new float[TOWER_number+1];
69 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = tower_eta_edges[i];
70
71 const float tower_dphi[40] = {
72 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 10,
73 10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; // temporary object
74 TOWER_dphi = new float[TOWER_number];
75 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = tower_dphi[i];
76
77
78 // Thresholds for reconstructed objetcs
79 PTCUT_elec = 10.0;
80 PTCUT_muon = 10.0;
81 PTCUT_jet = 20.0;
82 PTCUT_gamma = 10.0;
83 PTCUT_taujet = 10.0;
84
85 // General jet variable
86 JET_coneradius = 0.7; // generic jet radius ; not for tau's !!!
87 JET_jetalgo = 1; // 1 for Cone algorithm, 2 for MidPoint algorithm, 3 for SIScone algorithm, 4 for kt algorithm
88 JET_seed = 1.0; // minimum seed to start jet reconstruction
89
90 // Tagging definition
91 BTAG_b = 40;
92 BTAG_mistag_c = 10;
93 BTAG_mistag_l = 1;
94
95 // FLAGS
96 FLAG_bfield = 1; //1 to run the bfield propagation else 0
97 FLAG_vfd = 1; //1 to run the very forward detectors else 0
98 FLAG_trigger = 1; //1 to run the trigger selection else 0
99 FLAG_frog = 1; //1 to run the FROG event display
100
101 // In case BField propagation allowed
102 TRACK_radius = 129; //radius of the BField coverage
103 TRACK_length = 300; //length of the BField coverage
104 TRACK_bfield_x = 0; //X composant of the BField
105 TRACK_bfield_y = 0; //Y composant of the BField
106 TRACK_bfield_z = 3.8; //Z composant of the BField
107
108 // In case Very forward detectors allowed
109 VFD_min_calo_vfd = 5.2; // very forward calorimeter (if any) like CASTOR
110 VFD_max_calo_vfd = 6.6;
111 VFD_min_zdc = 8.3;
112 VFD_s_zdc = 140; // distance of the Zero Degree Calorimeter, from the Interaction poin, in [m]
113
114 RP_220_s = 220; // distance of the RP to the IP, in meters
115 RP_220_x = 0.002; // distance of the RP to the beam, in meters
116 RP_420_s = 420; // distance of the RP to the IP, in meters
117 RP_420_x = 0.004; // distance of the RP to the beam, in meters
118
119 // In case FROG event display allowed
120 NEvents_Frog = 10;
121
122 //********************************************
123 //jet stuffs not defined in the input datacard
124 //********************************************
125
126 JET_overlap = 0.75;
127 // MidPoint algorithm definition
128 JET_M_coneareafraction = 0.25;
129 JET_M_maxpairsize = 2;
130 JET_M_maxiterations = 100;
131 // Define Cone algorithm.
132 JET_C_adjacencycut = 2;
133 JET_C_maxiterations = 100;
134 JET_C_iratch = 1;
135 //Define SISCone algorithm.
136 JET_S_npass = 0;
137 JET_S_protojet_ptmin= 0.0;
138
139 //For Tau-jet definition
140 TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
141 TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
142 TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
143 TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
144
145 PT_QUARKS_MIN = 2.0 ; // minimal pt needed by quarks to do b-tag
146
147}
148
149//------------------------------------------------------------------------------
150void RESOLution::ReadDataCard(const string datacard) {
151
152 string temp_string;
153 istringstream curstring;
154
155 ifstream fichier_a_lire(datacard.c_str());
156 if(!fichier_a_lire.good()) {
157 cout << datacard << "Datadard " << datacard << " not found, use default values" << endl;
158 return;
159 }
160
161 while (getline(fichier_a_lire,temp_string)) {
162 curstring.clear(); // needed when using several times istringstream::str(string)
163 curstring.str(temp_string);
164 string varname;
165 float value; int ivalue;
166
167 if(strstr(temp_string.c_str(),"#")) { }
168 else if(strstr(temp_string.c_str(),"CEN_max_tracker")) {curstring >> varname >> value; CEN_max_tracker = value;}
169 else if(strstr(temp_string.c_str(),"CEN_max_calo_cen")) {curstring >> varname >> value; CEN_max_calo_cen = value;}
170 else if(strstr(temp_string.c_str(),"CEN_max_calo_fwd")) {curstring >> varname >> value; CEN_max_calo_fwd = value;}
171 else if(strstr(temp_string.c_str(),"CEN_max_mu")) {curstring >> varname >> value; CEN_max_mu = value;}
172
173 else if(strstr(temp_string.c_str(),"VFD_min_calo_vfd")) {curstring >> varname >> value; VFD_min_calo_vfd = value;}
174 else if(strstr(temp_string.c_str(),"VFD_max_calo_vfd")) {curstring >> varname >> value; VFD_max_calo_vfd = value;}
175 else if(strstr(temp_string.c_str(),"VFD_min_zdc")) {curstring >> varname >> value; VFD_min_zdc = value;}
176 else if(strstr(temp_string.c_str(),"VFD_s_zdc")) {curstring >> varname >> value; VFD_s_zdc = value;}
177
178 else if(strstr(temp_string.c_str(),"RP_220_s")) {curstring >> varname >> value; RP_220_s = value;}
179 else if(strstr(temp_string.c_str(),"RP_220_x")) {curstring >> varname >> value; RP_220_x = value;}
180 else if(strstr(temp_string.c_str(),"RP_420_s")) {curstring >> varname >> value; RP_420_s = value;}
181 else if(strstr(temp_string.c_str(),"RP_420_x")) {curstring >> varname >> value; RP_420_x = value;}
182
183 else if(strstr(temp_string.c_str(),"ELG_Scen")) {curstring >> varname >> value; ELG_Scen = value;}
184 else if(strstr(temp_string.c_str(),"ELG_Ncen")) {curstring >> varname >> value; ELG_Ncen = value;}
185 else if(strstr(temp_string.c_str(),"ELG_Ccen")) {curstring >> varname >> value; ELG_Ccen = value;}
186 else if(strstr(temp_string.c_str(),"ELG_Sfwd")) {curstring >> varname >> value; ELG_Sfwd = value;}
187 else if(strstr(temp_string.c_str(),"ELG_Cfwd")) {curstring >> varname >> value; ELG_Cfwd = value;}
188 else if(strstr(temp_string.c_str(),"ELG_Nfwd")) {curstring >> varname >> value; ELG_Nfwd = value;}
189 else if(strstr(temp_string.c_str(),"HAD_Shcal")) {curstring >> varname >> value; HAD_Shcal = value;}
190 else if(strstr(temp_string.c_str(),"HAD_Nhcal")) {curstring >> varname >> value; HAD_Nhcal = value;}
191 else if(strstr(temp_string.c_str(),"HAD_Chcal")) {curstring >> varname >> value; HAD_Chcal = value;}
192 else if(strstr(temp_string.c_str(),"HAD_Shf")) {curstring >> varname >> value; HAD_Shf = value;}
193 else if(strstr(temp_string.c_str(),"HAD_Nhf")) {curstring >> varname >> value; HAD_Nhf = value;}
194 else if(strstr(temp_string.c_str(),"HAD_Chf")) {curstring >> varname >> value; HAD_Chf = value;}
195 else if(strstr(temp_string.c_str(),"MU_SmearPt")) {curstring >> varname >> value; MU_SmearPt = value;}
196
197 else if(strstr(temp_string.c_str(),"TRACK_radius")) {curstring >> varname >> ivalue;TRACK_radius = ivalue;}
198 else if(strstr(temp_string.c_str(),"TRACK_length")) {curstring >> varname >> ivalue;TRACK_length = ivalue;}
199 else if(strstr(temp_string.c_str(),"TRACK_bfield_x")) {curstring >> varname >> value; TRACK_bfield_x = value;}
200 else if(strstr(temp_string.c_str(),"TRACK_bfield_y")) {curstring >> varname >> value; TRACK_bfield_y = value;}
201 else if(strstr(temp_string.c_str(),"TRACK_bfield_z")) {curstring >> varname >> value; TRACK_bfield_z = value;}
202 else if(strstr(temp_string.c_str(),"FLAG_bfield")) {curstring >> varname >> ivalue; FLAG_bfield = ivalue;}
203 else if(strstr(temp_string.c_str(),"TRACK_ptmin")) {curstring >> varname >> value; TRACK_ptmin = value;}
204 else if(strstr(temp_string.c_str(),"TRACK_eff")) {curstring >> varname >> ivalue;TRACK_eff = ivalue;}
205
206 else if(strstr(temp_string.c_str(),"TOWER_number")) {curstring >> varname >> ivalue;TOWER_number = ivalue;}
207 else if(strstr(temp_string.c_str(),"TOWER_eta_edges")){
208 curstring >> varname; for(unsigned int i=0; i<TOWER_number+1; i++) {curstring >> value; TOWER_eta_edges[i] = value;} }
209 else if(strstr(temp_string.c_str(),"TOWER_dphi")){
210 curstring >> varname; for(unsigned int i=0; i<TOWER_number; i++) {curstring >> value; TOWER_dphi[i] = value;} }
211
212 else if(strstr(temp_string.c_str(),"PTCUT_elec")) {curstring >> varname >> value; PTCUT_elec = value;}
213 else if(strstr(temp_string.c_str(),"PTCUT_muon")) {curstring >> varname >> value; PTCUT_muon = value;}
214 else if(strstr(temp_string.c_str(),"PTCUT_jet")) {curstring >> varname >> value; PTCUT_jet = value;}
215 else if(strstr(temp_string.c_str(),"PTCUT_gamma")) {curstring >> varname >> value; PTCUT_gamma = value;}
216 else if(strstr(temp_string.c_str(),"PTCUT_taujet")) {curstring >> varname >> value; PTCUT_taujet = value;}
217
218 else if(strstr(temp_string.c_str(),"JET_coneradius")) {curstring >> varname >> value; JET_coneradius = value;}
219 else if(strstr(temp_string.c_str(),"JET_jetalgo")) {curstring >> varname >> ivalue;JET_jetalgo = ivalue;}
220 else if(strstr(temp_string.c_str(),"JET_seed")) {curstring >> varname >> value; JET_seed = value;}
221
222 else if(strstr(temp_string.c_str(),"BTAG_b")) {curstring >> varname >> ivalue;BTAG_b = ivalue;}
223 else if(strstr(temp_string.c_str(),"BTAG_mistag_c")) {curstring >> varname >> ivalue;BTAG_mistag_c = ivalue;}
224 else if(strstr(temp_string.c_str(),"BTAG_mistag_l")) {curstring >> varname >> ivalue;BTAG_mistag_l = ivalue;}
225
226 else if(strstr(temp_string.c_str(),"FLAG_vfd")) {curstring >> varname >> ivalue; FLAG_vfd = ivalue;}
227 else if(strstr(temp_string.c_str(),"FLAG_trigger")) {curstring >> varname >> ivalue; FLAG_trigger = ivalue;}
228 else if(strstr(temp_string.c_str(),"FLAG_frog")) {curstring >> varname >> ivalue; FLAG_frog = ivalue;}
229 else if(strstr(temp_string.c_str(),"NEvents_Frog")) {curstring >> varname >> ivalue; NEvents_Frog = ivalue;}
230 }
231
232 //jet stuffs not defined in the input datacard
233 JET_overlap = 0.75;
234 // MidPoint algorithm definition
235 JET_M_coneareafraction = 0.25;
236 JET_M_maxpairsize = 2;
237 JET_M_maxiterations = 100;
238 // Define Cone algorithm.
239 JET_C_adjacencycut = 2;
240 JET_C_maxiterations = 100;
241 JET_C_iratch = 1;
242 //Define SISCone algorithm.
243 JET_S_npass = 0;
244 JET_S_protojet_ptmin= 0.0;
245
246 //For Tau-jet definition
247 TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
248 TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
249 TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
250 TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
251
252}
253
254void RESOLution::Logfile(string LogName) {
255 //void RESOLution::Logfile(string outputfilename) {
256
257 ofstream f_out(LogName.c_str());
258
259 f_out<<"#*********************************************************************"<<"\n";
260 f_out<<"# *"<<"\n";
261 f_out<<"# ---- DELPHES release 1.0 ---- *"<<"\n";
262 f_out<<"# *"<<"\n";
263 f_out<<"# A Fast Simulator for general purpose LHC detector *"<<"\n";
264 f_out<<"# Written by S. Ovyn and X. Rouby *"<<"\n";
265 f_out<<"# severine.ovyn@uclouvain.be *"<<"\n";
266 f_out<<"# *"<<"\n";
267 f_out<<"# http: *"<<"\n";
268 f_out<<"# *"<<"\n";
269 f_out<<"# Center for Particle Physics and Phenomenology (CP3) *"<<"\n";
270 f_out<<"# Universite Catholique de Louvain (UCL) *"<<"\n";
271 f_out<<"# Louvain-la-Neuve, Belgium *"<<"\n";
272 f_out<<"# *"<<"\n";
273 f_out<<"#....................................................................*"<<"\n";
274 f_out<<"# *"<<"\n";
275 f_out<<"# This package uses: *"<<"\n";
276 f_out<<"# FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] *"<<"\n";
277 f_out<<"# Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] *"<<"\n";
278 f_out<<"# ExRootAnalysis *"<<"\n";
279 f_out<<"# *"<<"\n";
280 f_out<<"#....................................................................*"<<"\n";
281 f_out<<"# *"<<"\n";
282 f_out<<"# This file contains all the running parameters (detector and cuts) *"<<"\n";
283 f_out<<"# necessary to reproduce the detector simulation *"<<"\n";
284 f_out<<"# *"<<"\n";
285 f_out<<"#....................................................................*"<<"\n";
286 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
287 f_out<<"* *"<<"\n";
288 f_out<<"#******************************** *"<<"\n";
289 f_out<<"# Central detector caracteristics *"<<"\n";
290 f_out<<"#******************************** *"<<"\n";
291 f_out<<"* *"<<"\n";
292 f_out << left << setw(30) <<"* Maximum tracking system: "<<""
293 << left << setw(10) <<CEN_max_tracker <<""<< right << setw(15)<<"*"<<"\n";
294 f_out << left << setw(30) <<"* Maximum central calorimeter: "<<""
295 << left << setw(10) <<CEN_max_calo_cen <<""<< right << setw(15)<<"*"<<"\n";
296 f_out << left << setw(30) <<"* Maximum forward calorimeter: "<<""
297 << left << setw(10) <<CEN_max_calo_fwd <<""<< right << setw(15)<<"*"<<"\n";
298 f_out << left << setw(30) <<"* Muon chambers coverage: "<<""
299 << left << setw(10) <<CEN_max_mu <<""<< right << setw(15)<<"*"<<"\n";
300 f_out<<"* *"<<"\n";
301 if(FLAG_vfd==1){
302 f_out<<"#********************************** *"<<"\n";
303 f_out<<"# Very forward detector switches on *"<<"\n";
304 f_out<<"#********************************** *"<<"\n";
305 f_out<<"* *"<<"\n";
306 f_out << left << setw(55) <<"* Minimum very forward calorimeter: "<<""
307 << left << setw(5) <<VFD_min_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
308 f_out << left << setw(55) <<"* Maximum very forward calorimeter: "<<""
309 << left << setw(5) <<VFD_max_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
310 f_out << left << setw(55) <<"* Minimum coverage zero_degree calorimeter "<<""
311 << left << setw(5) <<VFD_min_zdc <<""<< right << setw(10)<<"*"<<"\n";
312 f_out << left << setw(55) <<"* Distance of the ZDC to the IP, in meters: "<<""
313 << left << setw(5) <<VFD_s_zdc <<""<< right << setw(10)<<"*"<<"\n";
314 f_out << left << setw(55) <<"* Distance of the RP to the IP, in meters: "<<""
315 << left << setw(5) <<RP_220_s <<""<< right << setw(10)<<"*"<<"\n";
316 f_out << left << setw(55) <<"* Distance of the RP to the beam, in meters: "<<""
317 << left << setw(5) <<RP_220_x <<""<< right << setw(10)<<"*"<<"\n";
318 f_out << left << setw(55) <<"* Distance of the RP to the IP, in meters: "<<""
319 << left << setw(5) <<RP_420_s <<""<< right << setw(10)<<"*"<<"\n";
320 f_out << left << setw(55) <<"* Distance of the RP to the beam, in meters: "<<""
321 << left << setw(5) <<RP_420_x <<""<< right << setw(10)<<"*"<<"\n";
322 f_out<<"* *"<<"\n";
323 }
324 else {
325 f_out<<"#*********************************** *"<<"\n";
326 f_out<<"# Very forward detector switches off *"<<"\n";
327 f_out<<"#*********************************** *"<<"\n";
328 f_out<<"* *"<<"\n";
329 }
330 f_out<<"#************************************ *"<<"\n";
331 f_out<<"# Electromagnetic smearing parameters *"<<"\n";
332 f_out<<"#************************************ *"<<"\n";
333 f_out<<"* *"<<"\n";
334 //# \sigma/E = C + N/E + S/\sqrt{E}
335 f_out << left << setw(30) <<"* S term for central ECAL: "<<""
336 << left << setw(30) <<ELG_Scen <<""<< right << setw(10)<<"*"<<"\n";
337 f_out << left << setw(30) <<"* N term for central ECAL: "<<""
338 << left << setw(30) <<ELG_Ncen <<""<< right << setw(10)<<"*"<<"\n";
339 f_out << left << setw(30) <<"* C term for central ECAL: "<<""
340 << left << setw(30) <<ELG_Ccen <<""<< right << setw(10)<<"*"<<"\n";
341 f_out << left << setw(30) <<"* S term for forward ECAL: "<<""
342 << left << setw(30) <<ELG_Sfwd <<""<< right << setw(10)<<"*"<<"\n";
343 f_out << left << setw(30) <<"* N term for forward ECAL: "<<""
344 << left << setw(30) <<ELG_Nfwd <<""<< right << setw(10)<<"*"<<"\n";
345 f_out << left << setw(30) <<"* C term for forward ECAL: "<<""
346 << left << setw(30) <<ELG_Cfwd <<""<< right << setw(10)<<"*"<<"\n";
347 f_out<<"* *"<<"\n";
348 f_out<<"#***************************** *"<<"\n";
349 f_out<<"# Hadronic smearing parameters *"<<"\n";
350 f_out<<"#***************************** *"<<"\n";
351 f_out<<"* *"<<"\n";
352 f_out << left << setw(30) <<"* S term for central HCAL: "<<""
353 << left << setw(30) <<HAD_Shcal <<""<< right << setw(10)<<"*"<<"\n";
354 f_out << left << setw(30) <<"* N term for central HCAL: "<<""
355 << left << setw(30) <<HAD_Nhcal <<""<< right << setw(10)<<"*"<<"\n";
356 f_out << left << setw(30) <<"* C term for central HCAL: "<<""
357 << left << setw(30) <<HAD_Chcal <<""<< right << setw(10)<<"*"<<"\n";
358 f_out << left << setw(30) <<"* S term for forward HCAL: "<<""
359 << left << setw(30) <<HAD_Shf <<""<< right << setw(10)<<"*"<<"\n";
360 f_out << left << setw(30) <<"* N term for forward HCAL: "<<""
361 << left << setw(30) <<HAD_Nhf <<""<< right << setw(10)<<"*"<<"\n";
362 f_out << left << setw(30) <<"* C term for forward HCAL: "<<""
363 << left << setw(30) <<HAD_Chf <<""<< right << setw(10)<<"*"<<"\n";
364 f_out<<"* *"<<"\n";
365 f_out<<"#************************* *"<<"\n";
366 f_out<<"# Muon smearing parameters *"<<"\n";
367 f_out<<"#************************* *"<<"\n";
368 f_out<<"* *"<<"\n";
369 f_out << left << setw(55) <<"* PT resolution for muons : "<<""
370 << left << setw(5) <<MU_SmearPt <<""<< right << setw(10)<<"*"<<"\n";
371 f_out<<"* *"<<"\n";
372 if(FLAG_bfield==1){
373 f_out<<"#*************************** *"<<"\n";
374 f_out<<"# Magnetic field switches on *"<<"\n";
375 f_out<<"#*************************** *"<<"\n";
376 f_out<<"* *"<<"\n";
377 f_out << left << setw(55) <<"* Radius of the BField coverage: "<<""
378 << left << setw(5) <<TRACK_radius <<""<< right << setw(10)<<"*"<<"\n";
379 f_out << left << setw(55) <<"* Length of the BField coverage: "<<""
380 << left << setw(5) <<TRACK_length <<""<< right << setw(10)<<"*"<<"\n";
381 f_out << left << setw(55) <<"* BField X component: "<<""
382 << left << setw(5) <<TRACK_bfield_x <<""<< right << setw(10)<<"*"<<"\n";
383 f_out << left << setw(55) <<"* BField Y component: "<<""
384 << left << setw(5) <<TRACK_bfield_y <<""<< right << setw(10)<<"*"<<"\n";
385 f_out << left << setw(55) <<"* BField Z component: "<<""
386 << left << setw(5) <<TRACK_bfield_z <<""<< right << setw(10)<<"*"<<"\n";
387 f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
388 << left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
389 f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
390 << left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
391 f_out<<"* *"<<"\n";
392 }
393 else {
394 f_out<<"#**************************** *"<<"\n";
395 f_out<<"# Magnetic field switches off *"<<"\n";
396 f_out<<"#**************************** *"<<"\n";
397 f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
398 << left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
399 f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
400 << left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
401 f_out<<"* *"<<"\n";
402 }
403 f_out<<"#******************** *"<<"\n";
404 f_out<<"# Calorimetric Towers *"<<"\n";
405 f_out<<"#******************** *"<<"\n";
406 f_out << left << setw(55) <<"* Number of calorimetric towers in eta, for eta>0: "<<""
407 << left << setw(5) << TOWER_number <<""<< right << setw(10)<<"*"<<"\n";
408 f_out << left << setw(55) <<"* Tower edges in eta, for eta>0: "<<"" << right << setw(15)<<"*"<<"\n";
409 f_out << "* ";
410 for (unsigned int i=0; i<TOWER_number+1; i++) {
411 f_out << left << setw(7) << TOWER_eta_edges[i];
412 if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "* ";
413 }
414 for (unsigned int i=(TOWER_number+1)%9; i<9; i++) f_out << left << setw(7) << "";
415 f_out << right << setw(3)<<"*"<<"\n";
416 f_out << left << setw(55) <<"* Tower sizes in phi, for eta>0 [degree]:"<<"" << right << setw(15)<<"*"<<"\n";
417 f_out << "* ";
418 for (unsigned int i=0; i<TOWER_number; i++) {
419 f_out << left << setw(7) << TOWER_dphi[i];
420 if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "* ";
421 }
422 for (unsigned int i=(TOWER_number)%9; i<9; i++) f_out << left << setw(7) << "";
423 f_out << right << setw(3)<<"*"<<"\n";
424 f_out<<"* *"<<"\n";
425 f_out<<"#******************* *"<<"\n";
426 f_out<<"# Minimum pT's [GeV] *"<<"\n";
427 f_out<<"#******************* *"<<"\n";
428 f_out<<"* *"<<"\n";
429 f_out << left << setw(40) <<"* Minimum pT for electrons: "<<""
430 << left << setw(20) <<PTCUT_elec <<""<< right << setw(10)<<"*"<<"\n";
431 f_out << left << setw(40) <<"* Minimum pT for muons: "<<""
432 << left << setw(20) <<PTCUT_muon <<""<< right << setw(10)<<"*"<<"\n";
433 f_out << left << setw(40) <<"* Minimum pT for jets: "<<""
434 << left << setw(20) <<PTCUT_jet <<""<< right << setw(10)<<"*"<<"\n";
435 f_out << left << setw(40) <<"* Minimum pT for Tau-jets: "<<""
436 << left << setw(20) <<PTCUT_taujet <<""<< right << setw(10)<<"*"<<"\n";
437 f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
438 << left << setw(20) <<PTCUT_gamma <<""<< right << setw(10)<<"*"<<"\n";
439 f_out<<"* *"<<"\n";
440 f_out<<"#*************** *"<<"\n";
441 f_out<<"# Jet definition *"<<"\n";
442 f_out<<"#*************** *"<<"\n";
443 f_out<<"* *"<<"\n";
444 f_out<<"* Six algorithms are currently available: *"<<"\n";
445 f_out<<"* - 1) CDF cone algorithm, *"<<"\n";
446 f_out<<"* - 2) CDF MidPoint algorithm, *"<<"\n";
447 f_out<<"* - 3) SIScone algorithm, *"<<"\n";
448 f_out<<"* - 4) kt algorithm, *"<<"\n";
449 f_out<<"* - 5) Cambrigde/Aachen algorithm, *"<<"\n";
450 f_out<<"* - 6) Anti-kt algorithm. *"<<"\n";
451 f_out<<"* *"<<"\n";
452 f_out<<"* You have chosen *"<<"\n";
453 switch(JET_jetalgo) {
454 default:
455 case 1: {
456 f_out<<"* CDF JetClu jet algorithm with parameters: *"<<"\n";
457 f_out << left << setw(40) <<"* - Seed threshold: "<<""
458 << left << setw(10) <<JET_seed <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
459 f_out << left << setw(40) <<"* - Cone radius: "<<""
460 << left << setw(10) <<JET_coneradius <<""<< right << setw(20)<<"*"<<"\n";
461 f_out << left << setw(40) <<"* - Adjacency cut: "<<""
462 << left << setw(10) <<JET_C_adjacencycut <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
463 f_out << left << setw(40) <<"* - Max iterations: "<<""
464 << left << setw(10) <<JET_C_maxiterations <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
465 f_out << left << setw(40) <<"* - Iratch: "<<""
466 << left << setw(10) <<JET_C_iratch <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
467 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
468 << left << setw(10) <<JET_overlap <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
469 }
470 break;
471 case 2: {
472 f_out<<"* CDF midpoint jet algorithm with parameters: *"<<"\n";
473 f_out << left << setw(40) <<"* - Seed threshold: "<<""
474 << left << setw(20) <<JET_seed <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
475 f_out << left << setw(40) <<"* - Cone radius: "<<""
476 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
477 f_out << left << setw(40) <<"* - Cone area fraction:"<<""
478 << left << setw(20) <<JET_M_coneareafraction <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
479 f_out << left << setw(40) <<"* - Maximum pair size: "<<""
480 << left << setw(20) <<JET_M_maxpairsize <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
481 f_out << left << setw(40) <<"* - Max iterations: "<<""
482 << left << setw(20) <<JET_M_maxiterations <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
483 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
484 << left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
485 }
486 break;
487 case 3: {
488 f_out <<"* SISCone jet algorithm with parameters: *"<<"\n";
489 f_out << left << setw(40) <<"* - Cone radius: "<<""
490 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
491 f_out << left << setw(40) <<"* - Overlap threshold: "<<""
492 << left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
493 f_out << left << setw(40) <<"* - Number pass max: "<<""
494 << left << setw(20) <<JET_S_npass <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
495 f_out << left << setw(40) <<"* - Minimum pT for protojet: "<<""
496 << left << setw(20) <<JET_S_protojet_ptmin <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
497 }
498 break;
499 case 4: {
500 f_out <<"* KT jet algorithm with parameters: *"<<"\n";
501 f_out << left << setw(40) <<"* - Cone radius: "<<""
502 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
503 }
504 break;
505 case 5: {
506 f_out <<"* Cambridge/Aachen jet algorithm with parameters: *"<<"\n";
507 f_out << left << setw(40) <<"* - Cone radius: "<<""
508 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
509 }
510 break;
511 case 6: {
512 f_out <<"* Anti-kt jet algorithm with parameters: *"<<"\n";
513 f_out << left << setw(40) <<"* - Cone radius: "<<""
514 << left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
515 }
516 break;
517 }
518 f_out<<"* *"<<"\n";
519 f_out<<"#****************************** *"<<"\n";
520 f_out<<"# Tau-jet definition parameters *"<<"\n";
521 f_out<<"#****************************** *"<<"\n";
522 f_out<<"* *"<<"\n";
523 f_out << left << setw(45) <<"* Cone radius for calorimeter tagging: "<<""
524 << left << setw(5) <<TAU_energy_scone <<""<< right << setw(20)<<"*"<<"\n";
525 f_out << left << setw(45) <<"* Fraction of energy in the small cone: "<<""
526 << left << setw(5) <<TAU_energy_frac*100 <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
527 f_out << left << setw(45) <<"* Cone radius for tracking tagging: "<<""
528 << left << setw(5) <<TAU_track_scone <<""<< right << setw(20)<<"*"<<"\n";
529 f_out << left << setw(45) <<"* Minimum track pT [GeV]: "<<""
530 << left << setw(5) <<TAU_track_pt <<""<< right << setw(20)<<"*"<<"\n";
531 f_out<<"* *"<<"\n";
532 f_out<<"#*************************** *"<<"\n";
533 f_out<<"# B-tagging efficiencies [%] *"<<"\n";
534 f_out<<"#*************************** *"<<"\n";
535 f_out<<"* *"<<"\n";
536 f_out << left << setw(50) <<"* Efficiency to tag a \"b\" as a b-jet: "<<""
537 << left << setw(10) <<BTAG_b <<""<< right << setw(10)<<"*"<<"\n";
538 f_out << left << setw(50) <<"* Efficiency to mistag a c-jet as a b-jet: "<<""
539 << left << setw(10) <<BTAG_mistag_c <<""<< right << setw(10)<<"*"<<"\n";
540 f_out << left << setw(50) <<"* Efficiency to mistag a light jet as a b-jet: "<<""
541 << left << setw(10) <<BTAG_mistag_l <<""<< right << setw(10)<<"*"<<"\n";
542 f_out<<"* *"<<"\n";
543 f_out<<"* *"<<"\n";
544 f_out<<"#....................................................................*"<<"\n";
545 f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
546
547}
548
549// **********Provides the smeared TLorentzVector for the electrons********
550// Smears the electron energy, and changes the 4-momentum accordingly
551// different smearing if the electron is central (eta < 2.5) or forward
552void RESOLution::SmearElectron(TLorentzVector &electron) {
553 // the 'electron' variable will be changed by the function
554 float energy = electron.E(); // before smearing
555 float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
556
557 if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
558 energyS = gRandom->Gaus(energy, sqrt(
559 pow(ELG_Ncen,2) +
560 pow(ELG_Ccen*energy,2) +
561 pow(ELG_Scen*sqrt(energy),2) ));
562 }
563 if(fabs(electron.Eta()) > CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_fwd){
564 energyS = gRandom->Gaus(energy, sqrt(
565 pow(ELG_Nfwd,2) +
566 pow(ELG_Cfwd*energy,2) +
567 pow(ELG_Sfwd*sqrt(energy),2) ) );
568 }
569 electron.SetPtEtaPhiE(energyS/cosh(electron.Eta()), electron.Eta(), electron.Phi(), energyS);
570 if(electron.E() < 0)electron.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
571}
572
573
574// **********Provides the smeared TLorentzVector for the muons********
575// Smears the muon pT and changes the 4-momentum accordingly
576void RESOLution::SmearMu(TLorentzVector &muon) {
577 // the 'muon' variable will be changed by the function
578 float pt = muon.Pt(); // before smearing
579 float ptS=pt;
580
581 if(fabs(muon.Eta()) < CEN_max_mu )
582 {
583 ptS = gRandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
584 }
585 muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
586
587 if(muon.E() < 0)muon.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
588}
589
590
591// **********Provides the smeared TLorentzVector for the hadrons********
592// Smears the hadron 4-momentum
593void RESOLution::SmearHadron(TLorentzVector &hadron, const float frac)
594 // the 'hadron' variable will be changed by the function
595 // the 'frac' variable describes the long-living particles. Should be 0.7 for K0S and Lambda, 1. otherwise
596{
597 float energy = hadron.E(); // before smearing
598 float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
599 float energy_ecal = (1.0 - frac)*energy; // electromagnetic calorimeter
600 float energy_hcal = frac*energy; // hadronic calorimeter
601 // frac takes into account the decay of long-living particles, that decay in the calorimeters
602 // some of the particles decay mostly in the ecal, some mostly in the hcal
603
604 float energyS1,energyS2;
605 if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
606 energyS1 = gRandom->Gaus(energy_hcal, sqrt(
607 pow(HAD_Nhcal,2) +
608 pow(HAD_Chcal*energy_hcal,2) +
609 pow(HAD_Shcal*sqrt(energy_hcal),2) )) ;
610
611
612 energyS2 = gRandom->Gaus(energy_ecal, sqrt(
613 pow(ELG_Ncen,2) +
614 pow(ELG_Ccen*energy_ecal,2) +
615 pow(ELG_Scen*sqrt(energy_ecal),2) ) );
616
617 energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
618 }
619 if(abs(hadron.Eta()) > CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_fwd){
620 energyS = gRandom->Gaus(energy, sqrt(
621 pow(HAD_Nhf,2) +
622 pow(HAD_Chf*energy,2) +
623 pow(HAD_Shf*sqrt(energy),2) ));
624}
625
626
627
628 hadron.SetPtEtaPhiE(energyS/cosh(hadron.Eta()),hadron.Eta(), hadron.Phi(), energyS);
629
630 if(hadron.E() < 0)hadron.SetPxPyPzE(0,0,0,0);
631}
632
633//******************************************************************************************
634
635void RESOLution::SortedVector(vector<ParticleUtil> &vect)
636{
637 int i,j = 0;
638 TLorentzVector tmp;
639 bool en_desordre = true;
640 int entries=vect.size();
641 for(i = 0 ; (i < entries) && en_desordre; i++)
642 {
643 en_desordre = false;
644 for(j = 1 ; j < entries - i ; j++)
645 {
646 if ( vect[j].Pt() > vect[j-1].Pt() )
647 {
648 ParticleUtil tmp = vect[j-1];
649 vect[j-1] = vect[j];
650 vect[j] = tmp;
651 en_desordre = true;
652 }
653 }
654 }
655}
656
657// **********Provides the energy in the cone of radius TAU_CONE_ENERGY for the tau identification********
658// to be taken into account, a calo tower should
659// 1) have a transverse energy \f$ E_T = \sqrt{E_X^2 + E_Y^2} \f$ above a given threshold
660// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
661double RESOLution::EnergySmallCone(const vector<PhysicsTower> &towers, const float eta, const float phi) {
662 double Energie=0;
663 for(unsigned int i=0; i < towers.size(); i++) {
664 if(towers[i].fourVector.pt() < JET_seed) continue;
665 if((DeltaR(phi,eta,towers[i].fourVector.phi(),towers[i].fourVector.eta()) < TAU_energy_scone)) {
666 Energie += towers[i].fourVector.E;
667 }
668 }
669 return Energie;
670}
671
672
673// **********Provides the number of tracks in the cone of radius TAU_CONE_TRACKS for the tau identification********
674// to be taken into account, a track should
675// 1) avec a transverse momentum \$f p_T \$ above a given threshold
676// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
677// IMPORTANT REMARK !!!!!
678// previously, the argument 'phi' was before the argument 'eta'
679// this has been changed for consistency with the other functions
680// double check your running code that uses NumTracks !
681unsigned int RESOLution::NumTracks(const vector<TLorentzVector> &tracks, const float pt_track, const float eta, const float phi) {
682 unsigned int numtrack=0;
683 for(unsigned int i=0; i < tracks.size(); i++) {
684 if((tracks[i].Pt() < pt_track )||
685 (DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()) > TAU_track_scone)
686 )continue;
687 numtrack++;
688 }
689 return numtrack;
690}
691
692
693//*** Returns the PID of the particle with the highest energy, in a cone with a radius CONERADIUS and an axis (eta,phi) *********
694//used by Btaggedjet
695///// Attention : bug removed => CONERADIUS/2 -> CONERADIUS !!
696int RESOLution::Bjets(const TSimpleArray<TRootGenParticle> &subarray, const float eta, const float phi) {
697 float emax=0;
698 int Ppid=0;
699 if(subarray.GetEntries()>0) {
700 for(int i=0; i < subarray.GetEntries();i++) { // should have pt>PT_JETMIN and a small cone radius (r<CONE_JET)
701 float genDeltaR = DeltaR(subarray[i]->Phi,subarray[i]->Eta,phi,eta);
702 if(genDeltaR < JET_coneradius && subarray[i]->E > emax) {
703 emax=subarray[i]->E;
704 Ppid=abs(subarray[i]->PID);
705 }
706 }
707 }
708 return Ppid;
709}
710
711
712//******************** Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************
713bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootGenParticle> &subarray) {
714 if( rand()%100 < (BTAG_b+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
715 else if( rand()%100 < (BTAG_mistag_c+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
716 else if( rand()%100 < (BTAG_mistag_l+1) && Bjets(subarray,JET.Eta(),JET.Phi())!=0) return true; // b-tag of light jets is 1%
717 return false;
718}
719
720//***********************Isolation criteria***********************
721//****************************************************************
722bool RESOLution::Isolation(Float_t phi,Float_t eta,const vector<TLorentzVector> &tracks,float PT_TRACK2)
723{
724 bool isolated = false;
725 Float_t deltar=5000.; // Initial value; should be high; no further repercussion
726 // loop on all final charged particles, with p_t >2, close enough from the electron
727 for(unsigned int i=0; i < tracks.size(); i++)
728 {
729 if(tracks[i].Pt() < PT_TRACK2)continue;
730 Float_t genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()); // slower to evaluate
731 if(
732 (genDeltaR > deltar) ||
733 (genDeltaR==0)
734 ) continue ;
735 deltar=genDeltaR;
736 }
737 if(deltar > 0.5)isolated = true; // returns the closest distance
738 return isolated;
739}
740
741
742 //********** returns a segmented value for eta and phi, for calo towers *****
743void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
744 iEta = -100;
745 int index=-100;
746 for (unsigned int i=1; i< TOWER_number+1; i++) {
747 if(fabs(eta)>TOWER_eta_edges[i-1] && fabs(eta)<TOWER_eta_edges[i]) {
748 iEta = (eta>0) ? TOWER_eta_edges[i-1] : -TOWER_eta_edges[i];
749 index = i-1;
750 //cout << setw(15) << left << eta << "\t" << iEta << endl;
751 break;
752 }
753 }
754 if(index==-100) return;
755 iPhi = -100;
756 float dphi = TOWER_dphi[index]*PI/180.;
757 for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
758 float low = -PI+(i-1)*dphi;
759 float high= low+dphi;
760 if(phi > low && phi < high ){
761 iPhi = low;
762 break;
763 }
764 }
765 if (phi > PI-dphi) iPhi = PI-dphi;
766}
767
768//**************************** Returns the delta Phi ****************************
769float DeltaPhi(const float phi1, const float phi2) {
770 float deltaphi=phi1-phi2; // in here, -PI < phi < PI
771 if(fabs(deltaphi) > PI) deltaphi=2.*PI-fabs(deltaphi);// put deltaphi between 0 and PI
772 else deltaphi=fabs(deltaphi);
773
774 return deltaphi;
775}
776
777//**************************** Returns the delta R****************************
778float DeltaR(const float phi1, const float eta1, const float phi2, const float eta2) {
779 return sqrt(pow(DeltaPhi(phi1,phi2),2) + pow(eta1-eta2,2));
780}
781
782int sign(const int myint) {
783 if (myint >0) return 1;
784 else if (myint <0) return -1;
785 else return 0;
786}
787
788int sign(const float myfloat) {
789 if (myfloat >0) return 1;
790 else if (myfloat <0) return -1;
791 else return 0;
792}
793
794int Charge(int pid)
795{
796 int charge;
797 if(
798 (pid == pGAMMA) ||
799 (pid == pPI0) ||
800 (pid == pK0L) ||
801 (pid == pN) ||
802 (pid == pSIGMA0) ||
803 (pid == pDELTA0) ||
804 (pid == pK0S) // not charged particles : invisible by tracker
805 )
806 charge = 0;
807 else charge = (sign(pid));
808 return charge;
809
810}
Note: See TracBrowser for help on using the repository browser.