---

LeptonMatrixManager.cxx

Go to the documentation of this file.

#include "AtlfastAlgs/LeptonMatrixManager.h"

#include "AtlfastAlgs/LeptonBinData.h"
#include "AtlfastAlgs/CorrelatedData.h"

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

  //-----------------------------------------------------------
  // PUBLIC: Constructor
  //-----------------------------------------------------------
  LeptonMatrixManager::LeptonMatrixManager(string config, int randSeed,
                                       MsgStream& log):
    m_config(config), m_log(log){
    m_file = "Atlfast_id_resol_mu_";
    m_file.append(config);
    m_file.append(".dat");
    m_correlatedData = new CorrelatedData(randSeed);
    this->initialise();
    m_log << MSG::INFO << "Constructed LeptonMatrixManager with configuration "
           << m_config <<endreq;
  }
  
  
  //-----------------------------------------------------------
  // PUBLIC: Destructor
  //-----------------------------------------------------------
  
  LeptonMatrixManager::~LeptonMatrixManager(){
    map<BinID, IBinData*>::iterator iter = m_binData.begin();
    map<BinID, IBinData*>::iterator end = m_binData.end();
      for (; iter!=end;++iter)
      {
      delete (iter->second);  // Is this right(m_binData=vector<*>)????
      }

  }
  
  //------------------------------------------------------------
  // PRIVATE: initialise : read file and construct data bins
  //------------------------------------------------------------
  
  void LeptonMatrixManager::initialise(){
    // open file
    ifstream input;
    input.open( m_file.c_str() );
    if (input) {
      m_log << MSG::INFO 
             << "LeptonMatrixManager: file " 
             << m_file
             << " open." 
             << endreq;
      
      double data;
      int nEtaVals;
      int nPTVals;
      
      // prepare containers
      vector<double> etaBins;
      vector<double> pTBins;
      
      // the are four correlations between the 5 variables
      vector< pair<double,double> > correlOrder; 
      correlOrder.push_back(pair<double,double>(3,5)); // phi / pt
      correlOrder.push_back(pair<double,double>(1,5)); // d0 / pt
      correlOrder.push_back(pair<double,double>(1,3)); // d0 / phi
      correlOrder.push_back(pair<double,double>(2,4)); // z0 / cotTheta
      
      
      // read eta and phi values and store them
      input >> nEtaVals;
      for (int i=0; i < nEtaVals; i++) {
        input >> data;
        etaBins.push_back(data);
      }
      
      input >> nPTVals;
      for (int i=0; i < nPTVals; i++) {
        input >> data;
        pTBins.push_back(data);
      }
      
      //========================================
      // make 9 resolution parameters for matrix
      HepMatrix etaPMatrix(nEtaVals, nPTVals,0);
      //stores nine parameters for each eta/pt bin:
      vector<HepMatrix> resParameters(9, etaPMatrix); 
      
      // read sigmas and correlations
      for (int iEta=0; iEta < nEtaVals; iEta++) {
        for (int iPT=0; iPT < nPTVals; iPT++) {
          //sigmas (diagonals)
          for (int i = 0; i < 5; i++) {
            input >> data;
            resParameters[i][iEta][iPT] = data*data;
          }
          // correlations
          vector< pair<double,double> >::iterator iOrder = correlOrder.begin();
          for (int i = 5; i < 9; i++) {
            input >> data;
            resParameters[i][iEta][iPT] = 
              data * sqrt(resParameters[(iOrder->first)-1][iEta][iPT]) * 
              sqrt(resParameters[(iOrder->second)-1][iEta][iPT]);
            iOrder++;
          }
        }
      }// end of make ===========================
      
      
      // make bins
      HepMatrix twoSquare(2,2,0);
      vector<HepMatrix> interpol(9, twoSquare);
      double iBin = 0;    
      
      for (int iEta = 0; iEta < (nEtaVals-1); iEta++) {
        for (int iPT = 0; iPT < (nPTVals-1); iPT++) {
          
          for (int i=0; i<9; i++) {
            interpol[i] = resParameters[i].sub(iEta+1,iEta+2,iPT+1, iPT+2);
          }
          
          BinID id(iBin, etaBins[iEta], etaBins[iEta+1],
                   pTBins[iPT], pTBins[iPT+1]);
          // make bin
          IBinData* binData = new LeptonBinData(id, interpol);      
          // enter bin data into map
          m_binData[id] = binData;
          iBin++;
        }
      }
      HeaderPrinter hp("Atlfast LeptonTrack Smearer:", m_log);
      hp.add("Data file                ", m_file);        
      hp.add("configuration            ", m_config);        
      hp.add("Number of Eta Bins       ", nEtaVals-1);
      hp.add("Eta Min                  ", etaBins.front());
      hp.add("Eta Max                  ", etaBins.back());
      hp.add("Number of Pt  Bins       ", nPTVals-1);
      hp.add("pT Min                   ", pTBins.front());
      hp.add("pT Max                   ", pTBins.back());
      hp.print();
      
    }else{ 
      m_log << MSG::WARNING 
             << "LeptonMatrixManager: no data file ( "
             <<m_file<<" )!!!!" 
             << endreq;
    }
    
  }
  
  //-----------------------------------------------------------
  // PUBLIC: getVariables : returns the correlated smear variables
  //                        from the sigma matrix corresponding
  //                        to a given track trajectory
  //-----------------------------------------------------------
  vector<double> LeptonMatrixManager::getVariables(const TrackTrajectory& track, 
                                                 HepMatrix& returnSigma) const{
    vector<double> variables;
    HepMatrix sigma;
    IBinData* binData = getBinData(track);
    sigma = binData->getMatrix(track);
    // do Dcorset and Dcorgen to get smear parameters

    variables = m_correlatedData->generate(m_correlatedData->root(sigma));
    scale(sigma);
    returnSigma = sigma;
    scale(variables);
    return variables;
  }
  
  //----------------------------------------------------------
  // Private: Scale variables and matrix
  //----------------------------------------------------------
  void LeptonMatrixManager::scale(HepMatrix& sigma) const{
    sigma(1,1) /= 1e8;
    sigma(2,2) /= 1e8;
    sigma(3,3) /= 1e6;
    sigma(4,4) /= 1e6;
    sigma(5,5) /= 1e6;
    sigma(1,3) = sigma(3,1) /= 1e7;
    sigma(1,5) = sigma(5,1) /= 1e7;
    sigma(2,4) = sigma(4,2) /= 1e7;
    sigma(3,5) = sigma(5,3) /= 1e6;
  }
  
  void LeptonMatrixManager::scale(vector<double>& variables) const{
    variables[0] /= 1e4;
    variables[1] /= 1e4;
    variables[2] /= 1e3;
    variables[3] /= 1e3;
    variables[4] /= 1e3;
  }
  
  
  //-----------------------------------------------------------
  // Private: getBinID : returns the BinID of bin corresponding
  //                     to a given track trajectory
  //-----------------------------------------------------------
  IBinData* LeptonMatrixManager::getBinData(const TrackTrajectory& track) const{

    TrackParameters traj = track.parameters();
    vector<double> etaPT;
    double eta = abs( traj.eta() ) ;
    double pT= abs( traj.pT() );

    // validity check
    double etaLow = ( (m_binData.begin())->first ).low(0);
    double etaHigh = ( (m_binData.rbegin())->first ).high(0);
    double pTLow = ( (m_binData.begin())->first ).low(1);
    double pTHigh = ( (m_binData.rbegin())->first ).high(1);
    if ( eta < etaLow ) eta = etaLow;
    if ( eta > etaHigh) eta = etaHigh;  
    if ( pT < pTLow ) pT = pTLow;
    if ( pT > pTHigh) pT = pTHigh;  

    // find BinID
    etaPT.push_back( eta );
    etaPT.push_back( pT );
    
    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(etaPT)) return binIter->second;
    } 
    // OOPS! couldn't fin bin
    cout << "WARNING: LeptonMatrixManager - No bin; eta "<< eta <<", pT "<<pT<<endl;
    return (m_binData.begin())->second;
  }
  
  
}// namespace

---