File indexing completed on 2025-01-21 01:40:03
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.getParameter<bool>("useLegacyError")),
0016 maxChgOneMIP(conf.getParameter<double>("maxChgOneMIP")),
0017 m_algo(useLegacyError ? Algo::legacy : (maxChgOneMIP < 0 ? Algo::mergeCK : Algo::chargeCK)) {
0018 mLC_P[0] = conf.getParameter<double>("mLC_P0");
0019 mLC_P[1] = conf.getParameter<double>("mLC_P1");
0020 mLC_P[2] = conf.getParameter<double>("mLC_P2");
0021
0022 mHC_P[SiStripDetId::TIB - 3][0] = conf.getParameter<double>("mTIB_P0");
0023 mHC_P[SiStripDetId::TIB - 3][1] = conf.getParameter<double>("mTIB_P1");
0024 mHC_P[SiStripDetId::TID - 3][0] = conf.getParameter<double>("mTID_P0");
0025 mHC_P[SiStripDetId::TID - 3][1] = conf.getParameter<double>("mTID_P1");
0026 mHC_P[SiStripDetId::TOB - 3][0] = conf.getParameter<double>("mTOB_P0");
0027 mHC_P[SiStripDetId::TOB - 3][1] = conf.getParameter<double>("mTOB_P1");
0028 mHC_P[SiStripDetId::TEC - 3][0] = conf.getParameter<double>("mTEC_P0");
0029 mHC_P[SiStripDetId::TEC - 3][1] = conf.getParameter<double>("mTEC_P1");
0030 }
0031
0032 void StripCPEfromTrackAngle::fillPSetDescription(edm::ParameterSetDescription& desc) {
0033 desc.add<bool>("useLegacyError", true);
0034 desc.add<double>("maxChgOneMIP", -6000.);
0035 desc.add<double>("mLC_P0", -.326);
0036 desc.add<double>("mLC_P1", .618);
0037 desc.add<double>("mLC_P2", .300);
0038 desc.add<double>("mTIB_P0", -.742);
0039 desc.add<double>("mTIB_P1", .202);
0040 desc.add<double>("mTID_P0", -1.026);
0041 desc.add<double>("mTID_P1", .253);
0042 desc.add<double>("mTOB_P0", -1.427);
0043 desc.add<double>("mTOB_P1", .433);
0044 desc.add<double>("mTEC_P0", -1.885);
0045 desc.add<double>("mTEC_P1", .471);
0046 }
0047
0048 float StripCPEfromTrackAngle::stripErrorSquared(const unsigned N,
0049 const float uProj,
0050 const SiStripDetId::SubDetector loc) const {
0051 auto fun = [&](float x) -> float { return mLC_P[0] * x * vdt::fast_expf(-x * mLC_P[1]) + mLC_P[2]; };
0052 auto uerr = (N <= 4) ? fun(uProj) : mHC_P[loc - 3][0] + float(N) * mHC_P[loc - 3][1];
0053 return uerr * uerr;
0054 }
0055
0056 float StripCPEfromTrackAngle::legacyStripErrorSquared(const unsigned N, const float uProj) const {
0057 if UNLIKELY ((float(N) - uProj) > 3.5f)
0058 return float(N * N) / 12.f;
0059 else {
0060 static constexpr float P1 = -0.339;
0061 static constexpr float P2 = 0.90;
0062 static constexpr float P3 = 0.279;
0063 const float uerr = P1 * uProj * vdt::fast_expf(-uProj * P2) + P3;
0064 return uerr * uerr;
0065 }
0066 }
0067
0068 void StripCPEfromTrackAngle::localParameters(AClusters const& clusters,
0069 ALocalValues& retValues,
0070 const GeomDetUnit& det,
0071 const LocalTrajectoryParameters& ltp) const {
0072 auto const& par = getAlgoParam(det, ltp);
0073 auto const& p = par.p;
0074 auto loc = par.loc;
0075 auto corr = par.corr;
0076 auto afp = par.afullProjection;
0077
0078 auto fill = [&](unsigned int i, float uerr2) {
0079 const float strip = clusters[i]->barycenter() + corr;
0080 retValues[i].first = p.topology->localPosition(strip, ltp.vector());
0081 retValues[i].second = p.topology->localError(strip, uerr2, ltp.vector());
0082 };
0083
0084 switch (m_algo) {
0085 case Algo::chargeCK:
0086 for (auto i = 0U; i < clusters.size(); ++i) {
0087 auto dQdx = siStripClusterTools::chargePerCM(*clusters[i], ltp, p.invThickness);
0088 auto N = clusters[i]->amplitudes().size();
0089 auto uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0090 fill(i, uerr2);
0091 }
0092 break;
0093 case Algo::legacy:
0094 for (auto i = 0U; i < clusters.size(); ++i) {
0095 auto N = clusters[i]->amplitudes().size();
0096 auto uerr2 = legacyStripErrorSquared(N, afp);
0097 fill(i, uerr2);
0098 }
0099 break;
0100 case Algo::mergeCK:
0101 for (auto i = 0U; i < clusters.size(); ++i) {
0102 auto N = clusters[i]->amplitudes().size();
0103 auto uerr2 = clusters[i]->isMerged() ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0104 fill(i, uerr2);
0105 }
0106 break;
0107 }
0108 }
0109
0110 StripClusterParameterEstimator::LocalValues StripCPEfromTrackAngle::localParameters(const SiStripCluster& cluster,
0111 AlgoParam const& par) const {
0112 auto const& p = par.p;
0113 auto const& ltp = par.ltp;
0114 auto loc = par.loc;
0115 auto corr = par.corr;
0116 auto afp = par.afullProjection;
0117
0118 float uerr2 = 0;
0119
0120 auto N = cluster.amplitudes().size();
0121
0122 switch (m_algo) {
0123 case Algo::chargeCK: {
0124 auto dQdx = siStripClusterTools::chargePerCM(cluster, ltp, p.invThickness);
0125 uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0126 } break;
0127 case Algo::legacy:
0128 uerr2 = legacyStripErrorSquared(N, afp);
0129 break;
0130 case Algo::mergeCK:
0131 uerr2 = cluster.isMerged() ? legacyStripErrorSquared(N, afp) : stripErrorSquared(N, afp, loc);
0132 break;
0133 }
0134
0135 const float strip = cluster.barycenter() + corr;
0136
0137 return std::make_pair(p.topology->localPosition(strip, ltp.vector()),
0138 p.topology->localError(strip, uerr2, ltp.vector()));
0139 }
0140
0141 StripClusterParameterEstimator::LocalValues StripCPEfromTrackAngle::localParameters(
0142 const SiStripCluster& cluster, const GeomDetUnit& det, const LocalTrajectoryParameters& ltp) const {
0143 auto const& par = getAlgoParam(det, ltp);
0144 return localParameters(cluster, par);
0145 }
