TrackingRecHitProducer/plugins/FastTrackerRecHitMatcher.cc

0001 // system include files
0002 #include <cfloat>
0003 #include <memory>
0004
0005 // framework stuff
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 // fast tracker rechits
0013 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0014 #include "DataFormats/TrackerRecHit2D/interface/FastMatchedTrackerRecHit.h"
0015 #include "DataFormats/TrackerRecHit2D/interface/FastProjectedTrackerRecHit.h"
0016 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHitCollection.h"
0017
0018 // geometry stuff
0019 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0020 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0021 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0022 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0023 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0024 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0025
0026 // sim stuff
0027 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0028 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0029
0030 class FastTrackerRecHitMatcher : public edm::stream::EDProducer<> {
0031 public:
0032   explicit FastTrackerRecHitMatcher(const edm::ParameterSet&);
0033   ~FastTrackerRecHitMatcher() override { ; }
0034   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0035
0036 private:
0037   void produce(edm::Event&, const edm::EventSetup&) override;
0038
0039   // ---------- typedefs -----------------------------
0040   typedef std::pair<LocalPoint, LocalPoint> StripPosition;
0041
0042   // ----------internal functions ---------------------------
0043
0044   // create projected hit
0045   std::unique_ptr<FastTrackerRecHit> projectOnly(const FastSingleTrackerRecHit* originalRH,
0046                                                  const GeomDet* monoDet,
0047                                                  const GluedGeomDet* gluedDet,
0048                                                  LocalVector& ldir) const;
0049
0050   // create matched hit
0051   std::unique_ptr<FastTrackerRecHit> match(const FastSingleTrackerRecHit* monoRH,
0052                                            const FastSingleTrackerRecHit* stereoRH,
0053                                            const GluedGeomDet* gluedDet,
0054                                            LocalVector& trackdirection,
0055                                            bool stereLayerFirst) const;
0056
0057   StripPosition project(const GeomDetUnit* det,
0058                         const GluedGeomDet* glueddet,
0059                         const StripPosition& strip,
0060                         const LocalVector& trackdirection) const;
0061
0062   inline const FastSingleTrackerRecHit* _cast2Single(const FastTrackerRecHit* recHit) const {
0063     if (!recHit->isSingle()) {
0064       throw cms::Exception("FastTrackerRecHitMatcher")
0065           << "all rechits in simHit2RecHitMap must be instances of FastSingleTrackerRecHit. recHit's rtti: "
0066           << recHit->rtti() << std::endl;
0067     }
0068     return dynamic_cast<const FastSingleTrackerRecHit*>(recHit);
0069   }
0070
0071   // ----------member data ---------------------------
0072   edm::EDGetTokenT<edm::PSimHitContainer> simHitsToken;
0073   edm::EDGetTokenT<FastTrackerRecHitRefCollection> simHit2RecHitMapToken;
0074   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryESToken;
0075 };
0076
0077 FastTrackerRecHitMatcher::FastTrackerRecHitMatcher(const edm::ParameterSet& iConfig)
0078     : trackerGeometryESToken(esConsumes()) {
0079   simHitsToken = consumes<edm::PSimHitContainer>(iConfig.getParameter<edm::InputTag>("simHits"));
0080   simHit2RecHitMapToken =
0081       consumes<FastTrackerRecHitRefCollection>(iConfig.getParameter<edm::InputTag>("simHit2RecHitMap"));
0082
0083   produces<FastTrackerRecHitCollection>();
0084   produces<FastTrackerRecHitRefCollection>("simHit2RecHitMap");
0085 }
0086
0087 void FastTrackerRecHitMatcher::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0088   // services
0089   auto const& geometry = iSetup.getData(trackerGeometryESToken);
0090
0091   // input
0092   edm::Handle<edm::PSimHitContainer> simHits;
0093   iEvent.getByToken(simHitsToken, simHits);
0094
0095   edm::Handle<FastTrackerRecHitRefCollection> simHit2RecHitMap;
0096   iEvent.getByToken(simHit2RecHitMapToken, simHit2RecHitMap);
0097
0098   // output
0099   auto output_recHits = std::make_unique<FastTrackerRecHitCollection>();
0100   auto output_simHit2RecHitMap =
0101       std::make_unique<FastTrackerRecHitRefCollection>(simHit2RecHitMap->size(), FastTrackerRecHitRef());
0102   edm::RefProd<FastTrackerRecHitCollection> output_recHits_refProd =
0103       iEvent.getRefBeforePut<FastTrackerRecHitCollection>();
0104
0105   bool skipNext = false;
0106   for (unsigned simHitCounter = 0; simHitCounter < simHits->size(); ++simHitCounter) {
0107     // skip hit in case it was matched to previous one
0108     if (skipNext) {
0109       skipNext = false;
0110       continue;
0111     }
0112     skipNext = false;
0113
0114     // get simHit and associated recHit
0115     const PSimHit& simHit = (*simHits)[simHitCounter];
0116     const FastTrackerRecHitRef& recHitRef = (*simHit2RecHitMap)[simHitCounter];
0117
0118     // skip simHits w/o associated recHit
0119     if (recHitRef.isNull())
0120       continue;
0121
0122     // cast
0123     const FastSingleTrackerRecHit* recHit = _cast2Single(recHitRef.get());
0124
0125     // get subdetector id
0126     DetId detid = recHit->geographicalId();
0127     unsigned int subdet = detid.subdetId();
0128
0129     // treat pixel hits
0130     if (subdet <= 2) {
0131       (*output_simHit2RecHitMap)[simHitCounter] = recHitRef;
0132     }
0133
0134     // treat strip hits
0135     else {
0136       StripSubdetector stripSubDetId(detid);
0137
0138       // treat regular regular strip hits
0139       if (!stripSubDetId.glued()) {
0140         (*output_simHit2RecHitMap)[simHitCounter] = recHitRef;
0141       }
0142
0143       // treat strip hits on glued layers
0144       else {
0145         // Obtain direction of simtrack at simhit in local coordinates of glued module
0146         //   - direction of simtrack at simhit, in coordindates of the single module
0147         LocalVector localSimTrackDir = simHit.localDirection();
0148         //   - transform to global coordinates
0149         GlobalVector globalSimTrackDir = recHit->det()->surface().toGlobal(localSimTrackDir);
0150         //   - transform to local coordinates of glued module
0151         const GluedGeomDet* gluedDet = (const GluedGeomDet*)geometry.idToDet(DetId(stripSubDetId.glued()));
0152         LocalVector gluedLocalSimTrackDir = gluedDet->surface().toLocal(globalSimTrackDir);
0153
0154         // check whether next hit is partner
0155         const FastSingleTrackerRecHit* partnerRecHit = nullptr;
0156         //      - there must be a next hit
0157         if (simHitCounter + 1 < simHits->size()) {
0158           const FastTrackerRecHitRef& nextRecHitRef = (*simHit2RecHitMap)[simHitCounter + 1];
0159           const PSimHit& nextSimHit = (*simHits)[simHitCounter + 1];
0160           //  - partner hit must not be null
0161           //  - simHit and partner simHit must belong to same simTrack
0162           //  - partner hit must be on the module glued to the module of the hit
0163           if ((!nextRecHitRef.isNull()) && simHit.trackId() == nextSimHit.trackId() &&
0164               StripSubdetector(nextRecHitRef->geographicalId()).partnerDetId() == detid.rawId()) {
0165             partnerRecHit = _cast2Single(nextRecHitRef.get());
0166             skipNext = true;
0167           }
0168         }
0169
0170         std::unique_ptr<FastTrackerRecHit> newRecHit(nullptr);
0171
0172         // if partner found: create a matched hit
0173         if (partnerRecHit) {
0174           newRecHit = match(stripSubDetId.stereo() ? partnerRecHit : recHit,
0175                             stripSubDetId.stereo() ? recHit : partnerRecHit,
0176                             gluedDet,
0177                             gluedLocalSimTrackDir,
0178                             stripSubDetId.stereo());
0179         }
0180         // else: create projected hit
0181         else {
0182           newRecHit = projectOnly(recHit, geometry.idToDet(detid), gluedDet, gluedLocalSimTrackDir);
0183         }
0184         output_recHits->push_back(std::move(newRecHit));
0185         (*output_simHit2RecHitMap)[simHitCounter] =
0186             FastTrackerRecHitRef(output_recHits_refProd, output_recHits->size() - 1);
0187       }
0188     }
0189   }
0190
0191   iEvent.put(std::move(output_recHits));
0192   iEvent.put(std::move(output_simHit2RecHitMap), "simHit2RecHitMap");
0193 }
0194
0195 std::unique_ptr<FastTrackerRecHit> FastTrackerRecHitMatcher::match(const FastSingleTrackerRecHit* monoRH,
0196                                                                    const FastSingleTrackerRecHit* stereoRH,
0197                                                                    const GluedGeomDet* gluedDet,
0198                                                                    LocalVector& trackdirection,
0199                                                                    bool stereoHitFirst) const {
0200   // stripdet = mono
0201   // partnerstripdet = stereo
0202   const GeomDetUnit* stripdet = gluedDet->monoDet();
0203   const GeomDetUnit* partnerstripdet = gluedDet->stereoDet();
0204   const StripTopology& topol = (const StripTopology&)stripdet->topology();
0205
0206   LocalPoint position;
0207
0208   // position of the initial and final point of the strip (RPHI cluster) in local strip coordinates
0209   MeasurementPoint RPHIpoint = topol.measurementPosition(monoRH->localPosition());
0210   MeasurementPoint RPHIpointini = MeasurementPoint(RPHIpoint.x(), -0.5);
0211   MeasurementPoint RPHIpointend = MeasurementPoint(RPHIpoint.x(), 0.5);
0212
0213   // position of the initial and final point of the strip in local coordinates (mono det)
0214   StripPosition stripmono = StripPosition(topol.localPosition(RPHIpointini), topol.localPosition(RPHIpointend));
0215
0216   // in case of no track hypothesis assume a track from the origin through the center of the strip
0217   if (trackdirection.mag2() < FLT_MIN) {
0218     LocalPoint lcenterofstrip = monoRH->localPosition();
0219     GlobalPoint gcenterofstrip = (stripdet->surface()).toGlobal(lcenterofstrip);
0220     GlobalVector gtrackdirection = gcenterofstrip - GlobalPoint(0, 0, 0);
0221     trackdirection = (gluedDet->surface()).toLocal(gtrackdirection);
0222   }
0223
0224   //project mono hit on glued det
0225   StripPosition projectedstripmono = project(stripdet, gluedDet, stripmono, trackdirection);
0226   const StripTopology& partnertopol = (const StripTopology&)partnerstripdet->topology();
0227
0228   //error calculation (the part that depends on mono RH only)
0229   LocalVector RPHIpositiononGluedendvector = projectedstripmono.second - projectedstripmono.first;
0230   double c1 = sin(RPHIpositiononGluedendvector.phi());
0231   double s1 = -cos(RPHIpositiononGluedendvector.phi());
0232   MeasurementError errormonoRH = topol.measurementError(monoRH->localPosition(), monoRH->localPositionError());
0233   double pitch = topol.localPitch(monoRH->localPosition());
0234   double sigmap12 = errormonoRH.uu() * pitch * pitch;
0235
0236   // position of the initial and final point of the strip (STEREO cluster)
0237   MeasurementPoint STEREOpoint = partnertopol.measurementPosition(stereoRH->localPosition());
0238
0239   MeasurementPoint STEREOpointini = MeasurementPoint(STEREOpoint.x(), -0.5);
0240   MeasurementPoint STEREOpointend = MeasurementPoint(STEREOpoint.x(), 0.5);
0241
0242   // position of the initial and final point of the strip in local coordinates (stereo det)
0243   StripPosition stripstereo(partnertopol.localPosition(STEREOpointini), partnertopol.localPosition(STEREOpointend));
0244
0245   //project stereo hit on glued det
0246   StripPosition projectedstripstereo = project(partnerstripdet, gluedDet, stripstereo, trackdirection);
0247
0248   //perform the matching
0249   //(x2-x1)(y-y1)=(y2-y1)(x-x1)
0250   AlgebraicMatrix22 m;
0251   AlgebraicVector2 c, solution;
0252   m(0, 0) = -(projectedstripmono.second.y() - projectedstripmono.first.y());
0253   m(0, 1) = (projectedstripmono.second.x() - projectedstripmono.first.x());
0254   m(1, 0) = -(projectedstripstereo.second.y() - projectedstripstereo.first.y());
0255   m(1, 1) = (projectedstripstereo.second.x() - projectedstripstereo.first.x());
0256   c(0) = m(0, 1) * projectedstripmono.first.y() + m(0, 0) * projectedstripmono.first.x();
0257   c(1) = m(1, 1) * projectedstripstereo.first.y() + m(1, 0) * projectedstripstereo.first.x();
0258   m.Invert();
0259   solution = m * c;
0260   position = LocalPoint(solution(0), solution(1));
0261
0262   //
0263   // temporary fix by tommaso
0264   //
0265
0266   LocalError tempError(100, 0, 100);
0267
0268   // calculate the error
0269   LocalVector stereopositiononGluedendvector = projectedstripstereo.second - projectedstripstereo.first;
0270   double c2 = sin(stereopositiononGluedendvector.phi());
0271   double s2 = -cos(stereopositiononGluedendvector.phi());
0272   MeasurementError errorstereoRH =
0273       partnertopol.measurementError(stereoRH->localPosition(), stereoRH->localPositionError());
0274   pitch = partnertopol.localPitch(stereoRH->localPosition());
0275   double sigmap22 = errorstereoRH.uu() * pitch * pitch;
0276   double diff = (c1 * s2 - c2 * s1);
0277   double invdet2 = 1 / (diff * diff);
0278   float xx = invdet2 * (sigmap12 * s2 * s2 + sigmap22 * s1 * s1);
0279   float xy = -invdet2 * (sigmap12 * c2 * s2 + sigmap22 * c1 * s1);
0280   float yy = invdet2 * (sigmap12 * c2 * c2 + sigmap22 * c1 * c1);
0281   LocalError error = LocalError(xx, xy, yy);
0282
0283   //Added by DAO to make sure y positions are zero.
0284   DetId det(monoRH->geographicalId());
0285   if (det.subdetId() > 2) {
0286     return std::make_unique<FastMatchedTrackerRecHit>(position, error, *gluedDet, *monoRH, *stereoRH, stereoHitFirst);
0287
0288   }
0289
0290   else {
0291     throw cms::Exception("FastTrackerRecHitMatcher") << "Matched Pixel!?";
0292   }
0293 }
0294
0295 FastTrackerRecHitMatcher::StripPosition FastTrackerRecHitMatcher::project(const GeomDetUnit* det,
0296                                                                           const GluedGeomDet* glueddet,
0297                                                                           const StripPosition& strip,
0298                                                                           const LocalVector& trackdirection) const {
0299   GlobalPoint globalpointini = (det->surface()).toGlobal(strip.first);
0300   GlobalPoint globalpointend = (det->surface()).toGlobal(strip.second);
0301
0302   // position of the initial and final point of the strip in glued local coordinates
0303   LocalPoint positiononGluedini = (glueddet->surface()).toLocal(globalpointini);
0304   LocalPoint positiononGluedend = (glueddet->surface()).toLocal(globalpointend);
0305
0306   //correct the position with the track direction
0307
0308   float scale = -positiononGluedini.z() / trackdirection.z();
0309
0310   LocalPoint projpositiononGluedini = positiononGluedini + scale * trackdirection;
0311   LocalPoint projpositiononGluedend = positiononGluedend + scale * trackdirection;
0312
0313   return StripPosition(projpositiononGluedini, projpositiononGluedend);
0314 }
0315
0316 std::unique_ptr<FastTrackerRecHit> FastTrackerRecHitMatcher::projectOnly(const FastSingleTrackerRecHit* originalRH,
0317                                                                          const GeomDet* monoDet,
0318                                                                          const GluedGeomDet* gluedDet,
0319                                                                          LocalVector& ldir) const {
0320   LocalPoint position(originalRH->localPosition().x(), 0., 0.);
0321   const BoundPlane& gluedPlane = gluedDet->surface();
0322   const BoundPlane& hitPlane = monoDet->surface();
0323
0324   double delta = gluedPlane.localZ(hitPlane.position());
0325
0326   LocalPoint lhitPos = gluedPlane.toLocal(monoDet->surface().toGlobal(position));
0327   LocalPoint projectedHitPos = lhitPos - ldir * delta / ldir.z();
0328
0329   LocalVector hitXAxis = gluedPlane.toLocal(hitPlane.toGlobal(LocalVector(1, 0, 0)));
0330   LocalError hitErr = originalRH->localPositionError();
0331
0332   if (gluedPlane.normalVector().dot(hitPlane.normalVector()) < 0) {
0333     // the two planes are inverted, and the correlation element must change sign
0334     hitErr = LocalError(hitErr.xx(), -hitErr.xy(), hitErr.yy());
0335   }
0336   LocalError rotatedError = hitErr.rotate(hitXAxis.x(), hitXAxis.y());
0337
0338   const GeomDetUnit* gluedMonoDet = gluedDet->monoDet();
0339   const GeomDetUnit* gluedStereoDet = gluedDet->stereoDet();
0340   int isMono = 0;
0341   int isStereo = 0;
0342
0343   if (monoDet->geographicalId() == gluedMonoDet->geographicalId())
0344     isMono = 1;
0345   if (monoDet->geographicalId() == gluedStereoDet->geographicalId())
0346     isStereo = 1;
0347   //Added by DAO to make sure y positions are zero and correct Mono or stereo Det is filled.
0348
0349   if ((isMono && isStereo) || (!isMono && !isStereo))
0350     throw cms::Exception("FastTrackerRecHitMatcher") << "Something wrong with DetIds.";
0351   return std::make_unique<FastProjectedTrackerRecHit>(projectedHitPos, rotatedError, *gluedDet, *originalRH);
0352 }
0353
0354 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0355 void FastTrackerRecHitMatcher::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0356   //The following says we do not know what parameters are allowed so do no validation
0357   // Please change this to state exactly what you do use, even if it is no parameters
0358   edm::ParameterSetDescription desc;
0359   desc.setUnknown();
0360   descriptions.addDefault(desc);
0361 }
0362
0363 //define this as a plug-in
0364 DEFINE_FWK_MODULE(FastTrackerRecHitMatcher);