SiStripRecHitConverter/interface/SiStripRecHitMatcher.h

0001 #ifndef RECOLOCALTRACKER_SISTRIPCLUSTERIZER_SISTRIPRECHITMATCH_H
0002 #define RECOLOCALTRACKER_SISTRIPCLUSTERIZER_SISTRIPRECHITMATCH_H
0003
0004 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0005 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0006 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0008 #include "DataFormats/DetId/interface/DetId.h"
0009 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0010 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0011
0012 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0013
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include "TrackingTools/TransientTrackingRecHit/interface/HelpertRecHit2DLocalPos.h"
0017
0018 class GluedGeomDet;
0019
0020 #include <cfloat>
0021 #include <functional>
0022 #include <memory>
0023 #include <vector>
0024
0025 class SiStripRecHitMatcher {
0026 public:
0027   // may become a template argument
0028   typedef SiStripMatchedRecHit2DCollectionNew::FastFiller CollectorMatched;
0029
0030   typedef SiStripRecHit2DCollectionNew::DetSet::const_iterator RecHitIterator;
0031   typedef std::vector<const SiStripRecHit2D*> SimpleHitCollection;
0032   typedef SimpleHitCollection::const_iterator SimpleHitIterator;
0033
0034   typedef std::function<void(SiStripMatchedRecHit2D const&)> Collector;
0035
0036   typedef std::pair<LocalPoint, LocalPoint> StripPosition;
0037
0038   SiStripRecHitMatcher(const edm::ParameterSet& conf);
0039   SiStripRecHitMatcher(const double theScale);
0040
0041   bool preFilter() const { return preFilter_; }
0042
0043   inline static float sigmaPitch(LocalPoint const& pos, LocalError const& err, const StripTopology& topol) {
0044     MeasurementError error = topol.measurementError(pos, err);
0045     auto pitch = topol.localPitch(pos);
0046     return error.uu() * pitch * pitch;
0047   }
0048
0049   // optimized matching iteration (the algo is the same, just recoded)
0050   template <typename MonoIterator, typename StereoIterator, typename CollectorHelper>
0051   void doubleMatch(MonoIterator monoRHiter,
0052                    MonoIterator monoRHend,
0053                    StereoIterator seconditer,
0054                    StereoIterator seconditerend,
0055                    const GluedGeomDet* gluedDet,
0056                    LocalVector trdir,
0057                    CollectorHelper& collectorHelper) const;
0058
0059   //match a single hit
0060   std::unique_ptr<SiStripMatchedRecHit2D> match(const SiStripRecHit2D* monoRH,
0061                                                 const SiStripRecHit2D* stereoRH,
0062                                                 const GluedGeomDet* gluedDet,
0063                                                 LocalVector trackdirection,
0064                                                 bool force) const;
0065
0066   // this is the one used by the RecHitConverter
0067   void match(const SiStripRecHit2D* monoRH,
0068              RecHitIterator begin,
0069              RecHitIterator end,
0070              CollectorMatched& collector,
0071              const GluedGeomDet* gluedDet,
0072              LocalVector trackdirection) const;
0073
0074   void match(const SiStripRecHit2D* monoRH,
0075              SimpleHitIterator begin,
0076              SimpleHitIterator end,
0077              CollectorMatched& collector,
0078              const GluedGeomDet* gluedDet,
0079              LocalVector trackdirection) const;
0080
0081   // project strip coordinates on Glueddet
0082
0083   StripPosition project(const GeomDetUnit* det,
0084                         const GluedGeomDet* glueddet,
0085                         StripPosition strip,
0086                         LocalVector trackdirection) const;
0087
0088   // needed by the obsolete version still in use on some architectures
0089   void match(const SiStripRecHit2D* monoRH,
0090              SimpleHitIterator begin,
0091              SimpleHitIterator end,
0092              std::vector<std::unique_ptr<SiStripMatchedRecHit2D>>& collector,
0093              const GluedGeomDet* gluedDet,
0094              LocalVector trackdirection) const;
0095
0096   void match(const SiStripRecHit2D* monoRH,
0097              SimpleHitIterator begin,
0098              SimpleHitIterator end,
0099              std::vector<SiStripMatchedRecHit2D*>& collector,
0100              const GluedGeomDet* gluedDet,
0101              LocalVector trackdirection) const;
0102
0103   /// the actual implementation
0104   void match(const SiStripRecHit2D* monoRH,
0105              SimpleHitIterator begin,
0106              SimpleHitIterator end,
0107              Collector& collector,
0108              const GluedGeomDet* gluedDet,
0109              LocalVector trackdirection) const;
0110
0111   float scale_;
0112   bool preFilter_ = false;
0113 };
0114
0115 #include "DataFormats/GeometryVector/interface/Basic3DVector.h"
0116 #if defined(USE_SSEVECT) || defined(USE_EXTVECT)
0117 #define DOUBLE_MATCH
0118 #endif
0119
0120 #ifdef DOUBLE_MATCH
0121
0122 #if defined(USE_SSEVECT)
0123 using mathSSE::Rot3F;
0124 using mathSSE::Vec2D;
0125 using mathSSE::Vec3D;
0126 using mathSSE::Vec3F;
0127 #endif
0128
0129 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0130 #include "TrackingTools/TransientTrackingRecHit/interface/HelpertRecHit2DLocalPos.h"
0131
0132 namespace matcherDetails {
0133
0134   struct StereoInfo {
0135     Vec2D c1vec;
0136     const SiStripRecHit2D* secondHit;
0137     float sigmap22;
0138     double m10, m11;
0139   };
0140
0141 }  // namespace matcherDetails
0142
0143 template <typename MonoIterator, typename StereoIterator, typename CollectorHelper>
0144 void SiStripRecHitMatcher::doubleMatch(MonoIterator monoRHiter,
0145                                        MonoIterator monoRHend,
0146                                        StereoIterator seconditer,
0147                                        StereoIterator seconditerend,
0148                                        const GluedGeomDet* gluedDet,
0149                                        LocalVector trdir,
0150                                        CollectorHelper& collectorHelper) const {
0151   using matcherDetails::StereoInfo;
0152
0153   typedef GloballyPositioned<float> ToGlobal;
0154   typedef typename GloballyPositioned<float>::ToLocal ToLocal;
0155
0156   // hits in both mono and stero
0157   // match
0158   bool notk = trdir.mag2() < float(FLT_MIN);
0159   // FIXME we shall find a faster approximation for trdir: not useful to compute it each time for each strip
0160
0161   // stripdet = mono
0162   // partnerstripdet = stereo
0163   const GeomDetUnit* stripdet = gluedDet->monoDet();
0164   const GeomDetUnit* partnerstripdet = gluedDet->stereoDet();
0165   const StripTopology& topol = (const StripTopology&)stripdet->topology();
0166   const StripTopology& partnertopol = (const StripTopology&)partnerstripdet->topology();
0167
0168   // toGlobal is fast,  toLocal is slow
0169   ToGlobal const& stripDetTrans = stripdet->surface();
0170   ToGlobal const& partnerStripDetTrans = partnerstripdet->surface();
0171   ToLocal gluedDetInvTrans(gluedDet->surface());
0172
0173   std::vector<StereoInfo> cache(std::distance(seconditer, seconditerend));
0174   //iterate on stereo rechits
0175   // fill cache with relevant info
0176   int cacheSize = 0;
0177   for (; seconditer != seconditerend; ++seconditer) {
0178     const SiStripRecHit2D& secondHit = CollectorHelper::stereoHit(seconditer);
0179
0180     float sigmap22 = sigmaPitch(secondHit.localPosition(), secondHit.localPositionError(), partnertopol);
0181     // assert (sigmap22>=0);
0182
0183     auto STEREOpointX = partnertopol.measurementPosition(secondHit.localPositionFast()).x();
0184     MeasurementPoint STEREOpointini(STEREOpointX, -0.5f);
0185     MeasurementPoint STEREOpointend(STEREOpointX, 0.5f);
0186
0187     LocalPoint locp1 = partnertopol.localPosition(STEREOpointini);
0188     LocalPoint locp2 = partnertopol.localPosition(STEREOpointend);
0189
0190     GlobalPoint globalpointini = partnerStripDetTrans.toGlobal(locp1);
0191     GlobalPoint globalpointend = partnerStripDetTrans.toGlobal(locp2);
0192
0193     // position of the initial and final point of the strip in glued local coordinates
0194     LocalPoint positiononGluedini = gluedDetInvTrans.toLocal(globalpointini);
0195     LocalPoint positiononGluedend = gluedDetInvTrans.toLocal(globalpointend);
0196
0197     // in case of no track hypothesis assume a track from the origin through the center of the strip
0198     if (notk) {
0199       const LocalPoint& lcenterofstrip = secondHit.localPositionFast();
0200       GlobalPoint gcenterofstrip = partnerStripDetTrans.toGlobal(lcenterofstrip);
0201       GlobalVector gtrackdirection = gcenterofstrip - GlobalPoint(0, 0, 0);
0202       trdir = gluedDetInvTrans.toLocal(gtrackdirection);
0203     }
0204
0205     Vec3F offset = trdir.basicVector().v * positiononGluedini.z() / trdir.z();
0206
0207     Vec3F ret1 = positiononGluedini.basicVector().v - offset;
0208     Vec3F ret2 = positiononGluedend.basicVector().v - offset;
0209
0210     double m10 = -(ret2[1] - ret1[1]);
0211     double m11 = ret2[0] - ret1[0];
0212     double dd = m11 * ret1[1] + m10 * ret1[0];
0213
0214     Vec2D c1vec{dd, dd};
0215
0216     // store
0217     StereoInfo info = {c1vec, &secondHit, sigmap22, m10, m11};
0218     cache[cacheSize++] = info;
0219   }
0220
0221   for (; monoRHiter != monoRHend; ++monoRHiter) {
0222     SiStripRecHit2D const& monoRH = CollectorHelper::monoHit(monoRHiter);
0223
0224     // position of the initial and final point of the strip (RPHI cluster) in local strip coordinates
0225     auto RPHIpointX = topol.measurementPosition(monoRH.localPositionFast()).x();
0226     MeasurementPoint RPHIpointini(RPHIpointX, -0.5f);
0227     MeasurementPoint RPHIpointend(RPHIpointX, 0.5f);
0228
0229     // position of the initial and final point of the strip in local coordinates (mono det)
0230     //StripPosition stripmono=StripPosition(topol.localPosition(RPHIpointini),topol.localPosition(RPHIpointend));
0231     LocalPoint locp1o = topol.localPosition(RPHIpointini);
0232     LocalPoint locp2o = topol.localPosition(RPHIpointend);
0233
0234     // in case of no track hypothesis assume a track from the origin through the center of the strip
0235     if (notk) {
0236       const LocalPoint& lcenterofstrip = monoRH.localPositionFast();
0237       GlobalPoint gcenterofstrip = stripDetTrans.toGlobal(lcenterofstrip);
0238       GlobalVector gtrackdirection = gcenterofstrip - GlobalPoint(0, 0, 0);
0239       trdir = gluedDetInvTrans.toLocal(gtrackdirection);
0240     }
0241
0242     //project mono hit on glued det
0243     //StripPosition projectedstripmono=project(stripdet,gluedDet,stripmono,trackdirection);
0244
0245     GlobalPoint globalpointini = stripDetTrans.toGlobal(locp1o);
0246     GlobalPoint globalpointend = stripDetTrans.toGlobal(locp2o);
0247
0248     // position of the initial and final point of the strip in glued local coordinates
0249     LocalPoint positiononGluedini = gluedDetInvTrans.toLocal(globalpointini);
0250     LocalPoint positiononGluedend = gluedDetInvTrans.toLocal(globalpointend);
0251
0252     Vec3F offset = trdir.basicVector().v * positiononGluedini.z() / trdir.z();
0253
0254     Vec3F projini = positiononGluedini.basicVector().v - offset;
0255     Vec3F projend = positiononGluedend.basicVector().v - offset;
0256
0257     // ret1o = ret1o + (trdir * (ret1o.getSimd(2) / trdirz));
0258     // ret2o = ret2o + (trdir * (ret2o.getSimd(2) / trdirz));
0259
0260     double m00 = -(projend[1] - projini[1]);  //-(projectedstripmono.second.y()-projectedstripmono.first.y());
0261     double m01 = (projend[0] - projini[0]);   // (projectedstripmono.second.x()-projectedstripmono.first.x());
0262     double c0 =
0263         m01 * projini[1] + m00 * projini[0];  //m01*projectedstripmono.first.y()   + m00*projectedstripmono.first.x();
0264
0265     Vec2D c0vec{c0, c0};
0266     Vec2D minv00{-m01, m00};
0267
0268     //error calculation (the part that depends on mono RH only)
0269     double c1 = -m00;
0270     double s1 = -m01;
0271     double l1 = 1. / (c1 * c1 + s1 * s1);
0272
0273     float sigmap12 = sigmaPitch(monoRH.localPosition(), monoRH.localPositionError(), topol);
0274     // float sigmap12 = monoRH.sigmaPitch();
0275     // assert(sigmap12>=0);
0276
0277     //float code
0278     float fc1(c1), fs1(s1);
0279     Vec3F scc1{fs1, fc1, fc1, 0.f};
0280     Vec3F ssc1{fs1, fs1, fc1, 0.f};
0281     const Vec3F cslsimd = scc1 * ssc1 * float(l1);
0282
0283     for (int i = 0; i != cacheSize; ++i) {
0284       StereoInfo const si = cache[i];
0285
0286       // match
0287       Vec2D minv10{si.m11, -si.m10};
0288       double mult = 1. / (m00 * si.m11 - m01 * si.m10);
0289       Vec2D resultmatmul = mult * (minv10 * c0vec + minv00 * si.c1vec);
0290
0291       Local2DPoint position(resultmatmul[0], resultmatmul[1]);
0292
0293       // LocalError tempError (100,0,100);
0294       if (!((gluedDet->surface()).bounds().inside(position, 10.f * scale_)))
0295         continue;
0296
0297       double c2 = -si.m10;
0298       double s2 = -si.m11;
0299       double l2 = 1. / (c2 * c2 + s2 * s2);
0300
0301       double diff = (c1 * s2 - c2 * s1);
0302       double invdet2 = 1. / (diff * diff * l1 * l2);
0303
0304       float fc2(c2), fs2(s2), fid2(invdet2);
0305       Vec3F invdet2simd{fid2, -fid2, fid2, 0.f};
0306       Vec3F ccssimd{fs2, fc2, fc2, 0.f};
0307       Vec3F csssimd{fs2, fs2, fc2, 0.f};
0308       Vec3F result = invdet2simd * (si.sigmap22 * cslsimd + sigmap12 * ccssimd * csssimd * float(l2));
0309
0310       LocalError error(result[0], result[1], result[2]);
0311
0312       if ((gluedDet->surface()).bounds().inside(position, error, scale_)) {  //if it is inside the gluedet bonds
0313
0314         //Change NSigmaInside in the configuration file to accept more hits
0315         //...and add it to the Rechit collection
0316
0317         collectorHelper.collector()(
0318             SiStripMatchedRecHit2D(LocalPoint(position), error, *gluedDet, &monoRH, si.secondHit));
0319       }
0320
0321     }  // loop on cache info
0322
0323     collectorHelper.closure(monoRHiter);
0324   }  // loop on mono hit
0325 }
0326
0327 #endif  //DOUBLE_MATCH
0328
0329 #endif  //  RECOLOCALTRACKER_SISTRIPCLUSTERIZER_SISTRIPRECHITMATCH_H