File indexing completed on 2023-03-17 11:19:47
0001 #include "RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEfromTrackAngle.h"
0002 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0003 #include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
0004
0005 #include "vdt/vdtMath.h"
0006
0007 StripCPEfromTrackAngle::StripCPEfromTrackAngle(edm::ParameterSet& conf,
0008 const MagneticField& mag,
0009 const TrackerGeometry& geom,
0010 const SiStripLorentzAngle& lorentz,
0011 const SiStripBackPlaneCorrection& backPlaneCorrection,
0012 const SiStripConfObject& confObj,
0013 const SiStripLatency& latency)
0014 : StripCPE(conf, mag, geom, lorentz, backPlaneCorrection, confObj, latency),
0015 useLegacyError(conf.existsAs<bool>("useLegacyError") ? conf.getParameter<bool>("useLegacyError") : true),
0016 maxChgOneMIP(conf.existsAs<float>("maxChgOneMIP") ? conf.getParameter<double>("maxChgOneMIP") : -6000.),
0017 m_algo(useLegacyError ? Algo::legacy : (maxChgOneMIP < 0 ? Algo::mergeCK : Algo::chargeCK)) {
0018 mLC_P[0] = conf.existsAs<double>("mLC_P0") ? conf.getParameter<double>("mLC_P0") : -.326;
0019 mLC_P[1] = conf.existsAs<double>("mLC_P1") ? conf.getParameter<double>("mLC_P1") : .618;
0020 mLC_P[2] = conf.existsAs<double>("mLC_P2") ? conf.getParameter<double>("mLC_P2") : .300;
0021
0022 mHC_P[SiStripDetId::TIB - 3][0] = conf.existsAs<double>("mTIB_P0") ? conf.getParameter<double>("mTIB_P0") : -.742;
0023 mHC_P[SiStripDetId::TIB - 3][1] = conf.existsAs<double>("mTIB_P1") ? conf.getParameter<double>("mTIB_P1") : .202;
0024 mHC_P[SiStripDetId::TID - 3][0] = conf.existsAs<double>("mTID_P0") ? conf.getParameter<double>("mTID_P0") : -1.026;
0025 mHC_P[SiStripDetId::TID - 3][1] = conf.existsAs<double>("mTID_P1") ? conf.getParameter<double>("mTID_P1") : .253;
0026 mHC_P[SiStripDetId::TOB - 3][0] = conf.existsAs<double>("mTOB_P0") ? conf.getParameter<double>("mTOB_P0") : -1.427;
0027 mHC_P[SiStripDetId::TOB - 3][1] = conf.existsAs<double>("mTOB_P1") ? conf.getParameter<double>("mTOB_P1") : .433;
0028 mHC_P[SiStripDetId::TEC - 3][0] = conf.existsAs<double>("mTEC_P0") ? conf.getParameter<double>("mTEC_P0") : -1.885;
0029 mHC_P[SiStripDetId::TEC - 3][1] = conf.existsAs<double>("mTEC_P1") ? conf.getParameter<double>("mTEC_P1") : .471;
0030 }
0031
0032 float StripCPEfromTrackAngle::stripErrorSquared(const unsigned N,
0033 const float uProj,
0034 const SiStripDetId::SubDetector loc) const {
0035 auto fun = [&](float x) -> float { return mLC_P[0] * x * vdt::fast_expf(-x * mLC_P[1]) + mLC_P[2]; };
0036 auto uerr = (N <= 4) ? fun(uProj) : mHC_P[loc - 3][0] + float(N) * mHC_P[loc - 3][1];
0037 return uerr * uerr;
0038 }
0039
0040 float StripCPEfromTrackAngle::legacyStripErrorSquared(const unsigned N, const float uProj) const {
0041 if UNLIKELY ((float(N) - uProj) > 3.5f)
0042 return float(N * N) / 12.f;
0043 else {
0044 static constexpr float P1 = -0.339;
0045 static constexpr float P2 = 0.90;
0046 static constexpr float P3 = 0.279;
0047 const float uerr = P1 * uProj * vdt::fast_expf(-uProj * P2) + P3;
0048 return uerr * uerr;
0049 }
0050 }
0051
0052 void StripCPEfromTrackAngle::localParameters(AClusters const& clusters,
0053 ALocalValues& retValues,
0054 const GeomDetUnit& det,
0055 const LocalTrajectoryParameters& ltp) const {
0056 auto const& par = getAlgoParam(det, ltp);
0057 auto const& p = par.p;
0058 auto loc = par.loc;
0059 auto corr = par.corr;
0060 auto afp = par.afullProjection;
0061
0062 auto fill = [&](unsigned int i, float uerr2) {
0063 const float strip = clusters[i]->barycenter() + corr;
0064 retValues[i].first = p.topology->localPosition(strip, ltp.vector());
0065 retValues[i].second = p.topology->localError(strip, uerr2, ltp.vector());
0066 };
0067
0068 switch (m_algo) {
0069 case Algo::chargeCK:
0070 for (auto i = 0U; i < clusters.size(); ++i) {
0071 auto dQdx = siStripClusterTools::chargePerCM(*clusters[i], ltp, p.invThickness);
0072 auto N = clusters[i]->amplitudes().size();
0073 auto uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0074 fill(i, uerr2);
0075 }
0076 break;
0077 case Algo::legacy:
0078 for (auto i = 0U; i < clusters.size(); ++i) {
0079 auto N = clusters[i]->amplitudes().size();
0080 auto uerr2 = legacyStripErrorSquared(N, afp);
0081 fill(i, uerr2);
0082 }
0083 break;
0084 case Algo::mergeCK:
0085 for (auto i = 0U; i < clusters.size(); ++i) {
0086 auto N = clusters[i]->amplitudes().size();
0087 auto uerr2 = clusters[i]->isMerged() ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0088 fill(i, uerr2);
0089 }
0090 break;
0091 }
0092 }
0093
0094 StripClusterParameterEstimator::LocalValues StripCPEfromTrackAngle::localParameters(const SiStripCluster& cluster,
0095 AlgoParam const& par) const {
0096 auto const& p = par.p;
0097 auto const& ltp = par.ltp;
0098 auto loc = par.loc;
0099 auto corr = par.corr;
0100 auto afp = par.afullProjection;
0101
0102 float uerr2 = 0;
0103
0104 auto N = cluster.amplitudes().size();
0105
0106 switch (m_algo) {
0107 case Algo::chargeCK: {
0108 auto dQdx = siStripClusterTools::chargePerCM(cluster, ltp, p.invThickness);
0109 uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0110 } break;
0111 case Algo::legacy:
0112 uerr2 = legacyStripErrorSquared(N, afp);
0113 break;
0114 case Algo::mergeCK:
0115 uerr2 = cluster.isMerged() ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0116 break;
0117 }
0118
0119 const float strip = cluster.barycenter() + corr;
0120
0121 return std::make_pair(p.topology->localPosition(strip, ltp.vector()),
0122 p.topology->localError(strip, uerr2, ltp.vector()));
0123 }
0124
0125 StripClusterParameterEstimator::LocalValues StripCPEfromTrackAngle::localParameters(
0126 const SiStripCluster& cluster, const GeomDetUnit& det, const LocalTrajectoryParameters& ltp) const {
0127 auto const& par = getAlgoParam(det, ltp);
0128 return localParameters(cluster, par);
0129 }
