Fork me on GitHub

Changeset 219 in svn for trunk/src/SmearUtil.cc


Ignore:
Timestamp:
Feb 2, 2009, 12:33:21 PM (16 years ago)
Author:
Xavier Rouby
Message:

JetUtils.cc

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/SmearUtil.cc

    r177 r219  
    1313
    1414
    15 #include "interface/SmearUtil.h"
     15#include "SmearUtil.h"
    1616#include "TRandom.h"
    1717
    1818#include <iostream>
     19#include <fstream>
    1920#include <sstream>
    20 #include <fstream>
    2121#include <iomanip>
    22 
    23 
    24 
    2522using namespace std;
     23
     24
     25ParticleUtil::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}
    2637
    2738//------------------------------------------------------------------------------
     
    147158}
    148159
     160
     161RESOLution::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
     265RESOLution& 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
    149374//------------------------------------------------------------------------------
    150375void RESOLution::ReadDataCard(const string datacard) {
     
    155380  ifstream fichier_a_lire(datacard.c_str());
    156381  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;
    158383    return;
    159384  }
     
    252477}
    253478
    254 void RESOLution::Logfile(string LogName) {
     479void RESOLution::Logfile(const string& LogName) {
    255480  //void RESOLution::Logfile(string outputfilename) {
    256481 
     
    617842  energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
    618843  }
    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){
    620845    energyS = gRandom->Gaus(energy, sqrt(
    621846                            pow(HAD_Nhf,2) +
     
    720945//***********************Isolation criteria***********************
    721946//****************************************************************
    722 bool RESOLution::Isolation(Float_t phi,Float_t eta,const vector<TLorentzVector> &tracks,float PT_TRACK2)
     947bool RESOLution::Isolation(const float phi, const float eta,const vector<TLorentzVector> &tracks, const float pt_second_track)
    723948{
    724949   bool isolated = false;
    725    Float_t deltar=5000.; // Initial value; should be high; no further repercussion
     950   float deltar=5000.; // Initial value; should be high; no further repercussion
    726951   // loop on all final charged particles, with p_t >2, close enough from the electron
    727952   for(unsigned int i=0; i < tracks.size(); i++)
    728953      {
    729        if(tracks[i].Pt() < PT_TRACK2)continue;
    730        Float_t genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta());    // slower to evaluate
     954       if(tracks[i].Pt() < pt_second_track)continue;
     955       float genDeltaR = DeltaR(phi,eta,tracks[i].Phi(),tracks[i].Eta());
    731956         if(
    732957             (genDeltaR > deltar) ||
     
    735960           deltar=genDeltaR;
    736961      }
    737    if(deltar > 0.5)isolated = true; // returns the closest distance
     962   if(deltar > 0.5) isolated = true;
    738963   return isolated;
    739964}
     
    769994float DeltaPhi(const float phi1, const float phi2) {
    770995  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        }
    772999  else deltaphi=fabs(deltaphi);
    7731000 
Note: See TracChangeset for help on using the changeset viewer.