JetSmearer.cxx

Go to the documentation of this file.

// JetSmearer.cxx
//
// Implementation of the Jet Smearer class
//
// Only the smear() method is overridden
//
//
// Authors: H.T. Phillips, P.Clarke, E.Richter-Was, P.Sherwood, R.Steward
//
//

#include "AtlfastAlgs/JetSmearer.h"
#include "AtlfastAlgs/PhotonSmearer.h"
#include <cmath>
#include <iostream>

#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RandGauss.h"

#include "CLHEP/Units/SystemOfUnits.h"

namespace Atlfast {

  HepLorentzVector JetSmearer::smear(const HepMC::GenParticle& particle){
    return smear(particle.momentum());
  }

  HepLorentzVector JetSmearer::smear (const HepLorentzVector& vec) {
    // do the smearing, copied verbatim (except for ROOT dependencies)
    // from Jetmaker codein Atlfast++
    //
    // This code has not otherwise been altered which is why it is all very
    // procedural
    //
    
    //........................................................
    //.....smear clusters energy
    //........................................................
    float sigma=0.;
    if ( vec.e() == 0 ) return vec;
    HepLorentzVector smearedVec(vec);
    
    
    //     parametrizes smearing for hadronic energy deposition
    //     parametrization from L. Poggioli
    //     pile-up added as in HADCALO
    
    float epileup = 0;
    float aa, bb, cc, dd;
    float abseta  = fabs(vec.pseudoRapidity());
    float rcone;
    //make sure momentum, energy is in GeV
    float sqrtene = sqrt(vec.e()/GeV);
    float pt      = vec.perp()/GeV;
    if (fabs(vec.pseudoRapidity()) <= m_BarrelForwardEta) {rcone=m_rconeb;}else{rcone=m_rconef;}
    if(rcone <= 0.4)                epileup = 0.4;
    if(rcone == 0.4)                epileup = 7.5;
    if(rcone > 0.4 && rcone <= 0.5) epileup = 12.0;
    if(rcone > 0.5 && rcone <= 0.7) epileup = 18.0;
    if(rcone > 0.7)                 epileup = 20.0;
    
    if(m_lumi <= 1) {
      while(1) {
        aa=randGauss()->fire();
        bb=randGauss()->fire();
        if(abseta < m_BarrelForwardEta) sigma = aa*0.5/sqrtene + bb*0.03;
        else            sigma = aa*1.0/sqrtene + bb*0.07;
        if(1.+sigma > .0) break;
      }
    } else if(m_lumi == 2) {
      while(1) {
        aa=randGauss()->fire();
        bb=randGauss()->fire();
        cc=randGauss()->fire();
        dd=randGauss()->fire();
        if(abseta < m_BarrelForwardEta) sigma = aa*0.5/sqrtene + bb*0.03  + cc*epileup/pt;
        else            sigma = aa*1.0/sqrtene + bb*0.07  + cc*epileup/pt;
        if(1.+sigma > .0) break;
      }
      
    }                                
    smearedVec.setPx(vec.px()*(1.0+sigma));
    smearedVec.setPy(vec.py()*(1.0+sigma));
    smearedVec.setPz(vec.pz()*(1.0+sigma));
    smearedVec.setE(vec.e()*(1.0+sigma));
    
    
    return smearedVec;
    
  }
  
  //cct: implement the setSmearParameters method
  int JetSmearer::setSmearParameters (const std::vector<double>& /*smearValues*/){
    return 0;
  }
  //cct: implement the setSmearParamSchema method
  int JetSmearer::setSmearParamSchema ( const int /*smearSchema*/){
    return 0;
  }
  
}//end of namespace bracket

---