TrackSmearer.cxx

Go to the documentation of this file.

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

#include "AtlfastAlgs/TrackMaker.h"
#include "AtlfastAlgs/MuonMatrixManager.h"
#include "AtlfastAlgs/PionMatrixManager.h"
#include "AtlfastAlgs/ElectronMatrixManager.h"

namespace Atlfast 
{
  using std::abs;
  using std::sin;
  using std::cos;
  using std::tan;
  using std::atan;
  using std::sqrt;
  using std::pow;
  
  //----------------------------------------
  // Constructor
  //----------------------------------------
  
  TrackSmearer::TrackSmearer( int randSeed, 
                              MsgStream& log,
                              std::string muonparamfile, 
                              std::string pionparamfile, 
                              std::string electronparamfile,
                              std::string bremparamfile ) :
    m_log(log), 
    m_muonParamFile(muonparamfile), 
    m_pionParamFile(pionparamfile),
    m_electronParamFile(electronparamfile),
    m_bremParamFile(bremparamfile)
  {

    makeMatrixManagers(randSeed);
    makeChargeService();
    
  }

  void TrackSmearer::makeMatrixManagers(int randSeed)
  {
    m_muonMatrixManager     = new MuonMatrixManager( m_muonParamFile, randSeed, m_log );
    m_pionMatrixManager     = new PionMatrixManager( m_pionParamFile, randSeed, m_log );
    m_electronMatrixManager = new ElectronMatrixManager( m_bremParamFile, m_electronParamFile, randSeed, m_log );
    
    m_log << MSG::INFO << "TrackSmearer constructed" << endreq;

    return;
  }
  
  void TrackSmearer::makeChargeService()
  {
    try
    {
      m_chargeService = new ChargeService() ;
    } 
    catch( std::string errMsg )
    {
      std::cout << "TrackSmearer: Error making ChargeService " << errMsg << std::endl;
    }
    catch( ... )
    {
      std::cout << "TrackSmearer: Unknown Error making ChargeService " << std::endl;
    }
    
    return;
  }

  //------------------------------------------
  // Destructor
  //------------------------------------------
  
  TrackSmearer::~TrackSmearer()
  {
    if (m_muonMatrixManager)     delete m_muonMatrixManager;
    if (m_pionMatrixManager)     delete m_pionMatrixManager;
    if (m_electronMatrixManager) delete m_electronMatrixManager;
    delete m_chargeService;
  }
  
  //------------------------------------------------------------------
  // PUBLIC: smear : takes track trajectory and smears its parameters
  //------------------------------------------------------------------
  Track TrackSmearer::smear( const Track& track ) const
  {

    vector<double> smearVariables ;   // Vector of correlated gaussian variables

    int pdg = abs( track.pdg_id() );
    IMatrixManager* correctManager = getCorrectMatrixManager(pdg);

    TrackTrajectory originalTrajectory = track.trajectory();
    HepSymMatrix sigma; // smear matrix for track
    smearVariables = correctManager->getVariables( originalTrajectory, sigma );
    
    // smear parameters with gaussian variables
    // NOTE: smear variables are delta(d0, z0, phi0, cot(theta0), q/pT)
    TrackParameters parameters = originalTrajectory.parameters();

    double impactParameter = parameters.impactParameter() + smearVariables[0];
    double zPerigee        = parameters.zPerigee() + smearVariables[1];
    Phi phi                = parameters.phi() + smearVariables[2];
    double cotTheta        = parameters.cotTheta() + smearVariables[3];
    double invPtCharge     = parameters.invPtCharge() + smearVariables[4];

    double curvature = originalTrajectory.curvature(); 
    double pT = abs( 1. / invPtCharge );
    Hep3Vector vec;
    vec.setX( pT * cos(phi) );
    vec.setY( pT * sin(phi) );
    vec.setZ( pT * cotTheta );
    
    // Convert results to object
    TrackTrajectory smearedTrajectory( TrackParameters( vec.pseudoRapidity(), 
                                                        phi, pT, impactParameter,
                                                        zPerigee, cotTheta,
                                                        invPtCharge ), 
                                       originalTrajectory.startPoint(),
                                       curvature );

    // sigma has to be converted into the representation used 
    // by Common Tracking, i.e. (d0, z0, phi0, theta0, q/p), 
    // when it is written to CBNTs or AODs.  
    HepSymMatrix sigmaConverted = convertToComTrk( track, sigma );
    
    Track newTrack( originalTrajectory, smearedTrajectory, track.truth(), sigmaConverted );

    return newTrack ;
  }


  IMatrixManager* TrackSmearer::getCorrectMatrixManager( int &pdg ) const 
  {
    IMatrixManager* correctManager = NULL;
    
    if ( pdg == 11 )  correctManager = m_electronMatrixManager;
    if ( pdg == 13 )  correctManager = m_muonMatrixManager;
    if ( pdg == 211 ) correctManager = m_pionMatrixManager;
    if ( correctManager == NULL )  correctManager = m_pionMatrixManager;

    return correctManager;
  }


  HepSymMatrix TrackSmearer::convertToComTrk( const Track& track, const HepSymMatrix& sigma ) const 
  {
    // Method to convert covariance matrix into the representation used 
    // by Common Tracking, i.e. (d0, z0, phi0, theta0, q/p), instead of 
    // internally used (d0, z0, phi0, cot(theta0), q/pT). Needed when  
    // matrix is written to CBNTs or AODs. 
    // Uses track parameters of original (MC) track to determine the 
    // Jacobian matrix.
    
    TrackTrajectory originalTrajectory = track.trajectory();
    TrackParameters parameters = originalTrajectory.parameters();
     
    double cotTheta    = parameters.cotTheta();
    double sinTheta    = sqrt( 1.0 / ( 1.0 + cotTheta * cotTheta ) );
    double invPtCharge = parameters.invPtCharge();
     
    HepSymMatrix conv(5,0);
    
    // USE the fact that from the sigma off-diagonals, the only non-zero 
    // elements are sigma(1,3), sigma(1,5), sigma(3,5), sigma(2,4)
     
    // cov(d0,d0)
    conv(1,1) = sigma(1,1);
    // cov(d0,phi0)
    conv(1,3) = conv(3,1) = sigma(1,3);
    // cov(d0,q/p)
    conv(1,5) = conv(5,1) = sigma(1,5) * sinTheta;
    // cov(phi0,phi0)
    conv(3,3) = sigma(3,3);
    // cov(phi0,q/p)
    conv(3,5) = conv(5,3) = sigma(3,5) * sinTheta;
    
    // cov(z0,z0)
    conv(2,2) = sigma(2,2);
    // cov(z0,theta0)
    conv(2,4) = conv(4,2) = -sigma(2,4) * sinTheta * sinTheta;
    // cov(theta0,theta0)
    conv(4,4) = sigma(4,4) * pow( sinTheta, 4 );
       
    // cov(z0,q/p)
    conv(2,5) = conv(5,2) = conv(2,4) * invPtCharge * sinTheta * cotTheta;
    // cov(theta0,q/p)
    conv(4,5) = conv(5,4) = conv(4,4) * invPtCharge * sinTheta * cotTheta;
    
    // cov(q/p,q/p)
    conv(5,5)  = sigma(5,5) + 
                 sigma(4,4) * pow( invPtCharge * sinTheta * sinTheta * cotTheta, 2 );  
    conv(5,5) *= sinTheta * sinTheta;             
    
    return conv;
  }

  
  double TrackSmearer::getPtResolution( const HepMC::GenParticle& particle )
  {
    int pdg = abs( particle.pdg_id() );
    IMatrixManager* correctManager = getCorrectMatrixManager(pdg);
    
    // Look up vertex
    Hep3Vector vertex = (particle.production_vertex())->position().getV();
    // Get 4-momentum and explicitly convert to 3-momentum
    Hep3Vector threeMomentum = particle.momentum().getV() ;
    // Obtain charge of particle
    double charge = m_chargeService->operator()( &particle );    
    // Magnetic field strength
    MagField *bfield = MagField::Instance();

    TrackTrajectory trTraj( threeMomentum, vertex, charge, bfield->strength() );
    
    // Get smearing matrix
    HepSymMatrix sigma; // smear matrix for track
    std::vector<double> smearVariables = correctManager->getVariables( trTraj, sigma );
    
    // return pT element
    // At some stage, need to develop a system for returning
    // element independent of track parameterisation, eg.
    // AtlfastMatrix->getPtRes()?

    double pt = ( particle.momentum() ).perp();

    // sigma(5,5) is resolution^2 for q/pT, 
    // sigma(pT) = pT^2 * sqrt( sigma(5,5) ) 
    // sigma(pT)/pT gives proportional resolution

    double ptResProp = pt * sqrt( sigma(5,5) );
    
    return ptResProp;
  }

} // end namespace

---