ElectronSmearer.cxx

Go to the documentation of this file.

// ElectronSmearer.cxx
//
// Implementation of the Electron Smearer class
//
// Only the smear() method is overridden
//
//
// Authors: H.T. Phillips, P.Clarke, E.Richter-Was, P.Sherwood, R.Steward
// Update:  Sabine Elles, Simon Dean
//
#include "AtlfastAlgs/ElectronSmearer.h"
#include "AtlfastAlgs/CorrectionFactorElectron.h"

#include <cmath>

#include "CLHEP/Vector/LorentzVector.h"


#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RandGauss.h"
#include "CLHEP/Units/SystemOfUnits.h"

namespace Atlfast {

  HepLorentzVector ElectronSmearer::smear(const HepMC::GenParticle& particle){
    return smear(particle.momentum());
  }
  
  HepLorentzVector ElectronSmearer::smear (const HepLorentzVector& avec) {
    
    double rpilup = 0.0;
    double aa, bb, sigph1, sigph, sigpu;
    //Make sure momentum, aenergy units are Gev
    double ene     = avec.e()/GeV ;
    double sqrtene = sqrt(ene);
    double pt      = avec.perp()/GeV ;
    double vEta    = fabs(avec.pseudoRapidity());
    HepLorentzVector bvec(0,0,0,0);
    
    if(m_smearParamSchema==1){ // TDR-style smearing
      
      // do the smearing, copied verbatim (except for ROOT dependencies)
      // from electronmaker codein Atlfast++
      //
      // This code has not otherwise been altered which is why it is all very
      // procedural
      //
      // Parametrisation was by L.Poggioli
      
      while (1) {
        aa=randGauss()->fire();
        //bb=randGauss()->fire();
        sigph1 = aa*m_smearParams[0]/sqrtene;
        if (1.0 + sigph1 <= 0.0) continue;
        sigph = sigph1;
        if (vEta < 1.4) {
          while (1) {
            aa=randGauss()->fire();
            bb=randGauss()->fire();
            sigph1 = aa*m_smearParams[1]/pt + bb*m_smearParams[2];
            if (1.0+sigph1 > 0) break;
          }
        } else {
          while (1) {
            aa=randGauss()->fire();
            bb=randGauss()->fire();
            sigph1 = aa*m_smearParams[3]*((m_smearParams[4]-vEta)+m_smearParams[5])/ene + bb*m_smearParams[6];
            //    sigph1 = aa*0.306*((2.4-vEta)+0.228)/ene + bb*0.007;
            if (1.0+sigph1 > 0) break;
          }
        }
        sigph += sigph1;
        if (m_lumi == 2) {
          if (vEta < 0.6) rpilup = 0.32;
          if (vEta > 0.6 && vEta < 1.4) rpilup = 0.295;
          if (vEta > 1.4) rpilup = 0.27;
          while (1) {
            aa=randGauss()->fire();
            //bb=randGauss()->fire();
            sigpu = aa*rpilup/pt;
            if (1.0+sigpu > 0) break;
          }
          sigph += sigpu;
        }
        if (1.0 + sigph > 0) break;
      }
      
      //now sigph is the electron "sigma" and we do the following a la Atlfast++ electron maker:
      bvec.set(avec.px()*(1.0+sigph), avec.py()*(1.0+sigph), avec.pz()*(1.0+sigph), avec.e() *(1.0+sigph));
      
    } else if (m_smearParamSchema==2){ // Smearing based on CSC sample tuning
    
      double vE = avec.e();

      double vLimEta=EtaElectron[NbEtaElectron-1]+(EtaElectron[NbEtaElectron-1]-EtaElectron[NbEtaElectron-2])*0.5;
    
      if( vEta>vLimEta )
        {
          m_log<<MSG::DEBUG<<"No calibration : "<<vEta<<endreq;
          return avec;
        }
    
      int iEnergy=0;
      m_log<<MSG::DEBUG<<vE<<" > "<<vEnergiesElectron[iEnergy]<<" ";
      while(iEnergy<NbEnergiesElectron&&vE>vEnergiesElectron[iEnergy]*1000)
        {
          m_log<<MSG::DEBUG<<vEnergiesElectron[iEnergy]<<" ";
          iEnergy++;
        }
      m_log<<MSG::DEBUG<<endreq;
      if(iEnergy==0)iEnergy++;
      double ratioEnergy=(vE*0.001-vEnergiesElectron[iEnergy-1])/(vEnergiesElectron[iEnergy]-vEnergiesElectron[iEnergy-1]);

      int iEta=0;
      m_log<<MSG::DEBUG<<vEta<<" > "<<EtaElectron[iEta]<<" ";
      while(iEta<NbEtaElectron&&vEta>EtaElectron[iEta])
        {
          m_log<<MSG::DEBUG<<EtaElectron[iEta]<<" ";
          iEta++;
        }
      m_log<<MSG::DEBUG<<endreq;

      if(iEta==0)iEta++;
      double ratioEta=(vEta-EtaElectron[iEta-1])/(EtaElectron[iEta]-EtaElectron[iEta-1]);

      int iPos;
      double sigma,sigma1,sigma2;

      iPos=iEnergy-1;
      sigma1=CorrFactorElectronSigma[iPos][iEta-1]+ratioEta*(CorrFactorElectronSigma[iPos][iEta]-CorrFactorElectronSigma[iPos][iEta-1]);
      iPos=iEnergy;
      sigma2=CorrFactorElectronSigma[iPos][iEta-1]+ratioEta*(CorrFactorElectronSigma[iPos][iEta]-CorrFactorElectronSigma[iPos][iEta-1]);

      sigma=sigma1+ratioEnergy*(sigma2-sigma1);

      m_log<<MSG::DEBUG<<" -> interpolation sigma "<<sigma1<<" "<<sigma2<<" "<<sigma<<endreq;

      double alpha=1.0;
      double vrandom=alpha;

      int iCmpt=int(randFlat()->fire()*20);
      for(int i=0; i<iCmpt; i++)
        vrandom=randGauss()->fire(1.0,sigma);
    
      m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesElectron[iEnergy]<<"   Eta : "<<vEnergiesElectron[iEnergy-1]<<" < "<<vE<<" < "<<vEnergiesElectron[iEnergy]<<endreq;
      m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesElectron[iEnergy]<<"   Eta : "<<EtaElectron[iEta-1]<<" < "<<vEta<<" < "<<EtaElectron[iEta]<<endreq;
      m_log<<MSG::DEBUG<<"       calib mean-sigma  "<<alpha<<" "<<sigma<<" -> "<<vrandom<<endreq;

      bvec.set(avec.px()*vrandom, avec.py()*vrandom, avec.pz()*vrandom, vE*vrandom);
    
      m_log<<MSG::DEBUG<<"-> HepLorentzVector calibration "<<vEta<<" "<<avec<<"  ->  "<<bvec<<endreq;

    }

    return bvec;
  } //end of smear method
  
  //cct: implement the setSmearParameters method
  int ElectronSmearer::setSmearParameters (const std::vector<double>& smearValues){
    m_smearParams=smearValues;
    return 0;
  }
  
  //cct: implement the setSmearParamSchema method
  int ElectronSmearer::setSmearParamSchema ( const int smearSchema){
    m_smearParamSchema=smearSchema;
    return 0;
  }
  
} // end of namespace bracket

---