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