---

MuonSmearer.cxx

Go to the documentation of this file.

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

#include "AtlfastAlgs/MuonSmearer.h"
#include <cmath>

#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Units/SystemOfUnits.h"  
#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RandGauss.h"
#include <iostream>

//-------------------------------------------------------------
// allow call to muon spectrometer fortran resolution function
extern "C" { 
  void resolumu_(float*, float*, float*, float*);
}
//-------------------------------------------------------------
namespace Atlfast {
  using std::abs;
  using std::pair;

  HepLorentzVector MuonSmearer::smear (const HepLorentzVector& avec) {
    // ATLFAST fortran routines   
    //        RESMUO(KEYMUO,KEYFUN,PT,ETA,PHI)
    //        RESOMU (this had been reduced to merely converting radians to degrees)
    // have been converted to C++ here, after restructuring the code.
    //        RESOLUMU (in Atlfast/mu_atlfast/resmuo.F) is left in fortran.

    float eta      = avec.pseudoRapidity();
    float pt       = avec.perp();
    float phi      = avec.phi()*deg;
    float sigmams;
    float sigmamuon;
    float sigmaid;
    float sigmatrack;
    float sigma=0.;
    pair<float, float> sigmapr;
    switch (m_keymuo) {
    case 1:
      while( ((sigmapr=resolms(pt, eta, phi)).first) <-1.) {}
      sigma=sigmapr.first;
      break;
    case 2:
      while( (sigma=resolid1(pt, eta))               <-1.) {}
      break;
    case 3:
      while( ((sigmapr=resolms(pt, eta, phi)).first) <-1.) {}
      sigmams    = sigmapr.first;
      sigmamuon  = sigmapr.second;
      while( ((sigmapr=resolid2(pt, eta)).first)     <-1.) {}
      sigmaid    = sigmapr.first;
      sigmatrack = sigmapr.second;

      sigma=resol1(sigmams, sigmamuon, sigmaid, sigmatrack);
      break;
    case 4:
      while( ((sigmapr=resolms(pt, eta, phi)).first) <-1.) {}
      sigmams    = sigmapr.first;
      sigmamuon  = sigmapr.second;
      while( ((sigmapr=resolid2(pt, eta)).first)     <-1.) {}
      sigmaid    = sigmapr.first;
      sigmatrack = sigmapr.second;

      sigma=resol2(sigmamuon, sigmatrack);
      break;
    case 5:
      sigma=resol3();
      break;
    case 6:
      sigma=resol4();
      
    }

    HepLorentzVector bvec(avec.px()/(1.0+sigma), avec.py()/(1.0+sigma), 
                          avec.pz()/(1.0+sigma), avec.e() /(1.0+sigma));
    return bvec;
  }

  pair<float, float> MuonSmearer::resolms (float pt, float eta, float phi) {
    float cpt=pt;
    float ceta=eta;
    float cphi=phi;
    float resol;
    resolumu_(&cpt, &ceta, &cphi, &resol);
    double aa=randGauss()->fire();
    float sigmamuon    = m_percent*resol;
    float sigmams      = aa*sigmamuon;
    pair<float, float>   sigmapr(sigmams, sigmamuon);

    return sigmapr;
  }

  float MuonSmearer::resolid1 (float pt, float eta) {
    double aa=randGauss()->fire();
    double bb=randGauss()->fire();
    float factor = 1. + pow(fabs(eta),m_magic2);
    float atrack = m_magic1*factor;
    float sigma = atrack*pt*aa + m_magic3*bb;
    sigma = sigma*factor;
      
    return sigma;
  }

  pair<float, float> MuonSmearer::resolid2 (float pt, float eta) {
    double aa=randGauss()->fire();
    double bb=randGauss()->fire();
    double factor = 1. + pow(fabs(eta),m_magic2);
    float atrack = m_magic1*factor;
    float sigmatr    = atrack*pt*aa + m_magic3*bb;
    float sigmatrack = sqrt( (atrack*pt)*(atrack*pt) + m_magic3*m_magic3);
    pair<float, float> sigmapr(sigmatr, sigmatrack);
      
    return sigmapr;
  }

  float MuonSmearer::resol1 ( float sigmams, float sigmamuon, 
                              float sigmaid, float sigmatrack) {
    double wmuon  = 1./(sigmamuon*sigmamuon);
    double wtrack = 1./(sigmatrack*sigmatrack);
    double wtot   = wmuon+wtrack;
    double corr   = (wmuon*(1.0+sigmams)+wtrack*(1.0+sigmaid))/wtot;
    float sigma   = corr-1.0;
      
    return sigma;
  }

  float MuonSmearer::resol2 (float sigmamuon, float sigmatrack) {
    double wmuon  = 1./(sigmamuon*sigmamuon);
    double wtrack = 1./(sigmatrack*sigmatrack);
    double wtot   = wmuon+wtrack;
    double aa     = randGauss()->fire();
    float sigma   = aa/sqrt(wtot);
      
    return sigma;
  }

  float MuonSmearer::resol3 () {
    float sigma  = m_magic4*randGauss()->fire();
      
    return sigma;
  }

  float MuonSmearer::resol4 () {
    float sigmams = m_magic4;
    float sigmaid = m_magic4;
    float wmuon   = 1.0/(sigmams*sigmams);
    float wtrak   = 1.0/(sigmaid*sigmaid);
    float wtot    = wmuon + wtrak;
    double aa     = randGauss()->fire();
    float sigma   = aa/sqrt(wtot);
      
    return sigma;
  }

} // end of namespace bracket

---