SiStripRecHitConverter/src/SiStripRecHitMatcher.cc

0001 // File: SiStripRecHitMatcher.cc
0002 // Description:  Matches into rechits
0003 // Author:  C.Genta
0004 #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripRecHitMatcher.h"
0005 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0006 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0008 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0009
0010 #include "TrackingTools/TransientTrackingRecHit/interface/HelpertRecHit2DLocalPos.h"
0011 #include <functional>
0012
0013 #include <DataFormats/GeometryCommonDetAlgo/interface/AlignmentPositionError.h>
0014
0015 SiStripRecHitMatcher::SiStripRecHitMatcher(const edm::ParameterSet& conf)
0016     : scale_(conf.getParameter<double>("NSigmaInside")),
0017       preFilter_(conf.existsAs<bool>("PreFilter") ? conf.getParameter<bool>("PreFilter") : false) {}
0018
0019 SiStripRecHitMatcher::SiStripRecHitMatcher(const double theScale) : scale_(theScale) {}
0020
0021 namespace {
0022   // FIXME for c++0X
0023   inline void pb1(std::vector<SiStripMatchedRecHit2D*>& v, SiStripMatchedRecHit2D* h) { v.push_back(h); }
0024   inline void pb2(SiStripRecHitMatcher::CollectorMatched& v, const SiStripMatchedRecHit2D& h) { v.push_back(h); }
0025
0026 }  // namespace
0027
0028 // needed by the obsolete version still in use on some architectures
0029 void SiStripRecHitMatcher::match(const SiStripRecHit2D* monoRH,
0030                                  SimpleHitIterator begin,
0031                                  SimpleHitIterator end,
0032                                  std::vector<std::unique_ptr<SiStripMatchedRecHit2D>>& collector,
0033                                  const GluedGeomDet* gluedDet,
0034                                  LocalVector trackdirection) const {
0035   std::vector<SiStripMatchedRecHit2D*> result;
0036   result.reserve(end - begin);
0037   match(monoRH, begin, end, result, gluedDet, trackdirection);
0038   for (std::vector<SiStripMatchedRecHit2D*>::iterator p = result.begin(); p != result.end(); p++)
0039     collector.emplace_back(*p);
0040 }
0041
0042 void SiStripRecHitMatcher::match(const SiStripRecHit2D* monoRH,
0043                                  SimpleHitIterator begin,
0044                                  SimpleHitIterator end,
0045                                  std::vector<SiStripMatchedRecHit2D*>& collector,
0046                                  const GluedGeomDet* gluedDet,
0047                                  LocalVector trackdirection) const {
0048   Collector result(
0049       std::bind(&pb1, std::ref(collector), std::bind(&SiStripMatchedRecHit2D::clone, std::placeholders::_1)));
0050   match(monoRH, begin, end, result, gluedDet, trackdirection);
0051 }
0052
0053 void SiStripRecHitMatcher::match(const SiStripRecHit2D* monoRH,
0054                                  SimpleHitIterator begin,
0055                                  SimpleHitIterator end,
0056                                  CollectorMatched& collector,
0057                                  const GluedGeomDet* gluedDet,
0058                                  LocalVector trackdirection) const {
0059   Collector result(std::bind(pb2, std::ref(collector), std::placeholders::_1));
0060   match(monoRH, begin, end, result, gluedDet, trackdirection);
0061 }
0062
0063 // this is the one used by the RecHitConverter
0064 void SiStripRecHitMatcher::match(const SiStripRecHit2D* monoRH,
0065                                  RecHitIterator begin,
0066                                  RecHitIterator end,
0067                                  CollectorMatched& collector,
0068                                  const GluedGeomDet* gluedDet,
0069                                  LocalVector trackdirection) const {
0070   // is this really needed now????
0071   SimpleHitCollection stereoHits;
0072   stereoHits.reserve(end - begin);
0073   for (RecHitIterator i = begin; i != end; ++i)
0074     stereoHits.push_back(&(*i));  // convert to simple pointer
0075
0076   return match(monoRH, stereoHits.begin(), stereoHits.end(), collector, gluedDet, trackdirection);
0077 }
0078
0079 // the "real implementation"
0080 void SiStripRecHitMatcher::match(const SiStripRecHit2D* monoRH,
0081                                  SimpleHitIterator begin,
0082                                  SimpleHitIterator end,
0083                                  Collector& collector,
0084                                  const GluedGeomDet* gluedDet,
0085                                  LocalVector trackdirection) const {
0086   // stripdet = mono
0087   // partnerstripdet = stereo
0088   const GeomDetUnit* stripdet = gluedDet->monoDet();
0089   const GeomDetUnit* partnerstripdet = gluedDet->stereoDet();
0090   const StripTopology& topol = (const StripTopology&)stripdet->topology();
0091
0092   // position of the initial and final point of the strip (RPHI cluster) in local strip coordinates
0093   double RPHIpointX = topol.measurementPosition(monoRH->localPositionFast()).x();
0094   MeasurementPoint RPHIpointini(RPHIpointX, -0.5);
0095   MeasurementPoint RPHIpointend(RPHIpointX, 0.5);
0096
0097   // position of the initial and final point of the strip in local coordinates (mono det)
0098   StripPosition stripmono = StripPosition(topol.localPosition(RPHIpointini), topol.localPosition(RPHIpointend));
0099
0100   // in case of no track hypothesis assume a track from the origin through the center of the strip
0101   if (trackdirection.mag2() < FLT_MIN) {
0102     const LocalPoint& lcenterofstrip = monoRH->localPositionFast();
0103     GlobalPoint gcenterofstrip = (stripdet->surface()).toGlobal(lcenterofstrip);
0104     GlobalVector gtrackdirection = gcenterofstrip - GlobalPoint(0, 0, 0);
0105     trackdirection = (gluedDet->surface()).toLocal(gtrackdirection);
0106   }
0107
0108   //project mono hit on glued det
0109   StripPosition projectedstripmono = project(stripdet, gluedDet, stripmono, trackdirection);
0110   const StripTopology& partnertopol = (const StripTopology&)partnerstripdet->topology();
0111
0112   double m00 = -(projectedstripmono.second.y() - projectedstripmono.first.y());
0113   double m01 = (projectedstripmono.second.x() - projectedstripmono.first.x());
0114   double c0 = m01 * projectedstripmono.first.y() + m00 * projectedstripmono.first.x();
0115
0116   //error calculation (the part that depends on mono RH only)
0117   //  LocalVector  RPHIpositiononGluedendvector=projectedstripmono.second-projectedstripmono.first;
0118   /*
0119   double l1 = 1./RPHIpositiononGluedendvector.perp2();
0120   double c1 = RPHIpositiononGluedendvector.y();
0121   double s1 =-RPHIpositiononGluedendvector.x();
0122   */
0123   double c1 = -m00;
0124   double s1 = -m01;
0125   double l1 = 1. / (c1 * c1 + s1 * s1);
0126
0127   double sigmap12 = sigmaPitch(monoRH->localPosition(), monoRH->localPositionError(), topol);
0128   // auto sigmap12 = monoRH->sigmaPitch();
0129   // assert(sigmap12>=0);
0130
0131   SimpleHitIterator seconditer;
0132
0133   for (seconditer = begin; seconditer != end; ++seconditer) {  //iterate on stereo rechits
0134
0135     // position of the initial and final point of the strip (STEREO cluster)
0136     double STEREOpointX = partnertopol.measurementPosition((*seconditer)->localPositionFast()).x();
0137     MeasurementPoint STEREOpointini(STEREOpointX, -0.5);
0138     MeasurementPoint STEREOpointend(STEREOpointX, 0.5);
0139
0140     // position of the initial and final point of the strip in local coordinates (stereo det)
0141     StripPosition stripstereo(partnertopol.localPosition(STEREOpointini), partnertopol.localPosition(STEREOpointend));
0142
0143     //project stereo hit on glued det
0144     StripPosition projectedstripstereo = project(partnerstripdet, gluedDet, stripstereo, trackdirection);
0145
0146     double m10 = -(projectedstripstereo.second.y() - projectedstripstereo.first.y());
0147     double m11 = (projectedstripstereo.second.x() - projectedstripstereo.first.x());
0148
0149     //perform the matching
0150     //(x2-x1)(y-y1)=(y2-y1)(x-x1)
0151     AlgebraicMatrix22 m;
0152     AlgebraicVector2 c;  // FIXME understand why moving this initializer out of the loop changes the output!
0153     m(0, 0) = m00;
0154     m(0, 1) = m01;
0155     m(1, 0) = m10;
0156     m(1, 1) = m11;
0157     c(0) = c0;
0158     c(1) = m11 * projectedstripstereo.first.y() + m10 * projectedstripstereo.first.x();
0159     m.Invert();
0160     AlgebraicVector2 solution = m * c;
0161     LocalPoint position(solution(0), solution(1));
0162
0163     /*
0164     {
0165       double m00 = -(projectedstripmono.second.y()-projectedstripmono.first.y());
0166       double m01 =  (projectedstripmono.second.x()-projectedstripmono.first.x());
0167       double m10 = -(projectedstripstereo.second.y()-projectedstripstereo.first.y());
0168       double m11 =  (projectedstripstereo.second.x()-projectedstripstereo.first.x());
0169       double c0  =  m01*projectedstripmono.first.y()   + m00*projectedstripmono.first.x();
0170       double c1  =  m11*projectedstripstereo.first.y() + m10*projectedstripstereo.first.x();
0171
0172       double invDet = 1./(m00*m11-m10*m01);
0173     }
0174     */
0175
0176     //
0177     // temporary fix by tommaso
0178     //
0179
0180     LocalError tempError(100, 0, 100);
0181     if (!((gluedDet->surface()).bounds().inside(position, tempError, scale_)))
0182       continue;
0183
0184     // then calculate the error
0185     /*
0186     LocalVector  stereopositiononGluedendvector=projectedstripstereo.second-projectedstripstereo.first;
0187     double l2 = 1./stereopositiononGluedendvector.perp2();
0188     double c2 = stereopositiononGluedendvector.y();
0189     double s2 =-stereopositiononGluedendvector.x();
0190     */
0191
0192     double c2 = -m10;
0193     double s2 = -m11;
0194     double l2 = 1. / (c2 * c2 + s2 * s2);
0195
0196     double sigmap22 = sigmaPitch((*seconditer)->localPosition(), (*seconditer)->localPositionError(), partnertopol);
0197     // auto sigmap22 = (*seconditer)->sigmaPitch();
0198     // assert(sigmap22>=0);
0199
0200     double diff = (c1 * s2 - c2 * s1);
0201     double invdet2 = 1 / (diff * diff * l1 * l2);
0202     float xx = invdet2 * (sigmap12 * s2 * s2 * l2 + sigmap22 * s1 * s1 * l1);
0203     float xy = -invdet2 * (sigmap12 * c2 * s2 * l2 + sigmap22 * c1 * s1 * l1);
0204     float yy = invdet2 * (sigmap12 * c2 * c2 * l2 + sigmap22 * c1 * c1 * l1);
0205     LocalError error(xx, xy, yy);
0206
0207     if ((gluedDet->surface()).bounds().inside(position, error, scale_)) {  //if it is inside the gluedet bonds
0208       //Change NSigmaInside in the configuration file to accept more hits
0209       //...and add it to the Rechit collection
0210
0211       const SiStripRecHit2D* secondHit = *seconditer;
0212       collector(SiStripMatchedRecHit2D(position, error, *gluedDet, monoRH, secondHit));
0213     }
0214   }
0215 }
0216
0217 SiStripRecHitMatcher::StripPosition SiStripRecHitMatcher::project(const GeomDetUnit* det,
0218                                                                   const GluedGeomDet* glueddet,
0219                                                                   StripPosition strip,
0220                                                                   LocalVector trackdirection) const {
0221   GlobalPoint globalpointini = (det->surface()).toGlobal(strip.first);
0222   GlobalPoint globalpointend = (det->surface()).toGlobal(strip.second);
0223
0224   // position of the initial and final point of the strip in glued local coordinates
0225   LocalPoint positiononGluedini = (glueddet->surface()).toLocal(globalpointini);
0226   LocalPoint positiononGluedend = (glueddet->surface()).toLocal(globalpointend);
0227
0228   //correct the position with the track direction
0229
0230   float scale = -positiononGluedini.z() / trackdirection.z();
0231
0232   LocalPoint projpositiononGluedini = positiononGluedini + scale * trackdirection;
0233   LocalPoint projpositiononGluedend = positiononGluedend + scale * trackdirection;
0234
0235   return StripPosition(projpositiononGluedini, projpositiononGluedend);
0236 }
0237
0238 //match a single hit
0239 std::unique_ptr<SiStripMatchedRecHit2D> SiStripRecHitMatcher::match(const SiStripRecHit2D* monoRH,
0240                                                                     const SiStripRecHit2D* stereoRH,
0241                                                                     const GluedGeomDet* gluedDet,
0242                                                                     LocalVector trackdirection,
0243                                                                     bool force) const {
0244   // stripdet = mono
0245   // partnerstripdet = stereo
0246   const GeomDetUnit* stripdet = gluedDet->monoDet();
0247   const GeomDetUnit* partnerstripdet = gluedDet->stereoDet();
0248   const StripTopology& topol = (const StripTopology&)stripdet->topology();
0249
0250   // position of the initial and final point of the strip (RPHI cluster) in local strip coordinates
0251   auto RPHIpointX = topol.measurementPosition(monoRH->localPositionFast()).x();
0252   MeasurementPoint RPHIpointini(RPHIpointX, -0.5f);
0253   MeasurementPoint RPHIpointend(RPHIpointX, 0.5f);
0254
0255   // position of the initial and final point of the strip in local coordinates (mono det)
0256   StripPosition stripmono = StripPosition(topol.localPosition(RPHIpointini), topol.localPosition(RPHIpointend));
0257
0258   // in case of no track hypothesis assume a track from the origin through the center of the strip
0259   if (trackdirection.mag2() < FLT_MIN) {
0260     const LocalPoint& lcenterofstrip = monoRH->localPositionFast();
0261     GlobalPoint gcenterofstrip = (stripdet->surface()).toGlobal(lcenterofstrip);
0262     GlobalVector gtrackdirection = gcenterofstrip - GlobalPoint(0, 0, 0);
0263     trackdirection = (gluedDet->surface()).toLocal(gtrackdirection);
0264   }
0265
0266   //project mono hit on glued det
0267   StripPosition projectedstripmono = project(stripdet, gluedDet, stripmono, trackdirection);
0268   const StripTopology& partnertopol = (const StripTopology&)partnerstripdet->topology();
0269
0270   double m00 = -(projectedstripmono.second.y() - projectedstripmono.first.y());
0271   double m01 = (projectedstripmono.second.x() - projectedstripmono.first.x());
0272   double c0 = m01 * projectedstripmono.first.y() + m00 * projectedstripmono.first.x();
0273
0274   //error calculation (the part that depends on mono RH only)
0275   //  LocalVector  RPHIpositiononGluedendvector=projectedstripmono.second-projectedstripmono.first;
0276   /*
0277   double l1 = 1./RPHIpositiononGluedendvector.perp2();
0278   double c1 = RPHIpositiononGluedendvector.y();
0279   double s1 =-RPHIpositiononGluedendvector.x();
0280   */
0281   double c1 = -m00;
0282   double s1 = -m01;
0283   double l1 = 1. / (c1 * c1 + s1 * s1);
0284
0285   double sigmap12 = sigmaPitch(monoRH->localPosition(), monoRH->localPositionError(), topol);
0286   // auto sigmap12 = monoRH->sigmaPitch();
0287   // assert(sigmap12>=0);
0288
0289   // position of the initial and final point of the strip (STEREO cluster)
0290   auto STEREOpointX = partnertopol.measurementPosition(stereoRH->localPositionFast()).x();
0291   MeasurementPoint STEREOpointini(STEREOpointX, -0.5f);
0292   MeasurementPoint STEREOpointend(STEREOpointX, 0.5f);
0293
0294   // position of the initial and final point of the strip in local coordinates (stereo det)
0295   StripPosition stripstereo(partnertopol.localPosition(STEREOpointini), partnertopol.localPosition(STEREOpointend));
0296
0297   //project stereo hit on glued det
0298   StripPosition projectedstripstereo = project(partnerstripdet, gluedDet, stripstereo, trackdirection);
0299
0300   double m10 = -(projectedstripstereo.second.y() - projectedstripstereo.first.y());
0301   double m11 = (projectedstripstereo.second.x() - projectedstripstereo.first.x());
0302
0303   //perform the matching
0304   //(x2-x1)(y-y1)=(y2-y1)(x-x1)
0305   AlgebraicMatrix22 m(ROOT::Math::SMatrixNoInit{});
0306   AlgebraicVector2 c(c0, m11 * projectedstripstereo.first.y() + m10 * projectedstripstereo.first.x());
0307   m(0, 0) = m00;
0308   m(0, 1) = m01;
0309   m(1, 0) = m10;
0310   m(1, 1) = m11;
0311   m.Invert();
0312   AlgebraicVector2 solution = m * c;
0313   Local2DPoint position(solution(0), solution(1));
0314
0315   if ((!force) && (!((gluedDet->surface()).bounds().inside(position, 10.f * scale_))))
0316     return std::unique_ptr<SiStripMatchedRecHit2D>(nullptr);
0317
0318   double c2 = -m10;
0319   double s2 = -m11;
0320   double l2 = 1. / (c2 * c2 + s2 * s2);
0321
0322   double sigmap22 = sigmaPitch(stereoRH->localPosition(), stereoRH->localPositionError(), partnertopol);
0323   // auto sigmap22 = stereoRH->sigmaPitch();
0324   // assert (sigmap22>0);
0325
0326   double diff = (c1 * s2 - c2 * s1);
0327   double invdet2 = 1 / (diff * diff * l1 * l2);
0328   float xx = invdet2 * (sigmap12 * s2 * s2 * l2 + sigmap22 * s1 * s1 * l1);
0329   float xy = -invdet2 * (sigmap12 * c2 * s2 * l2 + sigmap22 * c1 * s1 * l1);
0330   float yy = invdet2 * (sigmap12 * c2 * c2 * l2 + sigmap22 * c1 * c1 * l1);
0331   LocalError error(xx, xy, yy);
0332
0333   //if it is inside the gluedet bonds
0334   //Change NSigmaInside in the configuration file to accept more hits
0335   if (force || (gluedDet->surface()).bounds().inside(position, error, scale_))
0336     return std::make_unique<SiStripMatchedRecHit2D>(LocalPoint(position), error, *gluedDet, monoRH, stereoRH);
0337   return std::unique_ptr<SiStripMatchedRecHit2D>(nullptr);
0338 }