TrackProducer/src/GsfTrackProducerBase.cc

0001 #include "RecoTracker/TrackProducer/interface/GsfTrackProducerBase.h"
0002 // system include files
0003 #include <memory>
0004 // user include files
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "DataFormats/GsfTrackReco/interface/GsfTrackExtra.h"
0009 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0010 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0011 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0012 #include "TrackingTools/GsfTools/interface/GsfPropagatorAdapter.h"
0013 #include "TrackingTools/GsfTools/interface/MultiGaussianStateTransform.h"
0014 #include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
0015 #include "TrackingTools/GsfTools/interface/GaussianSumUtilities1D.h"
0016 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0017 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0018 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0019 #include "TrackingTools/PatternTools/interface/CollinearFitAtTM.h"
0020
0021 #include "TrackingTools/GsfTracking/interface/TrajGsfTrackAssociation.h"
0022
0023 #include "TrackingTools/GsfTools/interface/GetComponents.h"
0024
0025 #include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"
0026
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0029
0030 void GsfTrackProducerBase::putInEvt(edm::Event& evt,
0031                                     const Propagator* prop,
0032                                     const MeasurementTracker* measTk,
0033                                     std::unique_ptr<TrackingRecHitCollection>& selHits,
0034                                     std::unique_ptr<reco::GsfTrackCollection>& selTracks,
0035                                     std::unique_ptr<reco::TrackExtraCollection>& selTrackExtras,
0036                                     std::unique_ptr<reco::GsfTrackExtraCollection>& selGsfTrackExtras,
0037                                     std::unique_ptr<std::vector<Trajectory> >& selTrajectories,
0038                                     AlgoProductCollection& algoResults,
0039                                     TransientTrackingRecHitBuilder const* hitBuilder,
0040                                     const reco::BeamSpot& bs,
0041                                     const TrackerTopology* ttopo) {
0042   TrackingRecHitRefProd rHits = evt.getRefBeforePut<TrackingRecHitCollection>();
0043   reco::TrackExtraRefProd rTrackExtras = evt.getRefBeforePut<reco::TrackExtraCollection>();
0044   reco::GsfTrackExtraRefProd rGsfTrackExtras = evt.getRefBeforePut<reco::GsfTrackExtraCollection>();
0045   reco::GsfTrackRefProd rTracks = evt.getRefBeforePut<reco::GsfTrackCollection>();
0046
0047   edm::Ref<reco::TrackExtraCollection>::key_type idx = 0;
0048   edm::Ref<reco::GsfTrackExtraCollection>::key_type idxGsf = 0;
0049   edm::Ref<reco::GsfTrackCollection>::key_type iTkRef = 0;
0050   edm::Ref<std::vector<Trajectory> >::key_type iTjRef = 0;
0051   std::map<unsigned int, unsigned int> tjTkMap;
0052
0053   TSCBLBuilderNoMaterial tscblBuilder;
0054
0055   for (auto& i : algoResults) {
0056     Trajectory* theTraj = i.trajectory;
0057     if (trajectoryInEvent_) {
0058       selTrajectories->push_back(*theTraj);
0059       iTjRef++;
0060     }
0061
0062     reco::GsfTrack* theTrack = i.track;
0063     PropagationDirection seedDir = i.pDir;
0064
0065     LogDebug("TrackProducer") << "In GsfTrackProducerBase::putInEvt - seedDir=" << seedDir;
0066
0067     reco::GsfTrack t = *theTrack;
0068     selTracks->push_back(t);
0069     iTkRef++;
0070
0071     // Store indices in local map (starts at 0)
0072     if (trajectoryInEvent_)
0073       tjTkMap[iTjRef - 1] = iTkRef - 1;
0074
0075     //sets the outermost and innermost TSOSs
0076     TrajectoryStateOnSurface outertsos;
0077     TrajectoryStateOnSurface innertsos;
0078     unsigned int innerId, outerId;
0079
0080     // ---  NOTA BENE: the convention is to sort hits and measurements "along the momentum".
0081     // This is consistent with innermost and outermost labels only for tracks from LHC collision
0082     if (theTraj->direction() == alongMomentum) {
0083       outertsos = theTraj->lastMeasurement().updatedState();
0084       innertsos = theTraj->firstMeasurement().updatedState();
0085       outerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
0086       innerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
0087     } else {
0088       outertsos = theTraj->firstMeasurement().updatedState();
0089       innertsos = theTraj->lastMeasurement().updatedState();
0090       outerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
0091       innerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
0092     }
0093     //build the TrackExtra
0094     GlobalPoint v = outertsos.globalParameters().position();
0095     GlobalVector p = outertsos.globalParameters().momentum();
0096     math::XYZVector outmom(p.x(), p.y(), p.z());
0097     math::XYZPoint outpos(v.x(), v.y(), v.z());
0098     v = innertsos.globalParameters().position();
0099     p = innertsos.globalParameters().momentum();
0100     math::XYZVector inmom(p.x(), p.y(), p.z());
0101     math::XYZPoint inpos(v.x(), v.y(), v.z());
0102
0103     reco::TrackExtraRef teref = reco::TrackExtraRef(rTrackExtras, idx++);
0104     reco::GsfTrack& track = selTracks->back();
0105     track.setExtra(teref);
0106
0107     //======= I want to set the second hitPattern here =============
0108     if (theSchool.isValid()) {
0109       edm::Handle<MeasurementTrackerEvent> mte;
0110       evt.getByToken(mteSrc_, mte);
0111       // NavigationSetter setter( *theSchool );
0112       setSecondHitPattern(theTraj, track, prop, &*mte, ttopo);
0113     }
0114     //==============================================================
0115
0116     selTrackExtras->push_back(reco::TrackExtra(outpos,
0117                                                outmom,
0118                                                true,
0119                                                inpos,
0120                                                inmom,
0121                                                true,
0122                                                outertsos.curvilinearError(),
0123                                                outerId,
0124                                                innertsos.curvilinearError(),
0125                                                innerId,
0126                                                seedDir,
0127                                                theTraj->seedRef()));
0128
0129     reco::TrackExtra& tx = selTrackExtras->back();
0130
0131     // ---  NOTA BENE: the convention is to sort hits and measurements "along the momentum".
0132     // This is consistent with innermost and outermost labels only for tracks from LHC collisions
0133     reco::TrackExtra::TrajParams trajParams;
0134     reco::TrackExtra::Chi2sFive chi2s;
0135     Traj2TrackHits t2t;
0136     auto ih = selHits->size();
0137     t2t(*theTraj, *selHits, trajParams, chi2s);
0138     auto ie = selHits->size();
0139     tx.setHits(rHits, ih, ie - ih);
0140     tx.setTrajParams(std::move(trajParams), std::move(chi2s));
0141     for (; ih < ie; ++ih) {
0142       auto const& hit = (*selHits)[ih];
0143       track.appendHitPattern(hit, *ttopo);
0144     }
0145
0146     // ----
0147
0148     std::vector<reco::GsfTangent> tangents;
0149     const Trajectory::DataContainer& measurements = theTraj->measurements();
0150     if (measurements.size() > 2) {
0151       tangents.reserve(measurements.size() - 2);
0152       Trajectory::DataContainer::const_iterator ibegin, iend;
0153       int increment(0);
0154       if (theTraj->direction() == alongMomentum) {
0155         ibegin = measurements.begin() + 1;
0156         iend = measurements.end() - 1;
0157         increment = 1;
0158       } else {
0159         ibegin = measurements.end() - 2;
0160         iend = measurements.begin();
0161         increment = -1;
0162       }
0163       math::XYZPoint position;
0164       math::XYZVector momentum;
0165       Measurement1D deltaP;
0166       // only measurements on "mono" detectors
0167       for (Trajectory::DataContainer::const_iterator i = ibegin; i != iend; i += increment) {
0168         if (i->recHit().get()) {
0169           DetId detId(i->recHit()->geographicalId());
0170           if (detId.det() == DetId::Tracker) {
0171             int subdetId = detId.subdetId();
0172             if (subdetId == SiStripDetId::TIB || subdetId == SiStripDetId::TID || subdetId == SiStripDetId::TOB ||
0173                 subdetId == SiStripDetId::TEC) {
0174               if (SiStripDetId(detId).stereo())
0175                 continue;
0176             }
0177           }
0178         }
0179         bool valid = computeModeAtTM(*i, position, momentum, deltaP);
0180         if (valid) {
0181           tangents.push_back(reco::GsfTangent(position, momentum, deltaP));
0182         }
0183       }
0184     }
0185
0186     //build the GsfTrackExtra
0187     std::vector<reco::GsfComponent5D> outerStates;
0188     fillStates(outertsos, outerStates);
0189     std::vector<reco::GsfComponent5D> innerStates;
0190     fillStates(innertsos, innerStates);
0191
0192     reco::GsfTrackExtraRef terefGsf = reco::GsfTrackExtraRef(rGsfTrackExtras, idxGsf++);
0193     track.setGsfExtra(terefGsf);
0194     selGsfTrackExtras->push_back(reco::GsfTrackExtra(
0195         outerStates, outertsos.localParameters().pzSign(), innerStates, innertsos.localParameters().pzSign(), tangents));
0196
0197     if (innertsos.isValid()) {
0198       GsfPropagatorAdapter gsfProp(AnalyticalPropagator(innertsos.magneticField(), anyDirection));
0199       TransverseImpactPointExtrapolator tipExtrapolator(gsfProp);
0200       fillMode(track, innertsos, gsfProp, tipExtrapolator, tscblBuilder, bs);
0201     }
0202
0203     delete theTrack;
0204     delete theTraj;
0205   }
0206
0207   LogTrace("TrackingRegressionTest") << "========== TrackProducer Info ===================";
0208   LogTrace("TrackingRegressionTest") << "number of finalGsfTracks: " << selTracks->size();
0209   for (reco::GsfTrackCollection::const_iterator it = selTracks->begin(); it != selTracks->end(); it++) {
0210     LogTrace("TrackingRegressionTest") << "track's n valid and invalid hit, chi2, pt : " << it->found() << " , "
0211                                        << it->lost() << " , " << it->normalizedChi2() << " , " << it->pt() << " , "
0212                                        << it->eta();
0213   }
0214   LogTrace("TrackingRegressionTest") << "=================================================";
0215
0216   rTracks_ = evt.put(std::move(selTracks));
0217   evt.put(std::move(selTrackExtras));
0218   evt.put(std::move(selGsfTrackExtras));
0219   evt.put(std::move(selHits));
0220
0221   if (trajectoryInEvent_) {
0222     edm::OrphanHandle<std::vector<Trajectory> > rTrajs = evt.put(std::move(selTrajectories));
0223
0224     // Now Create traj<->tracks association map
0225     std::unique_ptr<TrajGsfTrackAssociationCollection> trajTrackMap(
0226         new TrajGsfTrackAssociationCollection(rTrajs, rTracks_));
0227     for (std::map<unsigned int, unsigned int>::iterator i = tjTkMap.begin(); i != tjTkMap.end(); i++) {
0228       edm::Ref<std::vector<Trajectory> > trajRef(rTrajs, (*i).first);
0229       edm::Ref<reco::GsfTrackCollection> tkRef(rTracks_, (*i).second);
0230       trajTrackMap->insert(edm::Ref<std::vector<Trajectory> >(rTrajs, (*i).first),
0231                            edm::Ref<reco::GsfTrackCollection>(rTracks_, (*i).second));
0232     }
0233     evt.put(std::move(trajTrackMap));
0234   }
0235 }
0236
0237 void GsfTrackProducerBase::fillStates(TrajectoryStateOnSurface tsos, std::vector<reco::GsfComponent5D>& states) const {
0238   reco::GsfComponent5D::ParameterVector pLocS;
0239   reco::GsfComponent5D::CovarianceMatrix cLocS;
0240   GetComponents comps(tsos);
0241   auto const& components = comps();
0242   states.reserve(components.size());
0243   for (auto const& st : components)
0244     states.emplace_back(st.weight(), st.localParameters().vector(), st.localError().matrix());
0245 }
0246
0247 void GsfTrackProducerBase::fillMode(reco::GsfTrack& track,
0248                                     const TrajectoryStateOnSurface innertsos,
0249                                     const Propagator& gsfProp,
0250                                     const TransverseImpactPointExtrapolator& tipExtrapolator,
0251                                     TrajectoryStateClosestToBeamLineBuilder& tscblBuilder,
0252                                     const reco::BeamSpot& bs) const {
0253   // Get transverse impact parameter plane (from mean). This is a first approximation;
0254   // the mode is then extrapolated to the
0255   // final position closest to the beamline.
0256   GlobalPoint bsPos(bs.position().x() + (track.vz() - bs.position().z()) * bs.dxdz(),
0257                     bs.position().y() + (track.vz() - bs.position().z()) * bs.dydz(),
0258                     track.vz());
0259   TrajectoryStateOnSurface vtxTsos = tipExtrapolator.extrapolate(innertsos, bsPos);
0260   if (!vtxTsos.isValid())
0261     vtxTsos = innertsos;
0262   // extrapolate mixture
0263   vtxTsos = gsfProp.propagate(innertsos, vtxTsos.surface());
0264   if (!vtxTsos.isValid())
0265     return;  // failed (GsfTrack keeps mode = mean)
0266   // extract mode
0267   // build perigee parameters (for covariance to be stored)
0268   AlgebraicVector5 modeParameters;
0269   AlgebraicSymMatrix55 modeCovariance;
0270   // set parameters and variances for "mode" state (local parameters)
0271   for (unsigned int iv = 0; iv < 5; ++iv) {
0272     MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(vtxTsos, iv);
0273     GaussianSumUtilities1D utils(state1D);
0274     modeParameters(iv) = utils.mode().mean();
0275     modeCovariance(iv, iv) = utils.mode().variance();
0276     if (!utils.modeIsValid()) {
0277       // if mode calculation fails: use mean
0278       modeParameters(iv) = utils.mean();
0279       modeCovariance(iv, iv) = utils.variance();
0280     }
0281   }
0282   // complete covariance matrix
0283   // approximation: use correlations from mean
0284   const AlgebraicSymMatrix55& meanCovariance(vtxTsos.localError().matrix());
0285   for (unsigned int iv1 = 0; iv1 < 5; ++iv1) {
0286     for (unsigned int iv2 = 0; iv2 < iv1; ++iv2) {
0287       double cov12 = meanCovariance(iv1, iv2) * sqrt(modeCovariance(iv1, iv1) / meanCovariance(iv1, iv1) *
0288                                                      modeCovariance(iv2, iv2) / meanCovariance(iv2, iv2));
0289       modeCovariance(iv1, iv2) = modeCovariance(iv2, iv1) = cov12;
0290     }
0291   }
0292   TrajectoryStateOnSurface modeTsos(LocalTrajectoryParameters(modeParameters, vtxTsos.localParameters().pzSign()),
0293                                     LocalTrajectoryError(modeCovariance),
0294                                     vtxTsos.surface(),
0295                                     vtxTsos.magneticField(),
0296                                     vtxTsos.surfaceSide());
0297   TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*modeTsos.freeState(), bs);
0298   if (!tscbl.isValid())
0299     return;  // failed (GsfTrack keeps mode = mean)
0300   //
0301   // extract state at PCA and create momentum vector and covariance matrix
0302   //
0303   const FreeTrajectoryState& fts = tscbl.trackStateAtPCA();
0304   GlobalVector tscblMom = fts.momentum();
0305   reco::GsfTrack::Vector mom(tscblMom.x(), tscblMom.y(), tscblMom.z());
0306   reco::GsfTrack::CovarianceMatrixMode cov;
0307   const AlgebraicSymMatrix55& tscblCov = fts.curvilinearError().matrix();
0308   for (unsigned int iv1 = 0; iv1 < reco::GsfTrack::dimensionMode; ++iv1) {
0309     for (unsigned int iv2 = 0; iv2 < reco::GsfTrack::dimensionMode; ++iv2) {
0310       cov(iv1, iv2) = tscblCov(iv1, iv2);
0311     }
0312   }
0313   track.setMode(fts.charge(), mom, cov);
0314 }
0315
0316 void GsfTrackProducerBase::localParametersFromQpMode(const TrajectoryStateOnSurface tsos,
0317                                                      AlgebraicVector5& parameters,
0318                                                      AlgebraicSymMatrix55& covariance) const {
0319   //
0320   // parameters and errors from combined state
0321   //
0322   parameters = tsos.localParameters().vector();
0323   covariance = tsos.localError().matrix();
0324   //
0325   // mode for parameter 0 (q/p)
0326   //
0327   MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(tsos, 0));
0328   GaussianSumUtilities1D qpGS(qpState);
0329   if (!qpGS.modeIsValid())
0330     return;
0331   double qp = qpGS.mode().mean();
0332   double varQp = qpGS.mode().variance();
0333   //
0334   // replace q/p value and variance, rescale correlation terms
0335   //   (heuristic procedure - alternative would be mode in 5D ...)
0336   //
0337   double VarQpRatio = sqrt(varQp / covariance(0, 0));
0338   parameters(0) = qp;
0339   covariance(0, 0) = varQp;
0340   for (int i = 1; i < 5; ++i)
0341     covariance(i, 0) *= VarQpRatio;
0342 }
0343
0344 bool GsfTrackProducerBase::computeModeAtTM(const TrajectoryMeasurement& tm,
0345                                            reco::GsfTrackExtra::Point& position,
0346                                            reco::GsfTrackExtra::Vector& momentum,
0347                                            Measurement1D& deltaP) const {
0348   //
0349   // states
0350   //
0351   TrajectoryStateOnSurface const& fwdState = tm.forwardPredictedState();
0352   TrajectoryStateOnSurface const& bwdState = tm.backwardPredictedState();
0353   TrajectoryStateOnSurface const& upState = tm.updatedState();
0354
0355   if (!fwdState.isValid() || !bwdState.isValid() || !upState.isValid()) {
0356     return false;
0357   }
0358   //
0359   // position from mean, momentum from mode (in cartesian coordinates)
0360   //  following PF code
0361   //
0362   GlobalPoint pos = upState.globalPosition();
0363   position = reco::GsfTrackExtra::Point(pos.x(), pos.y(), pos.z());
0364   GlobalVector mom;
0365   bool result = multiTrajectoryStateMode::momentumFromModeCartesian(upState, mom);
0366   if (!result) {
0367     //     std::cout << "momentumFromModeCartesian failed" << std::endl;
0368     return false;
0369   }
0370   momentum = reco::GsfTrackExtra::Vector(mom.x(), mom.y(), mom.z());
0371   //
0372   // calculation from deltaP from fit to forward & backward predictions
0373   //  (momentum from mode) and hit
0374   //
0375   // prepare input parameter vectors and covariance matrices
0376   AlgebraicVector5 fwdPars = fwdState.localParameters().vector();
0377   AlgebraicSymMatrix55 fwdCov = fwdState.localError().matrix();
0378   localParametersFromQpMode(fwdState, fwdPars, fwdCov);
0379   AlgebraicVector5 bwdPars = bwdState.localParameters().vector();
0380   AlgebraicSymMatrix55 bwdCov = bwdState.localError().matrix();
0381   localParametersFromQpMode(bwdState, bwdPars, bwdCov);
0382   LocalPoint hitPos(0., 0., 0.);
0383   LocalError hitErr(-1., -1., -1.);
0384   if (tm.recHit()->isValid()) {
0385     hitPos = tm.recHit()->localPosition();
0386     hitErr = tm.recHit()->localPositionError();
0387   }
0388   CollinearFitAtTM2 collinearFit(fwdPars, fwdCov, bwdPars, bwdCov, hitPos, hitErr);
0389   deltaP = collinearFit.deltaP();
0390
0391   return true;
0392 }