Changeset 219 in svn for trunk/src/SmearUtil.cc
- Timestamp:
- Feb 2, 2009, 12:33:21 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/SmearUtil.cc
r177 r219 13 13 14 14 15 #include " interface/SmearUtil.h"15 #include "SmearUtil.h" 16 16 #include "TRandom.h" 17 17 18 18 #include <iostream> 19 #include <fstream> 19 20 #include <sstream> 20 #include <fstream>21 21 #include <iomanip> 22 23 24 25 22 using namespace std; 23 24 25 ParticleUtil::ParticleUtil(const TLorentzVector &genMomentum, int pid) { 26 _pid=pid; 27 _e = genMomentum.E(); 28 _px = genMomentum.Px(); 29 _py = genMomentum.Py(); 30 _pz = genMomentum.Pz(); 31 _pt = genMomentum.Pt(); 32 33 //_e, _px, _py, _pz, _pt; 34 //float _eta, _etaCalo, _phi, _phiCalo; 35 //int _pid; 36 } 26 37 27 38 //------------------------------------------------------------------------------ … … 147 158 } 148 159 160 161 RESOLution::RESOLution(const RESOLution & DET) { 162 // Detector characteristics 163 CEN_max_tracker = DET.CEN_max_tracker; 164 CEN_max_calo_cen = DET.CEN_max_calo_cen; 165 CEN_max_calo_fwd = DET.CEN_max_calo_fwd; 166 CEN_max_mu = DET.CEN_max_mu; 167 168 // Energy resolution for electron/photon 169 ELG_Scen = DET.ELG_Scen; 170 ELG_Ncen = DET.ELG_Ncen; 171 ELG_Ccen = DET.ELG_Ccen; 172 ELG_Cfwd = DET.ELG_Cfwd; 173 ELG_Sfwd = DET.ELG_Sfwd; 174 ELG_Nfwd = DET.ELG_Nfwd; 175 176 // Energy resolution for hadrons in ecal/hcal/hf 177 HAD_Shcal = DET.HAD_Shcal; 178 HAD_Nhcal = DET.HAD_Nhcal; 179 HAD_Chcal = DET.HAD_Chcal; 180 HAD_Shf = DET.HAD_Shf; 181 HAD_Nhf = DET.HAD_Nhf; 182 HAD_Chf = DET.HAD_Chf; 183 184 // Muon smearing 185 MU_SmearPt = DET.MU_SmearPt; 186 187 // Tracking efficiencies 188 TRACK_ptmin = DET.TRACK_ptmin; 189 TRACK_eff = DET.TRACK_eff; 190 191 // Calorimetric towers 192 TOWER_number = DET.TOWER_number; 193 TOWER_eta_edges = new float[TOWER_number+1]; 194 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i]; 195 196 TOWER_dphi = new float[TOWER_number]; 197 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i]; 198 199 // Thresholds for reconstructed objetcs 200 PTCUT_elec = DET.PTCUT_elec; 201 PTCUT_muon = DET.PTCUT_muon; 202 PTCUT_jet = DET.PTCUT_jet; 203 PTCUT_gamma = DET.PTCUT_gamma; 204 PTCUT_taujet = DET.PTCUT_taujet; 205 206 // General jet variable 207 JET_coneradius = DET.JET_coneradius; 208 JET_jetalgo = DET.JET_jetalgo; 209 JET_seed = DET.JET_seed; 210 211 // Tagging definition 212 BTAG_b = DET.BTAG_b; 213 BTAG_mistag_c = DET.BTAG_mistag_c; 214 BTAG_mistag_l = DET.BTAG_mistag_l; 215 216 // FLAGS 217 FLAG_bfield = DET.FLAG_bfield; 218 FLAG_vfd = DET.FLAG_vfd; 219 FLAG_trigger = DET.FLAG_trigger; 220 FLAG_frog = DET.FLAG_frog; 221 222 // In case BField propagation allowed 223 TRACK_radius = DET.TRACK_radius; 224 TRACK_length = DET.TRACK_length; 225 TRACK_bfield_x = DET.TRACK_bfield_x; 226 TRACK_bfield_y = DET.TRACK_bfield_y; 227 TRACK_bfield_z = DET.TRACK_bfield_z; 228 229 // In case Very forward detectors allowed 230 VFD_min_calo_vfd = DET.VFD_min_calo_vfd; 231 VFD_max_calo_vfd = DET.VFD_max_calo_vfd; 232 VFD_min_zdc = DET.VFD_min_zdc; 233 VFD_s_zdc = DET.VFD_s_zdc; 234 235 RP_220_s = DET.RP_220_s; 236 RP_220_x = DET.RP_220_x; 237 RP_420_s = DET.RP_420_s; 238 RP_420_x = DET.RP_420_x; 239 240 // In case FROG event display allowed 241 NEvents_Frog = DET.NEvents_Frog; 242 243 JET_overlap = DET.JET_overlap; 244 // MidPoint algorithm definition 245 JET_M_coneareafraction = DET.JET_M_coneareafraction; 246 JET_M_maxpairsize = DET.JET_M_maxpairsize; 247 JET_M_maxiterations = DET.JET_M_maxiterations; 248 // Define Cone algorithm. 249 JET_C_adjacencycut = DET.JET_C_adjacencycut; 250 JET_C_maxiterations = DET.JET_C_maxiterations; 251 JET_C_iratch = DET.JET_C_iratch; 252 //Define SISCone algorithm. 253 JET_S_npass = DET.JET_S_npass; 254 JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin; 255 256 //For Tau-jet definition 257 TAU_energy_scone = DET.TAU_energy_scone; 258 TAU_track_scone = DET.TAU_track_scone; 259 TAU_track_pt = DET.TAU_track_pt; 260 TAU_energy_frac = DET.TAU_energy_frac; 261 262 PT_QUARKS_MIN = DET.PT_QUARKS_MIN; 263 } 264 265 RESOLution& RESOLution::operator=(const RESOLution& DET) { 266 if(this==&DET) return *this; 267 // Detector characteristics 268 CEN_max_tracker = DET.CEN_max_tracker; 269 CEN_max_calo_cen = DET.CEN_max_calo_cen; 270 CEN_max_calo_fwd = DET.CEN_max_calo_fwd; 271 CEN_max_mu = DET.CEN_max_mu; 272 273 // Energy resolution for electron/photon 274 ELG_Scen = DET.ELG_Scen; 275 ELG_Ncen = DET.ELG_Ncen; 276 ELG_Ccen = DET.ELG_Ccen; 277 ELG_Cfwd = DET.ELG_Cfwd; 278 ELG_Sfwd = DET.ELG_Sfwd; 279 ELG_Nfwd = DET.ELG_Nfwd; 280 281 // Energy resolution for hadrons in ecal/hcal/hf 282 HAD_Shcal = DET.HAD_Shcal; 283 HAD_Nhcal = DET.HAD_Nhcal; 284 HAD_Chcal = DET.HAD_Chcal; 285 HAD_Shf = DET.HAD_Shf; 286 HAD_Nhf = DET.HAD_Nhf; 287 HAD_Chf = DET.HAD_Chf; 288 289 // Muon smearing 290 MU_SmearPt = DET.MU_SmearPt; 291 292 // Tracking efficiencies 293 TRACK_ptmin = DET.TRACK_ptmin; 294 TRACK_eff = DET.TRACK_eff; 295 296 // Calorimetric towers 297 TOWER_number = DET.TOWER_number; 298 TOWER_eta_edges = new float[TOWER_number+1]; 299 for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i]; 300 301 TOWER_dphi = new float[TOWER_number]; 302 for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i]; 303 304 // Thresholds for reconstructed objetcs 305 PTCUT_elec = DET.PTCUT_elec; 306 PTCUT_muon = DET.PTCUT_muon; 307 PTCUT_jet = DET.PTCUT_jet; 308 PTCUT_gamma = DET.PTCUT_gamma; 309 PTCUT_taujet = DET.PTCUT_taujet; 310 311 // General jet variable 312 JET_coneradius = DET.JET_coneradius; 313 JET_jetalgo = DET.JET_jetalgo; 314 JET_seed = DET.JET_seed; 315 316 // Tagging definition 317 BTAG_b = DET.BTAG_b; 318 BTAG_mistag_c = DET.BTAG_mistag_c; 319 BTAG_mistag_l = DET.BTAG_mistag_l; 320 321 // FLAGS 322 FLAG_bfield = DET.FLAG_bfield; 323 FLAG_vfd = DET.FLAG_vfd; 324 FLAG_trigger = DET.FLAG_trigger; 325 FLAG_frog = DET.FLAG_frog; 326 327 // In case BField propagation allowed 328 TRACK_radius = DET.TRACK_radius; 329 TRACK_length = DET.TRACK_length; 330 TRACK_bfield_x = DET.TRACK_bfield_x; 331 TRACK_bfield_y = DET.TRACK_bfield_y; 332 TRACK_bfield_z = DET.TRACK_bfield_z; 333 334 // In case Very forward detectors allowed 335 VFD_min_calo_vfd = DET.VFD_min_calo_vfd; 336 VFD_max_calo_vfd = DET.VFD_max_calo_vfd; 337 VFD_min_zdc = DET.VFD_min_zdc; 338 VFD_s_zdc = DET.VFD_s_zdc; 339 340 RP_220_s = DET.RP_220_s; 341 RP_220_x = DET.RP_220_x; 342 RP_420_s = DET.RP_420_s; 343 RP_420_x = DET.RP_420_x; 344 345 // In case FROG event display allowed 346 NEvents_Frog = DET.NEvents_Frog; 347 348 JET_overlap = DET.JET_overlap; 349 // MidPoint algorithm definition 350 JET_M_coneareafraction = DET.JET_M_coneareafraction; 351 JET_M_maxpairsize = DET.JET_M_maxpairsize; 352 JET_M_maxiterations = DET.JET_M_maxiterations; 353 // Define Cone algorithm. 354 JET_C_adjacencycut = DET.JET_C_adjacencycut; 355 JET_C_maxiterations = DET.JET_C_maxiterations; 356 JET_C_iratch = DET.JET_C_iratch; 357 //Define SISCone algorithm. 358 JET_S_npass = DET.JET_S_npass; 359 JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin; 360 361 //For Tau-jet definition 362 TAU_energy_scone = DET.TAU_energy_scone; 363 TAU_track_scone = DET.TAU_track_scone; 364 TAU_track_pt = DET.TAU_track_pt; 365 TAU_energy_frac = DET.TAU_energy_frac; 366 367 PT_QUARKS_MIN = DET.PT_QUARKS_MIN; 368 return *this; 369 } 370 371 372 373 149 374 //------------------------------------------------------------------------------ 150 375 void RESOLution::ReadDataCard(const string datacard) { … … 155 380 ifstream fichier_a_lire(datacard.c_str()); 156 381 if(!fichier_a_lire.good()) { 157 cout << datacard << "Datadard " << datacard << " not found, use default values" << endl;382 cout <<"** Datadard not found, use default values **" << endl; 158 383 return; 159 384 } … … 252 477 } 253 478 254 void RESOLution::Logfile( stringLogName) {479 void RESOLution::Logfile(const string& LogName) { 255 480 //void RESOLution::Logfile(string outputfilename) { 256 481 … … 617 842 energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0); 618 843 } 619 if( abs(hadron.Eta()) > CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_fwd){844 if(fabs(hadron.Eta()) > CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_fwd){ 620 845 energyS = gRandom->Gaus(energy, sqrt( 621 846 pow(HAD_Nhf,2) + … … 720 945 //***********************Isolation criteria*********************** 721 946 //**************************************************************** 722 bool RESOLution::Isolation( Float_t phi,Float_t eta,const vector<TLorentzVector> &tracks,float PT_TRACK2)947 bool RESOLution::Isolation(const float phi, const float eta,const vector<TLorentzVector> &tracks, const float pt_second_track) 723 948 { 724 949 bool isolated = false; 725 Float_t deltar=5000.; // Initial value; should be high; no further repercussion950 float deltar=5000.; // Initial value; should be high; no further repercussion 726 951 // loop on all final charged particles, with p_t >2, close enough from the electron 727 952 for(unsigned int i=0; i < tracks.size(); i++) 728 953 { 729 if(tracks[i].Pt() < PT_TRACK2)continue;730 Float_t genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()); // slower to evaluate954 if(tracks[i].Pt() < pt_second_track)continue; 955 float genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta()); 731 956 if( 732 957 (genDeltaR > deltar) || … … 735 960 deltar=genDeltaR; 736 961 } 737 if(deltar > 0.5) isolated = true; // returns the closest distance962 if(deltar > 0.5) isolated = true; 738 963 return isolated; 739 964 } … … 769 994 float DeltaPhi(const float phi1, const float phi2) { 770 995 float deltaphi=phi1-phi2; // in here, -PI < phi < PI 771 if(fabs(deltaphi) > PI) deltaphi=2.*PI-fabs(deltaphi);// put deltaphi between 0 and PI 996 if(fabs(deltaphi) > PI) { 997 deltaphi=2.*PI -fabs(deltaphi);// put deltaphi between 0 and PI 998 } 772 999 else deltaphi=fabs(deltaphi); 773 1000
Note:
See TracChangeset
for help on using the changeset viewer.