Fork me on GitHub

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

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

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