PhotonSmearer.cxx

Go to the documentation of this file.

// PhotonSmearer.cxx
//
// Implementation of the Photon 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/PhotonSmearer.h"
#include "AtlfastAlgs/CorrectionFactorPhoton.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 {
  using std::abs;
  
  HepLorentzVector PhotonSmearer::smear(const HepMC::GenParticle& particle){
    return smear(particle.momentum());
  }

  HepLorentzVector PhotonSmearer::smear (const HepLorentzVector& avec) {
    
    double rpilup = 0.0;
    double aa, bb, sigph1, sigph, sigpu;
    double vEta  = fabs(avec.pseudoRapidity());
    
    //make sure e, momentum in GeV
    double pt = avec.perp()/GeV;
    double ene = avec.e()/GeV;
    double sqrtene = sqrt(ene);
    
    // The HLV to be returned
    HepLorentzVector return_vec(0.,0.,0.,0.);
    if (!sqrtene || !ene || !pt) return return_vec;

    if (m_smearParamSchema==1){ // TDR-style smearing
      
      // do the smearing, copied verbatim (except for ROOT dependencies)
      // from photonmaker codein Atlfast++
      // The code is very procedural (even more so than Atlfast++!) and should be tidied later!
      //
      // Parametrisation was by L.Poggioli and F. Gianotti
      
      // do smearing of theta first
      HepLorentzVector bvec(avec);
      
      //this is equivalent to RESTHE in FORTRAN 
      double sigph2 = 0.0;
      while (1) {
        aa=randGauss()->fire();
        if (vEta < 0.8) sigph2                  = aa*m_smearParams[0]/sqrtene;
        //if (vEta < 0.8) sigph2                  = aa*0.065/sqrtene;
        if (vEta >= 0.8 && vEta < 1.4) sigph2 = aa*m_smearParams[1]/sqrtene;  
        //if (vEta >= 0.8 && vEta < 1.4) sigph2 = aa*0.050/sqrtene;
        if (vEta >= 1.4 && vEta < 2.5) sigph2 = aa*m_smearParams[2]/sqrtene; 
        //if (vEta >= 1.4 && vEta < 2.5) sigph2 = aa*0.040/sqrtene;
        if (vEta >= 2.5)                 sigph2 = 0;
        if (fabs(bvec.theta()) + sigph2 <= M_PI) break;
      }
      
      // sigph2 is now the offset for theta...
      bvec.setPz(bvec.e()*cos(bvec.theta()+sigph2));
      
      vEta  = fabs(bvec.pseudoRapidity()); //27/02/03 JPC possible bug fix
      
      //this is RESPHO in FORTRAN
      // now do the energy resolution note: use original energy, pt and phi, but new vEta
      //
      while (1) {
        while (1){
          aa=randGauss()->fire();
          sigph1 = aa*m_smearParams[3]/sqrtene;
          //sigph1 = aa*0.10/sqrtene;
          if (1.0 + sigph1 > 0.0) break;
        }
        sigph = sigph1;
        if (vEta < 1.4) {
          while (1) {
            aa=randGauss()->fire();
            bb=randGauss()->fire();
            sigph1 = aa*m_smearParams[4]/pt + bb*m_smearParams[5];
            //sigph1 = aa*0.245/pt + bb*0.007;
            if (1.0+sigph1 > 0) break;
          }
        } else {
          while (1) {
            aa=randGauss()->fire();
            bb=randGauss()->fire();
            sigph1 = aa*m_smearParams[6]*((m_smearParams[7]-vEta)+m_smearParams[8])/ene + bb*m_smearParams[9];
            // 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();
            sigpu = aa*rpilup/pt;
            if (1+sigpu > 0) break;
          }
          sigph += sigpu;
        }
        
        if (1.0 + sigph > 0) break;
      }
      
      
      // Now sigph is the photon "sigma" and we do the following a la Atlfast++ photon maker:
      // this seems to be OK for energy resolution. 
      //
      // It looks to me as though the functional form above for photons is actually the same as
      // electrons, just that one magic number is different. Suggest we revisit this and
      // exploit commonality.
      //
      
      HepLorentzVector cvec;
      cvec.setPx(avec.px()*(1.0+sigph));
      cvec.setPy(avec.py()*(1.0+sigph));
      cvec.setPz(bvec.pz()*(1.0+sigph));//using recalculated z momentum
      
      // nb from FORTRAN code, we go back to the original photon energy here
      cvec.setE( avec.e() *(1.0+sigph));
      //return cvec;
      //now do what the fortran does to get massless particle! YUCK!!!
      
      //if(std::abs (cvec.pz()/cvec.e()) > 1.){
      //        theta = acos(cvec.pz()/sqrt(cvec.px()*cvec.px()+cvec.py()*cvec.py()+cvec.pz()*cvec.pz()));
      //} else if (cvec.pz()>0) theta = 0.;
      //else if (cvec.pz()<0) theta = M_PI;
      
      double ETA = cvec.pseudoRapidity();//-log(std::abs(tan(.5*theta)));//PPHOT(3)
      double PT  = sqrt(cvec.px()*cvec.px()+cvec.py()*cvec.py());//PPHOT(5)
      double PHI = cvec.phi();//ANGLE(pxpho,pypho)//PPHOT(4)
      
      double x = PT*cos(PHI);
      double y = PT*sin(PHI);
      double z = PT*sinh(ETA);
      double t = PT*cosh(ETA);

      return_vec.setPx(x);
      return_vec.setPy(y);
      return_vec.setPz(z);
      return_vec.setE(t);
      
    } else if (m_smearParamSchema==2){ // Smearing based on CSC sample tuning

      double vE = avec.e();

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

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

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

      int iPos;
      double sigma,sigma1,sigma2;

      iPos=iEnergy-1;
      sigma1=CorrFactorPhotonSigma[iPos][iEta-1]+ratioEta*(CorrFactorPhotonSigma[iPos][iEta]-CorrFactorPhotonSigma[iPos][iEta-1]);
      iPos=iEnergy;
      sigma2=CorrFactorPhotonSigma[iPos][iEta-1]+ratioEta*(CorrFactorPhotonSigma[iPos][iEta]-CorrFactorPhotonSigma[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 : "<<vEnergiesPhoton[iEnergy]<<"   Eta : "<<vEnergiesPhoton[iEnergy-1]<<" < "<<vE<<" < "<<vEnergiesPhoton[iEnergy]<<endreq;
      m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesPhoton[iEnergy]<<"   Eta : "<<EtaPhoton[iEta-1]<<" < "<<vEta<<" < "<<EtaPhoton[iEta]<<endreq;
      m_log<<MSG::DEBUG<<"       calib mean-sigma  "<<alpha<<" "<<sigma<<" -> "<<vrandom<<endreq;

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

    }

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

---