---

BremFitter.cxx

Go to the documentation of this file.

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

#include "CLHEP/Random/JamesRandom.h"
#include "AtlfastAlgs/CorrelatedData.h"
#include "AtlfastAlgs/BremEtaBin.h"
#include "AtlfastAlgs/TestValue.h"


namespace Atlfast {
  
  
  BremFitter::BremFitter(int randSeed, MsgStream& log)/*:m_log(log)*/ {
    
    // HOW IT WORKS ---------------
    //  pT is scaled according to Eta of the track
    // There exist 3D Eta, R and PT bins.
    // For a track Eta Bin, there are n bremR values.
    // for each bremR value there are 50 p Values.
    // The bremR value is chosen at random.
    // The p Value within this bremR value is also chosen at random
    // The pT scale is calculated from this random number and the pT bin edges

    m_correlatedData = new CorrelatedData(randSeed);
    m_randomEngine = new HepJamesRandom(randSeed);
    //read files and make variables
    // open file
    m_file = "Atlfast_XKalInput.dat";
    ifstream input;
    input.open( m_file.c_str() );
    if (input) {
      
      log << MSG::INFO <<"BremFitter: File "<<m_file<<" open."<<endreq;
      log << setprecision(3);
      vector<double> empty;
      string skip;
      // read some parameters
      readSingleParameter(log, "RMin",input,m_rMin);
      readSingleParameter(log, "RMax",input,m_rMax);
      readSingleParameter(log, "N RBins",input,m_nRBins);
      readSingleParameter(log, "Eta Extent",input,m_etaExtent);
      readSingleParameter(log, "N Eta Bins",input,m_nEtaBins);
      
      // read BremR values
      double bremR;
      vector<double> bremRStore;
      input >> skip;
      log << MSG::DEBUG << "Reading " << m_nRBins*m_nEtaBins <<" "<< skip << " values" << endreq;
      for (int i=0; i < (m_nRBins*m_nEtaBins ); ++i) {
          input >> bremR;
          log << MSG::DEBUG << i << ": " << bremR << endreq;
          bremRStore.push_back(bremR);
      }

      // make displacements
      readDisplacements(log, input, bremRStore);

      // read more parameters
      readSingleParameter(log, "PMin",input,m_pMin);
      readSingleParameter(log, "PMax",input,m_pMax);
      readSingleParameter(log, "N PBins",input,m_nPBins);
      readSingleParameter(log, "P Scale 1",input,m_p1scal);
      readSingleParameter(log, "P Scale 2",input,m_p2scal);

      double pTScaleHigh;
      double pTScaleLow;
      int pTID=1, etaID=1;
      double pTValue;
      double step = (m_pMax - m_pMin) / (m_nPBins);
      vector<double>::iterator  bremRIter = bremRStore.begin();

      // make pt distribution EtaBin objects
      //====================================
      // there are nEtaBins bins
      for (int iEta=0; iEta < m_nEtaBins; iEta++){
        BremEtaBin* etaBin = new BremEtaBin(randSeed);

        // each bin has nRBins number of RBins
        for (int iR=0; iR < m_nRBins; iR++){

             // skip comments in file
             input >> skip;
             input >> skip;
             input >> skip;
             input >> skip;
             log << MSG::DEBUG << "Reading pT values " << skip << " at radial bin ";
             input >> skip;
             log << MSG::DEBUG << skip << ", eta bin ";
             input >> skip;
             log << MSG::DEBUG << skip << endreq;
             
             //-----
             BremRBin* rBin = new BremRBin(randSeed);
             pTScaleLow =0;
             pTValue = m_pMin;
          // each rBin has nPBins number of PTBins

          for (int iP=0; iP < m_nPBins; iP++){
            input >> pTScaleHigh;
            log << MSG::DEBUG << setw(10) << "("<<pTScaleLow<<" - "<<pTScaleHigh<<
              ", "<<pTID<<", "<< pTValue<<","<<step<<")"<<endreq;
            rBin->addBin(TestValue(++pTID, pTScaleHigh), BremPTScaleBin( pTScaleLow, pTScaleHigh, pTValue, step ) ); 
            pTValue += step;
            pTScaleLow = pTScaleHigh;
          }

          log << MSG::DEBUG << "Adding R Bin "<<etaID<<" with Brem Radius "<<*bremRIter<<endreq;
          etaBin->addBin( TestValue(++etaID, *bremRIter), rBin );
          ++bremRIter;

        }
        m_pTDistEtaBins.push_back(etaBin);
      }
      vector<BremEtaBin*>::const_iterator iter = m_pTDistEtaBins.begin();
      vector<BremEtaBin*>::const_iterator end = m_pTDistEtaBins.end();
      log << MSG::DEBUG << "Electron Brem PT Distribution Bins" << endreq;
      int ieta=0;
      for(;iter != end; ++iter) {
        ++ieta;
        (log) << MSG::DEBUG <<"Eta "<<ieta<<" "<< (**iter) << endreq;
      }
      log << MSG::INFO << "finished reading "<<m_file<<endreq;
    }
    else { log <<MSG::WARNING << "BremFitter: no data file ( "
               <<m_file<<" )!!!!" << endreq;}

    
    // done reading files!!    
  }
  
  TrackTrajectory BremFitter::doBremFit(const TrackTrajectory& traj) const{

    TrackParameters track = traj.parameters();
    int iEta;
    double oldQPT,newQPT, newD0, newPhi;
    double pScale, deltaP;
    //    double pStep;

    // get eta index and find pScale
    iEta = int(m_nEtaBins*abs(track.eta())/ m_etaExtent);
    if (iEta >= m_nEtaBins) iEta = m_nEtaBins-1;

    // -------- Get PSCALE Random --------------
    double random = m_randomEngine->flat();

    pScale = m_pTDistEtaBins[iEta]->calculatePScale(random);

    // scale pT
    oldQPT = track.pT() * traj.signOfCharge();
    newQPT = scalePT(oldQPT, pScale);

    // do a check
    if (newQPT != 0 && oldQPT != 0) deltaP = 1/newQPT - 1/oldQPT;
    else deltaP = 0;

    // generate phi and d0 displacements
    newD0 = scaleDPhi(iEta, random, oldQPT, track.impactParameter(),
                      m_d0Disp, deltaP);
    newPhi = scaleDPhi(iEta, random, oldQPT, track.phi(), m_phiDisp, deltaP);

    // make new track
    double newInvPtCharge = 1 / newQPT;

    TrackParameters newParams(track.eta(),
                              newPhi,
                              abs(newQPT), 
                              newD0,
                              track.zPerigee(),
                              track.cotTheta(), 
                              newInvPtCharge);

    TrackTrajectory bremTrack(newParams, traj.startPoint(), traj.curvature());
    return bremTrack;
  }
  
  double BremFitter::scaleDPhi(int iEta, double randBremR, double pT, double old,
                                const Displacement& axis,
                                double deltaP) const
  {
    double pscale, scale = 20; //GeV scale
    double random = m_correlatedData->generate(1);
    double numerator =  pT*pT + axis.off()*axis.off();
    double denominator = scale*scale + axis.off()*axis.off();
    pscale = pow( numerator/denominator, axis.pow()/2 );
    double mean = axis.average(iEta,randBremR);
    double rms = axis.rms(iEta, randBremR);
    double grad = pscale * (mean + random*rms);

    return old + grad*deltaP;  
  }
  
  double BremFitter::scalePT(double pT, double pscale) const
  { 
    // do some maths (soft pT scaling normalised to 20GeV)
    double scale  = 1 + m_p2scal*log(abs(pT)/20);
    pscale = m_p1scal * (1 + scale*(pscale-1));
    if (pscale < 0) pscale=1;
    return ( pscale * pT );
  }

  void BremFitter::readSingleParameter(MsgStream& s, string name, ifstream& input, double& parameter) {
        string skip;
        input >> skip;
        input >> parameter;
        //m_log << MSG::INFO << "Read " << skip << " with value " << parameter << endreq;
        s << MSG::DEBUG << "Loading " << name <<": Read " << skip << " with value " << parameter << endreq;
      }
  void BremFitter::readSingleParameter(MsgStream& s, string name, ifstream& input, int& parameter) {
        string skip;
        input >> skip;
        input >> parameter;
        //m_log << MSG::INFO << "Read " << skip << " with value " << parameter << endreq;
        s << MSG::DEBUG << "Loading " << name <<": Read " << skip << " with value " << parameter << endreq;
      }

  void BremFitter::readDisplacements(MsgStream& log, ifstream& input, vector<double> bremRStore) {

      double dtpow,pdoff,ftpow,pfoff;
      string skip;

      readSingleParameter(log, "D0 Power",input,dtpow);
      readSingleParameter(log, "D0 Offset",input,pdoff);
      readSingleParameter(log, "Phi Power",input,ftpow);
      readSingleParameter(log, "Phi Offset",input,pfoff);

      // make displacements
      // 4 x n(RBins) (corresponding to 5(eta) lots of the same n(RBins) ) values
      double data;
      vector< vector< map<double, double> > > axes;
      
      // read the 4 vectors(d0 mean, rms and phi mean, rms)
      for (int nVectors = 0; nVectors < 4; ++nVectors) {
      vector<double>::iterator  bremRIter = bremRStore.begin();
        input >> skip;
        log << MSG::DEBUG << "Reading " << m_nRBins << " " << skip << " values" << endreq;
        vector<double> vecData;
        // read parameters
        for (int iR=0; iR < (m_nRBins) ; ++iR){
          input >> data;
          log << MSG::DEBUG << iR << ": " << data << endreq;
          vecData.push_back(data);
        }
        // these n(RBins) parameters are the same for each eta bin, with different bremR values
        // run through each eta bin and make map for each
        vector< map<double,double> > vecDataMap;
        vector<double>::iterator iData, iDataEnd = vecData.end();
        
        // make n(eta) lots of data maps
        for (int iE=0; iE < (m_nEtaBins) ; ++iE){
          map<double,double> dataMap;
          iData = vecData.begin();
          for (; iData != iDataEnd; ++iData) {
            dataMap[*bremRIter] = *iData;
            ++bremRIter;
          }
          vecDataMap.push_back(dataMap);
        }
        
        axes.push_back(vecDataMap);
      }
      // Done all four
      
      m_d0Disp = Displacement(pdoff, dtpow, axes[0], axes[1]);
      m_phiDisp = Displacement(pfoff, ftpow, axes[2], axes[3]);
      //========================================
  }

  
}
//end of namespace

---