Fork me on GitHub

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

Last change on this file since 87 was 82, checked in by uid677, 16 years ago

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