---

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"


namespace Atlfast {

  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.;
    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,pt;
  float sqrtene = sqrt(vec.e());
  float abseta  = fabs(vec.pseudoRapidity());
  float rcone;
  pt = vec.perp();
  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;

  }

}//end of namespace bracket

---