Fork me on GitHub

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

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

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