Back to home page

Project CMSSW displayed by LXR

 
 

    


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; }