TkSeedGenerator/plugins/SeedFromConsecutiveHitsCreator.cc

0001 #include "SeedFromConsecutiveHitsCreator.h"
0002 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0003
0004 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0005 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/ESInputTag.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0013 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h"
0014 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0015 #include "FWCore/Utilities/interface/Likely.h"
0016
0017 namespace {
0018
0019   template <class T>
0020   inline T sqr(T t) {
0021     return t * t;
0022   }
0023
0024 }  // namespace
0025
0026 SeedFromConsecutiveHitsCreator::SeedFromConsecutiveHitsCreator(const edm::ParameterSet& cfg,
0027                                                                edm::ConsumesCollector&& iC)
0028     : thePropagatorLabel(cfg.getParameter<std::string>("propagator")),
0029       theBOFFMomentum(cfg.getParameter<double>("SeedMomentumForBOFF")),
0030       theOriginTransverseErrorMultiplier(cfg.getParameter<double>("OriginTransverseErrorMultiplier")),
0031       theMinOneOverPtError(cfg.getParameter<double>("MinOneOverPtError")),
0032       TTRHBuilder(cfg.getParameter<std::string>("TTRHBuilder")),
0033       mfName_(cfg.getParameter<std::string>("magneticField")),
0034       forceKinematicWithRegionDirection_(cfg.getParameter<bool>("forceKinematicWithRegionDirection")),
0035       trackerGeometryESToken_(iC.esConsumes()),
0036       propagatorESToken_(iC.esConsumes(edm::ESInputTag("", thePropagatorLabel))),
0037       magneticFieldESToken_(iC.esConsumes(edm::ESInputTag("", mfName_))),
0038       transientTrackingRecHitBuilderESToken_(iC.esConsumes(edm::ESInputTag("", TTRHBuilder))) {}
0039
0040 SeedFromConsecutiveHitsCreator::~SeedFromConsecutiveHitsCreator() {}
0041
0042 void SeedFromConsecutiveHitsCreator::fillDescriptions(edm::ParameterSetDescription& desc) {
0043   desc.add<std::string>("propagator", "PropagatorWithMaterialParabolicMf");
0044   desc.add<double>("SeedMomentumForBOFF", 5.0);
0045   desc.add<double>("OriginTransverseErrorMultiplier", 1.0);
0046   desc.add<double>("MinOneOverPtError", 1.0);
0047   desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0048   desc.add<std::string>("magneticField", "ParabolicMf");
0049   desc.add<bool>("forceKinematicWithRegionDirection", false);
0050 }
0051
0052 void SeedFromConsecutiveHitsCreator::init(const TrackingRegion& iregion,
0053                                           const edm::EventSetup& es,
0054                                           const SeedComparitor* ifilter) {
0055   region = &iregion;
0056   filter = ifilter;
0057   trackerGeometry_ = &es.getData(trackerGeometryESToken_);
0058   propagator_ = &es.getData(propagatorESToken_);
0059   magneticField_ = &es.getData(magneticFieldESToken_);
0060   nomField = magneticField_->nominalValue();
0061   isBOFF = (0 == nomField);
0062
0063   TransientTrackingRecHitBuilder const* transientTrackingRecHitBuilder =
0064       &es.getData(transientTrackingRecHitBuilderESToken_);
0065   auto builder = (TkTransientTrackingRecHitBuilder const*)(transientTrackingRecHitBuilder);
0066   cloner = (*builder).cloner();
0067 }
0068
0069 void SeedFromConsecutiveHitsCreator::makeSeed(TrajectorySeedCollection& seedCollection, const SeedingHitSet& hits) {
0070   if (hits.size() < 2)
0071     return;
0072
0073   GlobalTrajectoryParameters kine;
0074   if (!initialKinematic(kine, hits))
0075     return;
0076
0077   auto sin2Theta = kine.momentum().perp2() / kine.momentum().mag2();
0078
0079   CurvilinearTrajectoryError error = initialError(sin2Theta);
0080   FreeTrajectoryState fts(kine, error);
0081
0082   if (region->direction().x() != 0 &&
0083       forceKinematicWithRegionDirection_)  // a direction was given, check if it is an etaPhi region
0084   {
0085     const RectangularEtaPhiTrackingRegion* etaPhiRegion = dynamic_cast<const RectangularEtaPhiTrackingRegion*>(region);
0086     if (etaPhiRegion) {
0087       //the following completely reset the kinematics, perhaps it makes no sense and newKine=kine would do better
0088       GlobalVector direction = region->direction() / region->direction().mag();
0089       GlobalVector momentum = direction * fts.momentum().mag();
0090       GlobalPoint position = region->origin() + 5 * direction;
0091       GlobalTrajectoryParameters newKine(position, momentum, fts.charge(), &fts.parameters().magneticField());
0092
0093       auto ptMin = region->ptMin();
0094       CurvilinearTrajectoryError newError;  //zeroed
0095       auto& C = newError.matrix();
0096       constexpr float minC00 = 0.4f;
0097       C[0][0] = std::max(sin2Theta / sqr(ptMin), minC00);
0098       auto zErr = sqr(region->originZBound());
0099       auto transverseErr = sqr(region->originRBound());  // assume equal cxx cyy
0100       auto twiceDeltaLambda =
0101           std::atan(etaPhiRegion->tanLambdaRange().first) - std::atan(etaPhiRegion->tanLambdaRange().second);
0102       auto twiceDeltaPhi = etaPhiRegion->phiMargin().right() + etaPhiRegion->phiMargin().left();
0103       C[1][1] = twiceDeltaLambda * twiceDeltaLambda;  //2 sigma of what given in input
0104       C[2][2] = twiceDeltaPhi * twiceDeltaPhi;
0105       C[3][3] = transverseErr;
0106       C[4][4] = zErr * sin2Theta + transverseErr * (1.f - sin2Theta);
0107       fts = FreeTrajectoryState(newKine, newError);
0108     }
0109   }
0110
0111   buildSeed(seedCollection, hits, fts);
0112 }
0113
0114 bool SeedFromConsecutiveHitsCreator::initialKinematic(GlobalTrajectoryParameters& kine,
0115                                                       const SeedingHitSet& hits) const {
0116   SeedingHitSet::ConstRecHitPointer tth1 = hits[0];
0117   SeedingHitSet::ConstRecHitPointer tth2 = hits[1];
0118
0119   const GlobalPoint& vertexPos = region->origin();
0120
0121   FastHelix helix(tth2->globalPosition(), tth1->globalPosition(), vertexPos, nomField, magneticField_);
0122   if (helix.isValid()) {
0123     kine = helix.stateAtVertex();
0124   } else {
0125     GlobalVector initMomentum(tth2->globalPosition() - vertexPos);
0126     initMomentum *= (100.f / initMomentum.perp());
0127     kine = GlobalTrajectoryParameters(vertexPos, initMomentum, 1, magneticField_);
0128   }
0129
0130   if UNLIKELY (isBOFF && (theBOFFMomentum > 0)) {
0131     kine = GlobalTrajectoryParameters(
0132         kine.position(), kine.momentum().unit() * theBOFFMomentum, kine.charge(), magneticField_);
0133   }
0134   return (filter ? filter->compatible(hits, kine, helix) : true);
0135 }
0136
0137 CurvilinearTrajectoryError SeedFromConsecutiveHitsCreator::initialError(float sin2Theta) const {
0138   // Set initial uncertainty on track parameters, using only P.V. constraint and no hit
0139   // information.
0140   CurvilinearTrajectoryError newError;  // zeroed
0141   auto& C = newError.matrix();
0142
0143   // FIXME: minC00. Prevent apriori uncertainty in 1/P from being too small,
0144   // to avoid instabilities.
0145   // N.B. This parameter needs optimising ...
0146   // Probably OK based on quick study: KS 22/11/12.
0147   auto sin2th = sin2Theta;
0148   auto minC00 = sqr(theMinOneOverPtError);
0149   C[0][0] = std::max(sin2th / sqr(region->ptMin()), minC00);
0150   auto zErr = sqr(region->originZBound());
0151   auto transverseErr = sqr(theOriginTransverseErrorMultiplier * region->originRBound());
0152   C[1][1] = C[2][2] = 1.;  // no good reason. no bad reason....
0153   C[3][3] = transverseErr;
0154   C[4][4] = zErr * sin2th + transverseErr * (1.f - sin2th);
0155
0156   return newError;
0157 }
0158
0159 void SeedFromConsecutiveHitsCreator::buildSeed(TrajectorySeedCollection& seedCollection,
0160                                                const SeedingHitSet& hits,
0161                                                const FreeTrajectoryState& fts) const {
0162   // get updator
0163   KFUpdator updator;
0164
0165   // Now update initial state track using information from seed hits.
0166
0167   TrajectoryStateOnSurface updatedState;
0168   edm::OwnVector<TrackingRecHit> seedHits;
0169
0170   const TrackingRecHit* hit = nullptr;
0171   for (unsigned int iHit = 0; iHit < hits.size(); iHit++) {
0172     hit = hits[iHit]->hit();
0173     TrajectoryStateOnSurface state =
0174         (iHit == 0) ? propagator_->propagate(fts, trackerGeometry_->idToDet(hit->geographicalId())->surface())
0175                     : propagator_->propagate(updatedState, trackerGeometry_->idToDet(hit->geographicalId())->surface());
0176     if (!state.isValid())
0177       return;
0178
0179     SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
0180
0181     std::unique_ptr<BaseTrackerRecHit> newtth(refitHit(tth, state));
0182
0183     if (!checkHit(state, &*newtth))
0184       return;
0185
0186     updatedState = updator.update(state, *newtth);
0187     if (!updatedState.isValid())
0188       return;
0189
0190     seedHits.push_back(newtth.release());
0191   }
0192
0193   if (!hit)
0194     return;
0195
0196   PTrajectoryStateOnDet const& PTraj =
0197       trajectoryStateTransform::persistentState(updatedState, hit->geographicalId().rawId());
0198   seedCollection.emplace_back(PTraj, std::move(seedHits), alongMomentum);
0199 }
0200
0201 SeedingHitSet::RecHitPointer SeedFromConsecutiveHitsCreator::refitHit(SeedingHitSet::ConstRecHitPointer hit,
0202                                                                       const TrajectoryStateOnSurface& state) const {
0203   return (SeedingHitSet::RecHitPointer)(cloner(*hit, state));
0204 }
0205
0206 bool SeedFromConsecutiveHitsCreator::checkHit(const TrajectoryStateOnSurface& tsos,
0207                                               SeedingHitSet::ConstRecHitPointer hit) const {
0208   return (filter ? filter->compatible(tsos, hit) : true);
0209 }