File indexing completed on 2024-04-06 12:29:00
0001 #include "RecoTracker/TrackProducer/interface/GsfTrackProducerBase.h"
0002
0003 #include <memory>
0004
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
0072 if (trajectoryInEvent_)
0073 tjTkMap[iTjRef - 1] = iTkRef - 1;
0074
0075
0076 TrajectoryStateOnSurface outertsos;
0077 TrajectoryStateOnSurface innertsos;
0078 unsigned int innerId, outerId;
0079
0080
0081
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
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
0108 if (theSchool.isValid()) {
0109 edm::Handle<MeasurementTrackerEvent> mte;
0110 evt.getByToken(mteSrc_, mte);
0111
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
0132
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
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
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
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
0254
0255
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
0263 vtxTsos = gsfProp.propagate(innertsos, vtxTsos.surface());
0264 if (!vtxTsos.isValid())
0265 return;
0266
0267
0268 AlgebraicVector5 modeParameters;
0269 AlgebraicSymMatrix55 modeCovariance;
0270
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
0278 modeParameters(iv) = utils.mean();
0279 modeCovariance(iv, iv) = utils.variance();
0280 }
0281 }
0282
0283
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;
0300
0301
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
0321
0322 parameters = tsos.localParameters().vector();
0323 covariance = tsos.localError().matrix();
0324
0325
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
0335
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
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
0360
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
0368 return false;
0369 }
0370 momentum = reco::GsfTrackExtra::Vector(mom.x(), mom.y(), mom.z());
0371
0372
0373
0374
0375
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 }