File indexing completed on 2024-04-06 12:26:29
0001 #include "RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEgeometric.h"
0002 #include "RecoLocalTracker/SiStripRecHitConverter/interface/CrosstalkInversion.h"
0003 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0004 #include <numeric>
0005
0006 StripClusterParameterEstimator::LocalValues StripCPEgeometric::localParameters(
0007 const SiStripCluster& cluster, const GeomDetUnit& det, const LocalTrajectoryParameters& ltp) const {
0008 StripCPE::Param const& p = param(det);
0009
0010 const LocalPoint& pos = ltp.position();
0011 LocalVector track = ltp.momentum();
0012 track *= (track.z() < 0) ? fabs(p.thickness / track.z())
0013 : (track.z() > 0) ? -fabs(p.thickness / track.z())
0014 : p.maxLength / track.mag();
0015
0016 const float fullProjection = p.coveredStrips(track + p.drift, pos);
0017 stats_t<float> projection;
0018 {
0019 const float absProj = fabs(fullProjection);
0020 const float minProj = 2 * p.thickness * tan_diffusion_angle / p.topology->localPitch(pos);
0021 const float projection_rel_err2 = thickness_rel_err2 + p.pitch_rel_err2;
0022 projection = stats_t<float>::from_relative_uncertainty2(std::max(absProj, minProj), projection_rel_err2);
0023 }
0024
0025 const std::vector<stats_t<float> > Q =
0026 reco::InverseCrosstalkMatrix::unfold(cluster, xtalk1[static_cast<int>(p.moduleGeom)]);
0027 const stats_t<float> strip = cluster.firstStrip() + offset_from_firstStrip(Q, projection);
0028
0029 const float corrected =
0030 strip() - 0.5 * (1 - p.backplanecorrection) * fullProjection + 0.5 * p.coveredStrips(track, ltp.position());
0031 const float error2 = std::max(strip.error2(), minimum_uncertainty_squared);
0032
0033 return std::make_pair(p.topology->localPosition(corrected, ltp.vector()),
0034 p.topology->localError(corrected, error2, ltp.vector()));
0035 }
0036
0037 stats_t<float> StripCPEgeometric::offset_from_firstStrip(const std::vector<stats_t<float> >& Q,
0038 const stats_t<float>& proj) const {
0039 WrappedCluster wc(Q);
0040 if (useNPlusOne(wc, proj))
0041 wc.addSuppressedEdgeStrip();
0042 else
0043 while (useNMinusOne(wc, proj))
0044 wc.dropSmallerEdgeStrip();
0045
0046 if (proj() < wc.N - 2)
0047 return stats_t<float>(wc.middle(), pow(wc.N - proj(), 2) / 12.);
0048 if (wc.deformed())
0049 return stats_t<float>(wc.centroid()(), 1 / 12.);
0050 if (proj > wc.maxProjection())
0051 return stats_t<float>(wc.centroid()(), 1 / 12.);
0052
0053 if (ambiguousSize(wc, proj)) {
0054 const stats_t<float> probably = geometric_position(wc, proj);
0055 wc.dropSmallerEdgeStrip();
0056 const stats_t<float> maybe = geometric_position(wc, proj);
0057 return stats_t<float>(probably(),
0058 std::max(probably.error2(), float(maybe.error2() + pow(probably() - maybe(), 2) / 12)));
0059 }
0060 return geometric_position(wc, proj);
0061 }
0062
0063 stats_t<float> StripCPEgeometric::geometric_position(const StripCPEgeometric::WrappedCluster& wc,
0064 const stats_t<float>& proj) const {
0065 const stats_t<float> x = wc.middle() + 0.5 * proj * wc.eta();
0066 return wc.N == 1 ? stats_t<float>(x(), pow(1 - 0.82 * proj(), 2) / 12)
0067 : stats_t<float>(x(), scaling_squared * x.error2());
0068 }
0069
0070 inline bool StripCPEgeometric::useNPlusOne(const WrappedCluster& wc, const stats_t<float>& proj) const {
0071 return wc.maxProjection() < proj && proj() < wc.N + 1;
0072 }
0073
0074 inline bool StripCPEgeometric::useNMinusOne(const WrappedCluster& wc, const stats_t<float>& proj) const {
0075 if (proj() > wc.N - 1)
0076 return false;
0077 if (wc.smallerEdgeStrip() < 0)
0078 return true;
0079 if (proj() < wc.N - 3)
0080 return false;
0081 if (proj() < wc.N - 2)
0082 return true;
0083 if (wc.eta().sigmaFrom(0) < 3)
0084 return false;
0085
0086 WrappedCluster wcTest(wc);
0087 wcTest.dropSmallerEdgeStrip();
0088 if (proj >= wcTest.maxProjection())
0089 return false;
0090 if (wc.sign() * wc.eta()() > 1. / (wc.N - 1))
0091 return true;
0092
0093 return wc.smallerEdgeStrip().sigmaFrom(0) < noise_threshold;
0094 }
0095
0096 inline bool StripCPEgeometric::ambiguousSize(const WrappedCluster& wc, const stats_t<float>& proj) const {
0097 return proj() < wc.N - 1 && wc.smallerEdgeStrip()() > 0 && wc.smallerEdgeStrip().sigmaFrom(0) < maybe_noise_threshold;
0098 }
0099
0100 StripCPEgeometric::WrappedCluster::WrappedCluster(const std::vector<stats_t<float> >& Q)
0101 : N(Q.size() - 2), clusterFirst(Q.begin() + 1), first(clusterFirst) {}
0102
0103 inline void StripCPEgeometric::WrappedCluster::addSuppressedEdgeStrip() {
0104 if (eta().sigmaFrom(0) < 1) {
0105 first--;
0106 N += 2;
0107 } else if (*first > last()) {
0108 first--;
0109 N += 1;
0110 } else {
0111 N += 1;
0112 }
0113 }
0114
0115 inline void StripCPEgeometric::WrappedCluster::dropSmallerEdgeStrip() {
0116 if (*first < last()) {
0117 first++;
0118 N -= 1;
0119 } else if (last() < *first) {
0120 N -= 1;
0121 } else {
0122 first++;
0123 N -= 2;
0124 }
0125 }
0126
0127 inline float StripCPEgeometric::WrappedCluster::middle() const { return (first - clusterFirst) + N / 2.; }
0128
0129 inline stats_t<float> StripCPEgeometric::WrappedCluster::centroid() const {
0130 stats_t<float> sumXQ(0);
0131 for (std::vector<stats_t<float> >::const_iterator i = first; i < first + N; i++)
0132 sumXQ += (i - clusterFirst) * (*i);
0133 return sumXQ / sumQ() + 0.5;
0134 }
0135
0136 inline stats_t<float> StripCPEgeometric::WrappedCluster::sumQ() const {
0137 return accumulate(first, first + N, stats_t<float>(0));
0138 }
0139
0140 inline stats_t<float> StripCPEgeometric::WrappedCluster::eta() const { return (last() - *first) / sumQ(); }
0141
0142 inline bool StripCPEgeometric::WrappedCluster::deformed() const {
0143 return N > 2 && std::max((*first)(), last()()) > accumulate(first + 1, first + N - 1, stats_t<float>(0))() / (N - 2);
0144 }
0145
0146 inline stats_t<float> StripCPEgeometric::WrappedCluster::maxProjection() const {
0147 return N * (1 + sign() * eta()).inverse();
0148 }
0149
0150 inline stats_t<float> StripCPEgeometric::WrappedCluster::smallerEdgeStrip() const { return std::min(*first, last()); }
0151
0152 inline int StripCPEgeometric::WrappedCluster::sign() const { return (*first < last()) ? 1 : -1; }