---

PionMatrixManager.cxx

Go to the documentation of this file.

#include "AtlfastAlgs/PionMatrixManager.h"
#include <iostream>

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

  //-----------------------------------------------------------
  // PUBLIC: Constructor
  //-----------------------------------------------------------
  
  PionMatrixManager::PionMatrixManager(int randSeed, MsgStream log):m_randSeed(randSeed){
    log<< MSG::INFO << "Constructing PionMatrixManager"<< endreq;
    m_file = "Atlfast_PionResParam.dat";
    m_correlatedData = new CorrelatedData(randSeed);
    m_log = &log;
    this->initialise();
    *m_log << MSG::INFO << "Constructed PionMatrixManager"<< endreq;
  }
  
  
  //-----------------------------------------------------------
  // PUBLIC: Destructor
  //-----------------------------------------------------------
  
  PionMatrixManager::~PionMatrixManager(){
    
    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 PionMatrixManager::initialise(){
    // open file
    ifstream input;
    input.open( m_file.c_str() );
    
    if (input) {
      
      *m_log << MSG::INFO <<"PionMatrixManager: File "<<m_file<<" open."<<endreq;
      
      double pTMin, etaMin, etaMax, rtBoundary;
      int nEtaBoundaries;
      
      // read some parameters
      input >> pTMin;
      input >> etaMin;
      input >> etaMax;
      input >> nEtaBoundaries;
      input >> m_nRTBins;
      m_nEtaBins = nEtaBoundaries - 1;
      
      // construct vector<binboundaries>
      double etaStep = (etaMax - etaMin)/double(nEtaBoundaries);
      for (int i=0; i<m_nEtaBins+1; 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:
      // Core: C0(d0, z0, phi, cotTheta, 1/qPT), C1(d0, z0, phi, cotTheta, 1/qPT) etc.
      // Tails: C0(d0, z0, phi, cotTheta, 1/qPT), C1(d0, z0, phi, cotTheta, 1/qPT) etc.
      // Fractions: C0(d0, z0, phi, cotTheta, 1/qPT), C1(d0, z0, phi, cotTheta, 1/qPT) etc.
      // Correlations: C0(1,2,3,4), C1(1,2,3,4) etc.

      //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> >  coreC0(nEtaBoundaries, empty); 
        vector< vector<double> >  coreC1(nEtaBoundaries, empty);
        vector< vector<double> >  coreC2(nEtaBoundaries, empty);
        vector< vector<double> >  tailsC0(nEtaBoundaries, empty);
        vector< vector<double> >  tailsC1(nEtaBoundaries, empty);
        vector< vector<double> >  tailsC2(nEtaBoundaries, empty);
        vector< vector<double> >  fractionsC0(nEtaBoundaries, empty);
        vector< vector<double> >  fractionsC1(nEtaBoundaries, empty);
        vector< vector<double> >  correlationsC0(nEtaBoundaries, empty);
        vector< vector<double> >  correlationsC1(nEtaBoundaries, empty);
        vector< vector<double> >  correlationsC3(nEtaBoundaries, empty);
        
        //read eta bin number of values for each C and parameter
        // read core data
        fillVector(input, coreC0, 5);
        fillVector(input, coreC1, 5);
        fillVector(input, coreC2, 5);
        // read tail data
        fillVector(input, tailsC0, 5);
        fillVector(input, tailsC1, 5);
        fillVector(input, tailsC2, 5);
        // read fractions data
        fillVector(input, fractionsC0, 5);
        fillVector(input, fractionsC1, 5);
        // read correlations data
        fillVector(input, correlationsC0, 4);
        fillVector(input, correlationsC1, 4);
        fillVector(input, correlationsC3, 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; ++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*> core;
          vector<ParameterResolutions*> tails;
          vector<ParameterResolutions*> fractions;
          vector<ParameterResolutions*> correlations;
          for (int param=0; param < 5; param++) {
            vector<BinID> coreCoeffBins;
            coreCoeffBins.push_back(BinID(0, coreC0[i][param], coreC0[i+1][param]) );
            coreCoeffBins.push_back(BinID(0, coreC1[i][param], coreC1[i+1][param]) );
            coreCoeffBins.push_back(BinID(0, coreC2[i][param], coreC2[i+1][param]) );
            core.push_back(new ParameterResolutions(coreCoeffBins,etaLow, etaHigh ));

            vector<BinID> tailsCoeffBins;
            tailsCoeffBins.push_back(BinID(0, tailsC0[i][param], tailsC0[i+1][param]) );
            tailsCoeffBins.push_back(BinID(0, tailsC1[i][param], tailsC1[i+1][param]) );
            tailsCoeffBins.push_back(BinID(0, tailsC2[i][param], tailsC2[i+1][param]) );
            tails.push_back(new ParameterResolutions(tailsCoeffBins,etaLow, etaHigh));

            vector<BinID> fractionsCoeffBins;
            fractionsCoeffBins.push_back(BinID(0, fractionsC0[i][param], fractionsC0[i+1][param]) );
            fractionsCoeffBins.push_back(BinID(0, fractionsC1[i][param], fractionsC1[i+1][param]) );
            fractions.push_back(new ParameterResolutions(fractionsCoeffBins,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, 0, 0) );
            correlationsCoeffBins.push_back(BinID(0, correlationsC3[i][param], correlationsC3[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
          PionBinData* binData =  new PionBinData(rTEtaBin,
                                                  core,
                                                  tails,
                                                  fractions,
                                                  correlations,
                                                  m_randSeed);
          
          // enter bin data into map
          m_binData[rTEtaBin] = binData;
          
          //increment bin ID
          iBin++;
          
        }//end of eta bin loop
      } //end of rt bin loop
      
    }
    else { *m_log << MSG::INFO << "PionMatrixManager: no data file ( "
                  <<m_file<<" )!!!!" << endreq;}
    
    makeHeader();
    
  }
  
  void PionMatrixManager::makeHeader(){
    HeaderPrinter hp("Atlfast PionTrack Smearer:", *m_log);
    hp.add("Total number of Bins     ", int(m_binData.size()) );
    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 PionMatrixManager::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 matrix corresponding
  //                     to a given track trajectory
  //-----------------------------------------------------------
  
  vector<double> PionMatrixManager::getVariables(const TrackTrajectory& track,
                                                 HepMatrix& returnSigma) const{
    HepMatrix 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* PionMatrixManager::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
    cout << "WARNING: PionMatrixManager - No bin; rT "<< rT <<", eta "<<eta<<endl;
    return (m_binData.begin())->second;
  }
  
  
}// namespace

---