StripCPEfromTrackAngle.cc

Go to the documentation of this file.

#include "RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEfromTrackAngle.h"
#include "Geometry/CommonTopologies/interface/StripTopology.h"                                                           
#include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"

#include "vdt/vdtMath.h"

StripCPEfromTrackAngle::StripCPEfromTrackAngle( edm::ParameterSet & conf,
                          const MagneticField& mag,
                          const TrackerGeometry& geom,
                          const SiStripLorentzAngle& lorentz,
                          const SiStripBackPlaneCorrection& backPlaneCorrection,
                          const SiStripConfObject& confObj,
                          const SiStripLatency& latency)
  : StripCPE(conf, mag, geom, lorentz, backPlaneCorrection, confObj, latency )
  , useLegacyError(conf.existsAs<bool>("useLegacyError") ? conf.getParameter<bool>("useLegacyError") : true)
  , maxChgOneMIP(conf.existsAs<float>("maxChgOneMIP") ? conf.getParameter<double>("maxChgOneMIP") : -6000.),
    m_algo(useLegacyError ? Algo::legacy : ( maxChgOneMIP<0 ? Algo::mergeCK : Algo::chargeCK))
  {
    mLC_P[0] = conf.existsAs<double>("mLC_P0") ? conf.getParameter<double>("mLC_P0") : -.326;
    mLC_P[1] = conf.existsAs<double>("mLC_P1") ? conf.getParameter<double>("mLC_P1") :  .618;
    mLC_P[2] = conf.existsAs<double>("mLC_P2") ? conf.getParameter<double>("mLC_P2") :  .300;

    mHC_P[SiStripDetId::TIB - 3][0] = conf.existsAs<double>("mTIB_P0") ? conf.getParameter<double>("mTIB_P0") : -.742  ;
    mHC_P[SiStripDetId::TIB - 3][1] = conf.existsAs<double>("mTIB_P1") ? conf.getParameter<double>("mTIB_P1") :  .202  ;
    mHC_P[SiStripDetId::TID - 3][0] = conf.existsAs<double>("mTID_P0") ? conf.getParameter<double>("mTID_P0") : -1.026 ;
    mHC_P[SiStripDetId::TID - 3][1] = conf.existsAs<double>("mTID_P1") ? conf.getParameter<double>("mTID_P1") :  .253  ;
    mHC_P[SiStripDetId::TOB - 3][0] = conf.existsAs<double>("mTOB_P0") ? conf.getParameter<double>("mTOB_P0") : -1.427 ;
    mHC_P[SiStripDetId::TOB - 3][1] = conf.existsAs<double>("mTOB_P1") ? conf.getParameter<double>("mTOB_P1") :  .433  ;
    mHC_P[SiStripDetId::TEC - 3][0] = conf.existsAs<double>("mTEC_P0") ? conf.getParameter<double>("mTEC_P0") : -1.885 ;
    mHC_P[SiStripDetId::TEC - 3][1] = conf.existsAs<double>("mTEC_P1") ? conf.getParameter<double>("mTEC_P1") :  .471  ;
  }

float StripCPEfromTrackAngle::stripErrorSquared(const unsigned N, const float uProj, const SiStripDetId::SubDetector loc ) const {
  auto fun = [&] (float x)  -> float { return mLC_P[0]*x*vdt::fast_expf(-x*mLC_P[1])+mLC_P[2];};
  auto uerr = (N <= 4) ?  fun(uProj) :  mHC_P[loc-3][0]+float(N)*mHC_P[loc-3][1];
  return uerr*uerr;
}


float StripCPEfromTrackAngle::legacyStripErrorSquared(const unsigned N, const float uProj) const {
  if UNLIKELY( (float(N)-uProj) > 3.5f )
    return float(N*N)/12.f;
  else {
    static constexpr float P1=-0.339;
    static constexpr float P2=0.90;
    static constexpr float P3=0.279;
    const float uerr = P1*uProj*vdt::fast_expf(-uProj*P2)+P3;
    return uerr*uerr;
  }
}



void StripCPEfromTrackAngle::localParameters(AClusters const & clusters, ALocalValues & retValues, const GeomDetUnit& det, const LocalTrajectoryParameters & ltp) const {


  auto const & par = getAlgoParam(det,ltp);
  auto const & p = par.p;
  auto loc = par.loc;
  auto corr = par.corr;
  auto afp = par.afullProjection;

  auto fill = [&](unsigned int i, float uerr2) {
    const float strip = clusters[i]->barycenter() + corr;
    retValues[i].first  = p.topology->localPosition(strip, ltp.vector());
    retValues[i].second = p.topology->localError(strip, uerr2, ltp.vector());
  };

  
  switch (m_algo) {
  case Algo::chargeCK :
    for (auto i=0U; i< clusters.size(); ++i) {
      auto dQdx = siStripClusterTools::chargePerCM(*clusters[i], ltp, p.invThickness);
      auto N = clusters[i]->amplitudes().size();
      auto uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N, afp,loc );
      fill(i, uerr2);
    }
    break;
  case Algo::legacy :
    for (auto i=0U; i< clusters.size(); ++i)  {
      auto N = clusters[i]->amplitudes().size();
      auto uerr2 = legacyStripErrorSquared(N,afp);
      fill(i, uerr2);
    }
    break;
  case Algo::mergeCK :
    for (auto i=0U; i< clusters.size(); ++i) {
      auto N = clusters[i]->amplitudes().size();
      auto uerr2 = clusters[i]->isMerged() ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N,afp,loc );
      fill(i, uerr2);
    }
    break;
  }


}

StripClusterParameterEstimator::LocalValues
StripCPEfromTrackAngle::localParameters( const SiStripCluster& cluster, AlgoParam const & par) const {
  auto const & p = par.p;
  auto const & ltp = par.ltp;
  auto loc = par.loc;
  auto corr = par.corr;
  auto afp = par.afullProjection;

  float uerr2=0;
  
  auto N = cluster.amplitudes().size();

  switch (m_algo) {
  case Algo::chargeCK :
    {
      auto dQdx = siStripClusterTools::chargePerCM(cluster, ltp, p.invThickness);
      uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N, afp,loc );
    }
    break;
  case Algo::legacy :
    uerr2 = legacyStripErrorSquared(N,afp);
    break;
  case Algo::mergeCK :
    uerr2 = cluster.isMerged() ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N, afp,loc );
    break;
  }
  
  const float strip = cluster.barycenter() + corr;
  
  return std::make_pair( p.topology->localPosition(strip, ltp.vector()),
             p.topology->localError(strip, uerr2, ltp.vector()) );
}

  

StripClusterParameterEstimator::LocalValues 
StripCPEfromTrackAngle::localParameters( const SiStripCluster& cluster, const GeomDetUnit& det, const LocalTrajectoryParameters& ltp) const {
  
  auto const & par = getAlgoParam(det,ltp);
  return localParameters(cluster,par);
}