---

LeptonBinData.cxx

Go to the documentation of this file.

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

namespace Atlfast {
  using std::vector;

  //-----------------------------------------------
  // PUBLIC - Constructor
  //-----------------------------------------------
  //  
  //  LeptonBinData::LeptonBinData(
  //                       const BinID& id,
  //                       const vector<HepMatrix>& interpol
  //                       ):
  //    m_id(id),m_interpol(interpol){
  //  }
  
  //-----------------------------------------------
  // PUBLIC - void calculate
  // calculate matrix from data
  //-----------------------------------------------
  
  HepMatrix LeptonBinData::getMatrix(const TrackTrajectory& traj) const{    
    TrackParameters track = traj.parameters();
    HepMatrix sigma(5,5,0);
    vector<double> sigmaParams, par1, par2;
    
    double eta = track.eta();
    double etaLow = m_id.low(0);
    double etaHigh = m_id.high(0);
    double etaRatio = (abs(eta) - etaLow ) / (etaHigh - etaLow);
    double pT = track.pT();
    double pTLow = m_id.low(1);
    double pTHigh = m_id.high(1);

    //    cout << "LeptonBinData:  Getting Matrix for eta: " <<eta<<" and pT: " << pT << endl;
    //    cout << "eta ratio " << etaRatio << " with pt low, high " << pTLow<<","<<pTHigh<<endl;
    //    cout << "eta low, high " << etaLow <<","<< etaHigh;

    vector<HepMatrix>::const_iterator iter = m_interpol.begin();
    vector<HepMatrix>::const_iterator end = m_interpol.end();
    for ( ;iter != end; ++iter){
      //  cout << "interpol = " << *iter << endl;
      par1.push_back( 
                     (*iter)[0][0] + 
                     ( (*iter)[1][0] - (*iter)[0][0] )*etaRatio
                     );
      
      par2.push_back( 
                     (*iter)[0][1] + 
                     ( (*iter)[1][1] - (*iter)[0][1] )*etaRatio
                     );
      
      // cout << "etaRatio interpol parameters are " << par1.back() << " and " << par2.back() << endl;

      double a,b;
      b = ( par2.back() - par1.back() ) * 
        ( (pTHigh*pTLow) ) *
        ( (pTHigh*pTLow) ) /
        ( (pTLow*pTLow)-(pTHigh*pTHigh) ); 
      
      a = ( par2.back()*(pTHigh*pTHigh) -
            par1.back()*(pTLow*pTLow)
            ) / 
        (   (pTHigh*pTHigh)-(pTLow*pTLow)   );
      
      sigmaParams.push_back( a + b/(pT*pT) );
      //cout << "a = " << a<<"   b = "<<b<<endl;
      //cout << "a + b/pT^2 = " << sigmaParams.back() << endl;
    }
    
    for (int i = 0;i < 5; i++){ 
      sigma[i][i] = sigmaParams[i];
    }
    sigma(1,3) = sigma(3,1) = sigmaParams[7];
    sigma(1,5) = sigma(5,1) = sigmaParams[6];
    sigma(2,4) = sigma(4,2) = sigmaParams[8];
    sigma(3,5) = sigma(5,3) = sigmaParams[5];
    return sigma;
  }
  
  
}//end of namespace

---