MuonMatrixManager.cxx

Go to the documentation of this file.

#include "AtlfastAlgs/MuonMatrixManager.h"
#include "PathResolver/PathResolver.h"
#include <iostream>

namespace Atlfast
{
  using std::abs;
  using std::ifstream;

  //-----------------------------------------------------------
  // PUBLIC: Constructor
  //-----------------------------------------------------------
  MuonMatrixManager::MuonMatrixManager( string paramFile, int randSeed, MsgStream log ) : 
       m_file( paramFile ), 
       m_randSeed( randSeed )
  {
    m_log = &log;
    *m_log << MSG::INFO 
           << "Constructing MuonMatrixManager with parameter file " 
           << paramFile 
           << endreq;
    m_file = PathResolver::find_file( m_file, "DATAPATH" );
    m_correlatedData = new CorrelatedData(randSeed);
    this->initialise();
    *m_log << MSG::INFO << "Constructed MuonMatrixManager" << endreq;
  }
  
  
  //-----------------------------------------------------------
  // PUBLIC: Destructor
  //-----------------------------------------------------------
  MuonMatrixManager::~MuonMatrixManager()
  {
    delete m_correlatedData;

    map< BinID, IBinData* >::iterator iter = m_binData.begin();
    map< BinID, IBinData* >::iterator end  = m_binData.end();
    for ( ; iter != end; ++iter )  delete ( iter->second ); 
  }
  
  //------------------------------------------------------------
  // PRIVATE: initialise : read file and construct data bins
  //------------------------------------------------------------
  void MuonMatrixManager::initialise()
  {
    // open file
    ifstream input;
    input.open( m_file.c_str() );
    
    if (input) 
    {
      *m_log << MSG::INFO << "MuonMatrixManager: File " << m_file << " open." << endreq;
      
      double pTMin, etaMin, etaMax, rTBoundary;
      
      // read some parameters
      input >> pTMin;
      input >> etaMin;
      input >> etaMax;
      input >> m_nEtaBins;
      input >> m_nRTBins;
      
      // construct vector<binboundaries>
      double etaStep = ( etaMax - etaMin )/ double( m_nEtaBins );
      for ( int i = 0; i < m_nEtaBins; i++ )
      { 
        m_etaBoundaries.push_back( etaMin + etaStep/2. + double(i)*etaStep );
      }
      for ( int i = 0; i < m_nRTBins + 1; i++ )
      {
        input >> rTBoundary;
        m_rTBoundaries.push_back( rTBoundary );
      }
      
      // parameters are stored in rT bin blocks----
      // with parameters in format:
      // Gaussian sigmas: 
      //              C0(d0, z0, phi0, cot(theta0), q/pT), 
      //              C1(d0, z0, phi0, cot(theta0), q/pT), etc.
      // correlation coefficients rho: 
      //              C0(rho(d0,phi0), rho(d0,q/pT), rho(phi0,q/pT), rho(z0,cot(theta0))), 
      //              C1(rho(d0,phi0), rho(d0,q/pT), rho(phi0,q/pT), rho(z0,cot(theta0))), etc.
      //
      // coefficients C0,C1,C2,C3,C4 define pT-dependence of the respective quantities
      // according to f(pT) = C0 + C1/sqrt(pT) + C2/pT + C3/pT/sqrt(pT) + C4/pT^2

      // start bin id ints at zero
      int iBin = 0; 
      
      for ( int rt = 0; rt < m_nRTBins; rt++ )
      {
        // make vectors to hold all resolution parameters for this rt bin
        vector<double> empty;
        vector< vector<double> >  sigmaC0( m_nEtaBins, empty ); 
        vector< vector<double> >  sigmaC1( m_nEtaBins, empty );
        vector< vector<double> >  sigmaC2( m_nEtaBins, empty );
        vector< vector<double> >  sigmaC3( m_nEtaBins, empty );
        vector< vector<double> >  sigmaC4( m_nEtaBins, empty );
        vector< vector<double> >  correlationsC0( m_nEtaBins, empty );
        vector< vector<double> >  correlationsC1( m_nEtaBins, empty );
        vector< vector<double> >  correlationsC2( m_nEtaBins, empty );
        vector< vector<double> >  correlationsC3( m_nEtaBins, empty );
        vector< vector<double> >  correlationsC4( m_nEtaBins, empty );
        
        // read eta bin number of values for each C and parameter
        // read Gaussian sigmas
        fillVector( input, sigmaC0, 5 );
        fillVector( input, sigmaC1, 5 );
        fillVector( input, sigmaC2, 5 );
        fillVector( input, sigmaC3, 5 );
        fillVector( input, sigmaC4, 5 );
        // read correlations data
        fillVector( input, correlationsC0, 4 );
        fillVector( input, correlationsC1, 4 );
        fillVector( input, correlationsC2, 4 );
        fillVector( input, correlationsC3, 4 );
        fillVector( input, correlationsC4, 4 );
        
        // DATA READING FINISHED FOR THIS RT BIN
        // got all values, now make rT/eta bins and resolution objects
        for ( int i = 0; i < m_nEtaBins - 1; ++i ) 
        {         
          double etaLow  = m_etaBoundaries[i];
          double etaHigh = m_etaBoundaries[i+1];
          
          // make bin id with rT and eta boundaries
          BinID rTEtaBin( iBin, m_rTBoundaries[rt], m_rTBoundaries[rt+1], etaLow, etaHigh );
          
          //make parameter resolutions for each parameter
          vector<ParameterResolutions*> sigmas;
          vector<ParameterResolutions*> correlations;
          for ( int param = 0; param < 5; param++ ) 
          {
             vector<BinID> sigmaCoeffBins;
             sigmaCoeffBins.push_back( BinID( 0, sigmaC0[i][param], sigmaC0[i+1][param] ) );
             sigmaCoeffBins.push_back( BinID( 0, sigmaC1[i][param], sigmaC1[i+1][param] ) );
             sigmaCoeffBins.push_back( BinID( 0, sigmaC2[i][param], sigmaC2[i+1][param] ) );
             sigmaCoeffBins.push_back( BinID( 0, sigmaC3[i][param], sigmaC3[i+1][param] ) );
             sigmaCoeffBins.push_back( BinID( 0, sigmaC4[i][param], sigmaC4[i+1][param] ) );
             sigmas.push_back( new ParameterResolutions( sigmaCoeffBins, etaLow, etaHigh ) );
          }
          for ( int param = 0; param < 4; param++ ) 
          {
            vector<BinID> correlationsCoeffBins;
            correlationsCoeffBins.push_back( BinID( 0, correlationsC0[i][param], correlationsC0[i+1][param]) );
            correlationsCoeffBins.push_back( BinID( 0, correlationsC1[i][param], correlationsC1[i+1][param]) );
            correlationsCoeffBins.push_back( BinID( 0, correlationsC2[i][param], correlationsC2[i+1][param]) );
            correlationsCoeffBins.push_back( BinID( 0, correlationsC3[i][param], correlationsC3[i+1][param]) );
            correlationsCoeffBins.push_back( BinID( 0, correlationsC4[i][param], correlationsC4[i+1][param]) );
            correlations.push_back( new ParameterResolutions( correlationsCoeffBins, etaLow, etaHigh ) );
          }
          
          // Now make a bin data with edges rT and eta, containing the 
          // correlation data
          MuonBinData* binData = new MuonBinData( rTEtaBin, sigmas, correlations );
          
          // enter bin data into map
          m_binData[rTEtaBin] = binData;
          
          // increment bin ID
          iBin++;
          
        } // end of eta bin loop
      } // end of rT bin loop

      // close file
      input.close();
      
    }
    else  // no input file 
    { 
      *m_log << MSG::INFO 
             << "MuonMatrixManager: no data file ( " << m_file << " )!!!!" 
             << endreq;
    }
    
    makeHeader();
    
  }
  
  void MuonMatrixManager::makeHeader()
  {
    HeaderPrinter hp( "Atlfast MuonTrack Smearer:", *m_log );
    hp.add( "Total number of Bins     ", m_nRTBins * m_nEtaBins );
    hp.add( "rT Bin Min               ", m_rTBoundaries.front() );        
    hp.add( "rT Bin Max               ", m_rTBoundaries.back() );
    hp.add( "Number of rT Bins        ", m_nRTBins );
    hp.add( "Eta Bin Min              ", m_etaBoundaries.front() );        
    hp.add( "Eta Bin Max              ", m_etaBoundaries.back() );
    hp.add( "Number of Eta Bins       ", m_nEtaBins );
    hp.print();
  }
  
  //-----------------------------------------------------------
  // PRIVATE: fillVector
  //             reads n x M parameters into member variable
  //-----------------------------------------------------------
  void MuonMatrixManager::fillVector( ifstream& input, 
                                        vector< vector<double> >& myVector, 
                                        int M=5 )
  {
    double res;
    for ( int param = 0; param < M; param++ ) 
    { 
      vector< vector<double> >::iterator binIter = myVector.begin();
      for ( ; binIter != myVector.end(); binIter++ )
      {
        input >> res;
        binIter->push_back(res);
      } 
    }
  }
  
 
  //-----------------------------------------------------------
  // PUBLIC: getMatrix : returns the Sigma (=covariance) matrix 
  //                  corresponding to a given track trajectory
  //-----------------------------------------------------------
  vector<double> MuonMatrixManager::getVariables( const TrackTrajectory& track,
                                                    HepSymMatrix& returnSigma ) const
  {
    HepSymMatrix Sigma;
    vector<double> variables;

    IBinData* binData = getBinData(track);
    Sigma = binData->getMatrix(track);
    
    // do Dcorset and Dcorgen to get smear parameters
    variables = m_correlatedData->generate( m_correlatedData->root(Sigma) );
    
    returnSigma = Sigma;
    return variables;
  }
  
  
  //-----------------------------------------------------------
  // Private: getBinData : returns the IBinData of bin 
  //                 corresponding to a given track trajectory
  //-----------------------------------------------------------
  IBinData* MuonMatrixManager::getBinData( const TrackTrajectory& traj ) const
  {
    TrackParameters track = traj.parameters();

    vector<double> rTEta;
    double rT  = abs( traj.radius() );
    double eta = abs( track.eta() );

    // validity check
    double rTLow   = ( (m_binData.begin())->first ).low(0);
    double rTHigh  = ( (m_binData.rbegin())->first ).high(0);
    double etaLow  = ( (m_binData.begin())->first ).low(1);
    double etaHigh = ( (m_binData.rbegin())->first ).high(1);

    if ( rT < rTLow ) rT = rTLow;
    if ( rT > rTHigh) rT = rTHigh;  
    if ( eta < etaLow ) eta = etaLow;
    if ( eta > etaHigh) eta = etaHigh;  

    // find BinID
    rTEta.push_back(rT);
    rTEta.push_back(eta);

    map<BinID, IBinData*>::const_iterator binIter = m_binData.begin();
    map<BinID, IBinData*>::const_iterator binEnd  = m_binData.end();

    for ( ; binIter != binEnd; ++binIter )
    {
      if ( binIter->first.isInBin(rTEta) ) 
      {
        return binIter->second;
      }
    }
    
    // OOPS! couldn't fin bin
    *m_log << MSG::WARNING 
           << "WARNING: MuonMatrixManager - No bin; rT " << rT << ", eta " << eta 
           << endreq;

    return ( m_binData.begin() )->second;
  }
  
}// namespace

---