Fork me on GitHub

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

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

sorted vector
File size: 40.0 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/*
580void RESOLution::SortedVector()
581{
582 int i,j = 0;
583 TLorentzVector tmp;
584 bool en_desordre = true;
585 for(i = 0 ; (i < numjetT) && en_desordre; i++)
586 {
587 en_desordre = false;
588 for(j = 1 ; j < numjetT - i ; j++)
589 {
590 if ( Wjets[j].Eta() > Wjets[j-1].Eta() )
591 {
592 TLorentzVector tmp = Wjets[j-1];
593 Wjets[j-1] = Wjets[j];
594 Wjets[j] = tmp;
595 en_desordre = true;
596 }
597 }
598 }
599
600}*/
601
602//******************************************************************************************
603
604void RESOLution::SortedVector(vector<ParticleUtil> &vect)
605{
606 int i,j = 0;
607 TLorentzVector tmp;
608 bool en_desordre = true;
609 int entries=vect.size();
610 for(i = 0 ; (i < entries) && en_desordre; i++)
611 {
612 en_desordre = false;
613 for(j = 1 ; j < entries - i ; j++)
614 {
615 if ( vect[j].Pt() > vect[j-1].Pt() )
616 {
617 ParticleUtil tmp = vect[j-1];
618 vect[j-1] = vect[j];
619 vect[j] = tmp;
620 en_desordre = true;
621 }
622 }
623 }
624}
625
626// **********Provides the energy in the cone of radius TAU_CONE_ENERGY for the tau identification********
627// to be taken into account, a calo tower should
628// 1) have a transverse energy \f$ E_T = \sqrt{E_X^2 + E_Y^2} \f$ above a given threshold
629// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
630double RESOLution::EnergySmallCone(const vector<PhysicsTower> &towers, const float eta, const float phi) {
631 double Energie=0;
632 for(unsigned int i=0; i < towers.size(); i++) {
633 if(towers[i].fourVector.pt() < SEEDTHRESHOLD) continue;
634 if((DeltaR(phi,eta,towers[i].fourVector.phi(),towers[i].fourVector.eta()) < TAU_CONE_ENERGY)) {
635 Energie += towers[i].fourVector.E;
636 }
637 }
638 return Energie;
639}
640
641
642// **********Provides the number of tracks in the cone of radius TAU_CONE_TRACKS for the tau identification********
643// to be taken into account, a track should
644// 1) avec a transverse momentum \$f p_T \$ above a given threshold
645// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
646// IMPORTANT REMARK !!!!!
647// previously, the argument 'phi' was before the argument 'eta'
648// this has been changed for consistency with the other functions
649// double check your running code that uses NumTracks !
650unsigned int RESOLution::NumTracks(const vector<TLorentzVector> &tracks, const float pt_track, const float eta, const float phi) {
651 unsigned int numtrack=0;
652 for(unsigned int i=0; i < tracks.size(); i++) {
653 if((tracks[i].Pt() < pt_track )||
654 (DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()) > TAU_CONE_TRACKS)
655 )continue;
656 numtrack++;
657 }
658 return numtrack;
659}
660
661
662//*** Returns the PID of the particle with the highest energy, in a cone with a radius CONERADIUS and an axis (eta,phi) *********
663//used by Btaggedjet
664///// Attention : bug removed => CONERADIUS/2 -> CONERADIUS !!
665int RESOLution::Bjets(const TSimpleArray<TRootGenParticle> &subarray, const float eta, const float phi) {
666 float emax=0;
667 int Ppid=0;
668 if(subarray.GetEntries()>0) {
669 for(int i=0; i < subarray.GetEntries();i++) { // should have pt>PT_JETMIN and a small cone radius (r<CONE_JET)
670 float genDeltaR = DeltaR(subarray[i]->Phi,subarray[i]->Eta,phi,eta);
671 if(genDeltaR < CONERADIUS && subarray[i]->E > emax) {
672 emax=subarray[i]->E;
673 Ppid=abs(subarray[i]->PID);
674 }
675 }
676 }
677 return Ppid;
678}
679
680
681//******************** Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************
682bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootGenParticle> &subarray) {
683 if( rand()%100 < (TAGGING_B+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
684 else if( rand()%100 < (MISTAGGING_C+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
685 else if( rand()%100 < (MISTAGGING_L+1) && Bjets(subarray,JET.Eta(),JET.Phi())!=0) return true; // b-tag of light jets is 1%
686 return false;
687}
688
689//***********************Isolation criteria***********************
690//****************************************************************
691bool RESOLution::Isolation(Float_t phi,Float_t eta,const vector<TLorentzVector> &tracks,float PT_TRACK2)
692{
693 bool isolated = false;
694 Float_t deltar=5000.; // Initial value; should be high; no further repercussion
695 // loop on all final charged particles, with p_t >2, close enough from the electron
696 for(unsigned int i=0; i < tracks.size(); i++)
697 {
698 if(tracks[i].Pt() < PT_TRACK2)continue;
699 Float_t genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()); // slower to evaluate
700 if(
701 (genDeltaR > deltar) ||
702 (genDeltaR==0)
703 ) continue ;
704 deltar=genDeltaR;
705 }
706 if(deltar > 0.5)isolated = true; // returns the closest distance
707 return isolated;
708}
709
710
711 //********** returns a segmented value for eta and phi, for calo towers *****
712void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
713 iEta = -100;
714 int index=-100;
715 for (unsigned int i=1; i< NTOWERS+1; i++) {
716 if(fabs(eta)>TOWER_ETA_EDGES[i-1] && fabs(eta)<TOWER_ETA_EDGES[i]) {
717 iEta = (eta>0) ? TOWER_ETA_EDGES[i-1] : -TOWER_ETA_EDGES[i];
718 index = i-1;
719 //cout << setw(15) << left << eta << "\t" << iEta << endl;
720 break;
721 }
722 }
723 if(index==-100) return;
724 iPhi = -100;
725 float dphi = TOWER_DPHI[index]*PI/180.;
726 for (unsigned int i=1; i < 360/TOWER_DPHI[index]; i++ ) {
727 float low = -PI+(i-1)*dphi;
728 float high= low+dphi;
729 if(phi > low && phi < high ){
730 iPhi = low;
731 break;
732 }
733 }
734 if (phi > PI-dphi) iPhi = PI-dphi;
735}
736
737//**************************** Returns the delta Phi ****************************
738float DeltaPhi(const float phi1, const float phi2) {
739 float deltaphi=phi1-phi2; // in here, -PI < phi < PI
740 if(fabs(deltaphi) > PI) deltaphi=2.*PI-fabs(deltaphi);// put deltaphi between 0 and PI
741 else deltaphi=fabs(deltaphi);
742
743 return deltaphi;
744}
745
746//**************************** Returns the delta R****************************
747float DeltaR(const float phi1, const float eta1, const float phi2, const float eta2) {
748 return sqrt(pow(DeltaPhi(phi1,phi2),2) + pow(eta1-eta2,2));
749}
750
751int sign(const int myint) {
752 if (myint >0) return 1;
753 else if (myint <0) return -1;
754 else return 0;
755}
756
757int sign(const float myfloat) {
758 if (myfloat >0) return 1;
759 else if (myfloat <0) return -1;
760 else return 0;
761}
762
763int Charge(int pid)
764{
765 int charge;
766 if(
767 (pid == pGAMMA) ||
768 (pid == pPI0) ||
769 (pid == pK0L) ||
770 (pid == pN) ||
771 (pid == pSIGMA0) ||
772 (pid == pDELTA0) ||
773 (pid == pK0S) // not charged particles : invisible by tracker
774 )
775 charge = 0;
776 else charge = (sign(pid));
777 return charge;
778
779}
Note: See TracBrowser for help on using the repository browser.