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>

#include "PathResolver/PathResolver.h"

//-------------------------------------------------------------
namespace Atlfast {
  using std::abs;
  using std::pair;
  

  // MuonSmearer methods

  MuonSmearer::MuonSmearer(const int aseed, const int lumi, const int keymuo,
                           std::string &resolutionFile, MsgStream& log): 
    ISmearer(), DefaultSmearer(aseed), m_lumi(lumi), m_keymuo(keymuo)
  {    
    makeResolutionCalculator(resolutionFile);
   
    // Make a (default) TrackSmearer
    int randSeed = 12345;
    m_tracksmearer = new TrackSmearer(randSeed,log);
  }
  
  MuonSmearer::MuonSmearer(const int aseed, const int lumi, const int keymuo,
                           std::string &resolutionFile): 
    ISmearer(), DefaultSmearer(aseed), m_lumi(lumi), m_keymuo(keymuo)
  {    
    makeResolutionCalculator(resolutionFile);
  }

  MuonSmearer::~MuonSmearer()
  {
    delete m_muonResCalculator;
    delete m_tracksmearer;
  }
  
  void MuonSmearer::makeResolutionCalculator(std::string &resolutionFile){
    
    std::string filename =
      PathResolver::find_file(resolutionFile, "DATAPATH");
    std::cout << "MuonSmearer getting file " << filename << std::endl;
    
    if (filename!="")
      m_muonResCalculator = new MuonResolutionCalculator(filename.c_str());
    else
      std::cout << "Could not open muon resolution file!!" << std::endl;

    return;
    
  }
  
  HepLorentzVector MuonSmearer::smear(const HepLorentzVector& avec) {

    // ID smearing requires a GenParticle so this method is called for
    // spectrometer standalone muons only
    
    HepLorentzVector bvec(0,0,0,0);
    pair<double, double> sigmapr;
    if(m_smearParamSchema==1){
      while( ((sigmapr=resolMS(avec)).first) <-1.) {}
      double sigma = sigmapr.first;
      bvec = avec/(1.0+sigma);
    }
    return bvec; // bvec is zeroes if m_smearParamSchema!=1
  }

  HepLorentzVector MuonSmearer::smear(const HepMC::GenParticle& particle){
    
    HepLorentzVector avec = particle.momentum();
    
    double sigmaMS;
    double sigmamuon;
    double sigmaID;
    double sigmatrack;
    double sigma=0.;
    HepLorentzVector bvec(0,0,0,0);
    pair<double, double> sigmapr;

    if(m_smearParamSchema==1){
      switch (m_keymuo) {
      case 1:
        while( ((sigmapr=resolMS(avec)).first) <-1.) {}
        sigma = sigmapr.first;
        break;
      case 2:
        while( ((sigmapr=resolID(particle)).first)<-1.) {}
        sigma = sigmapr.first;
        break;
      case 3:
        while( ((sigmapr=resolMS(avec)).first) <-1.) {}
        sigmaMS    = sigmapr.first;
        sigmamuon  = sigmapr.second;
        while( ((sigmapr=resolID(particle)).first)<-1.) {}
        sigmaID    = sigmapr.first;
        sigmatrack = sigmapr.second;
        
        sigma=resol(sigmaMS, sigmamuon, sigmaID, sigmatrack);
        break;
      default:
        std::cout << "Invalid value for MuonSmearKey: " << m_keymuo 
                  << ", returning unsmeared muon vector" << std::endl;
        return avec;
        break;      
      }
      
      bvec = avec/(1.0+sigma);
      
      //std::cout << "Unsmeared muon: " << avec << ", smeared muon: " << bvec << std::endl;
      
    } 
    
    return bvec; // bvec is zeroes if m_smearParamSchema!=1
  }
  
  pair<double, double> MuonSmearer::resolMS(const HepLorentzVector& avec) {
    
    if (!m_muonResCalculator)
      std::cout << "No MuonResolutionCalculator, resolution set to 0.0" << std::endl;
    
    double sigmamuon = m_muonResCalculator ?
      m_muonResCalculator->calculateResolution(avec) : 0.;
    if (!sigmamuon){
      std::cout << "Zero spectrometer resolution!!!" << std::endl;
      pair<double, double>  sigmapr(0,0);
      return sigmapr;
    }
    
    double aa=randGauss()->fire();
    double sigmams   = aa*sigmamuon;
    pair<double, double> sigmapr(sigmams, sigmamuon);
    
    return sigmapr;
  }
  
  pair<double, double> MuonSmearer::resolID(const HepMC::GenParticle& particle) {

    if (!m_tracksmearer){
      std::cout << "No TrackSmearer in MuonSmearer, resolution set to 0.0" << std::endl;
      std::cout << "To get inner detector resolutions, you must instantiate"
                << " MuonSmearer with a MsgStream" << std::endl;
    }

    double sigmatrack = m_tracksmearer ? 
      m_tracksmearer->getPtResolution(particle) : 0.;

    if (!sigmatrack){
      std::cout << "Zero inner detector resolution!!!" << std::endl;
      pair<double, double>  sigmapr(0,0);
      return sigmapr;
    }

    double aa=randGauss()->fire();
    double sigmatr = aa*sigmatrack;
    pair<double, double> sigmapr(sigmatr, sigmatrack);
      
    return sigmapr;

  }

  double MuonSmearer::resol ( double sigmams, double sigmamuon, 
                              double sigmaid, double 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;
    double sigma   = corr-1.0;
    
    return sigma;
  }
  
  //cct: implement the setSmearParameters method
  int MuonSmearer::setSmearParameters (const std::vector<double>& smearValues){
    m_smearParams=smearValues;
    return 0;
  }
  //cct: implement the setSmearParamSchema method
  int MuonSmearer::setSmearParamSchema ( const int smearSchema){
    m_smearParamSchema=smearSchema;
    return 0;
  }

  
  
} // end of namespace bracket

---