TrackProducer/interface/TrackProducerBase.icc

0001 #include "RecoTracker/TrackProducer/interface/TrackProducerBase.h"
0002
0003 /// user include files

0004 #include "DataFormats/Common/interface/AssociationMap.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/EmptyGroupDescription.h"
0010 #include "FWCore/Utilities/interface/ESInputTag.h"
0011 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0012 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0013 #include "MagneticField/Engine/interface/MagneticField.h"
0014 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0015 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0016 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0017 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0018 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
0019 #include "TrackingTools/DetLayers/interface/GeometricSearchDet.h"
0020 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0021 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0022 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0023 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0024 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0025 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0026 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0027 #include "TrackingTools/TrackFitters/interface/TrajectoryFitterRecord.h"
0028
0029 //destructor

0030 template <class T>
0031 TrackProducerBase<T>::~TrackProducerBase() noexcept(false) {}
0032
0033 // member functions

0034
0035 template <class T>
0036 void TrackProducerBase<T>::initTrackProducerBase(const edm::ParameterSet& conf,
0037                                                  edm::ConsumesCollector cc,
0038                                                  const edm::EDGetToken& src) {
0039   conf_ = conf;
0040   src_ = src;
0041   bsSrc_ = cc.consumes(conf.getParameter<edm::InputTag>("beamSpot"));
0042   mteSrc_ = cc.consumes(conf.getParameter<edm::InputTag>("MeasurementTrackerEvent"));
0043
0044   trackGeomSrc_ = cc.esConsumes();
0045
0046   // 2014/02/11 mia:

0047   // we should get rid of the boolean parameter useSimpleMF,

0048   // and use only a string magneticField [instead of SimpleMagneticField]

0049   // or better an edm::ESInputTag (at the moment HLT does not handle ESInputTag)

0050   bool useSimpleMF = conf_.getParameter<bool>("useSimpleMF");
0051   std::string mfName = "";
0052   if (useSimpleMF) {
0053     mfName = conf_.getParameter<std::string>("SimpleMagneticField");
0054   }
0055   mfSrc_ = cc.esConsumes(edm::ESInputTag("", mfName));
0056
0057   std::string fitterName = conf_.getParameter<std::string>("Fitter");
0058   fitterSrc_ = cc.esConsumes(edm::ESInputTag("", fitterName));
0059
0060   std::string propagatorName = conf_.getParameter<std::string>("Propagator");
0061   propagatorSrc_ = cc.esConsumes(edm::ESInputTag("", propagatorName));
0062
0063   std::string builderName = conf_.getParameter<std::string>("TTRHBuilder");
0064   builderSrc_ = cc.esConsumes(edm::ESInputTag("", builderName));
0065
0066   //

0067   // get also the measurementTracker and the NavigationSchool

0068   // (they are necessary to fill in the secondary hit patterns)

0069   //

0070
0071   std::string theNavigationSchool = conf_.getParameter<std::string>("NavigationSchool");
0072   if (!theNavigationSchool.empty()) {
0073     useSchool_ = true;
0074     schoolSrc_ = cc.esConsumes(edm::ESInputTag("", theNavigationSchool));
0075     std::string measTkName = conf_.getParameter<std::string>("MeasurementTracker");
0076     measTkSrc_ = cc.esConsumes(edm::ESInputTag("", measTkName));
0077   }
0078 }
0079
0080 template <class T>
0081 void TrackProducerBase<T>::fillPSetDescription(edm::ParameterSetDescription& desc) {
0082   desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0083   desc.add<edm::InputTag>("MeasurementTrackerEvent", edm::InputTag("MeasurementTrackerEvent"));
0084   desc.add<bool>("useSimpleMF", false);
0085   desc.add<std::string>("SimpleMagneticField", "");
0086   desc.add<std::string>("Fitter", "KFFittingSmootherWithOutliersRejectionAndRK");
0087   desc.add<std::string>("Propagator", "RungeKuttaTrackerPropagator");
0088   desc.add<std::string>("TTRHBuilder", "WithAngleAndTemplate");
0089   desc.add<std::string>("NavigationSchool", "SimpleNavigationSchool");
0090   desc.add<std::string>("MeasurementTracker", "");
0091 }
0092
0093 // ------------ method called to produce the data  ------------

0094 template <class T>
0095 void TrackProducerBase<T>::getFromES(const edm::EventSetup& setup,
0096                                      edm::ESHandle<TrackerGeometry>& theG,
0097                                      edm::ESHandle<MagneticField>& theMF,
0098                                      edm::ESHandle<TrajectoryFitter>& theFitter,
0099                                      edm::ESHandle<Propagator>& thePropagator,
0100                                      edm::ESHandle<MeasurementTracker>& theMeasTk,
0101                                      edm::ESHandle<TransientTrackingRecHitBuilder>& theBuilder) {
0102   //

0103   //get geometry

0104   //

0105   LogDebug("TrackProducer") << "get geometry"
0106                             << "\n";
0107   theG = setup.getHandle(trackGeomSrc_);
0108   //

0109   //get magnetic field

0110   //

0111   LogDebug("TrackProducer") << "get magnetic field"
0112                             << "\n";
0113   theMF = setup.getHandle(mfSrc_);
0114
0115   //

0116   // get the fitter from the ES

0117   //

0118   LogDebug("TrackProducer") << "get the fitter from the ES"
0119                             << "\n";
0120   theFitter = setup.getHandle(fitterSrc_);
0121   //

0122   // get also the propagator

0123   //

0124   LogDebug("TrackProducer") << "get also the propagator"
0125                             << "\n";
0126   thePropagator = setup.getHandle(propagatorSrc_);
0127
0128   //

0129   // get the builder

0130   //

0131   LogDebug("TrackProducer") << "get also the TransientTrackingRecHitBuilder"
0132                             << "\n";
0133   theBuilder = setup.getHandle(builderSrc_);
0134
0135   //

0136   // get also the measurementTracker and the NavigationSchool

0137   // (they are necessary to fill in the secondary hit patterns)

0138   //

0139
0140   if (useSchool_) {
0141     LogDebug("TrackProducer") << "get a navigation school";
0142     theSchool = setup.getHandle(schoolSrc_);
0143     LogDebug("TrackProducer") << "get also the measTk"
0144                               << "\n";
0145     theMeasTk = setup.getHandle(measTkSrc_);
0146   } else {
0147     theSchool = edm::ESHandle<NavigationSchool>();    //put an invalid handle

0148     theMeasTk = edm::ESHandle<MeasurementTracker>();  //put an invalid handle

0149   }
0150 }
0151
0152 template <class T>
0153 void TrackProducerBase<T>::getFromEvt(edm::Event& theEvent,
0154                                       edm::Handle<TrackCandidateCollection>& theTCCollection,
0155                                       reco::BeamSpot& bs) {
0156   //

0157   //get the TrackCandidateCollection from the event

0158   //

0159   LogDebug("TrackProducer") << "get the TrackCandidateCollection from the event, source is "
0160                             << conf_.getParameter<edm::InputTag>("src") << "\n";
0161   theEvent.getByToken(src_, theTCCollection);
0162
0163   //get the BeamSpot

0164   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0165   theEvent.getByToken(bsSrc_, recoBeamSpotHandle);
0166   if (recoBeamSpotHandle.isValid())
0167     bs = *recoBeamSpotHandle;
0168   else
0169     edm::LogWarning("TrackProducerBase") << " BeamSpot is not valid";
0170 }
0171
0172 template <class T>
0173 void TrackProducerBase<T>::getFromEvt(edm::Event& theEvent,
0174                                       edm::Handle<TrackView>& theTCollection,
0175                                       reco::BeamSpot& bs) {
0176   //

0177   //get the TrackCollection from the event

0178   //

0179   LogDebug("TrackProducer") << "get the TrackCollection from the event, source is "
0180                             << conf_.getParameter<edm::InputTag>("src") << "\n";
0181   theEvent.getByToken(src_, theTCollection);
0182
0183   //get the BeamSpot

0184   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0185   theEvent.getByToken(bsSrc_, recoBeamSpotHandle);
0186   if (recoBeamSpotHandle.isValid())
0187     bs = *recoBeamSpotHandle;
0188   else
0189     edm::LogWarning("TrackProducerBase") << " BeamSpot is not valid";
0190 }
0191
0192 #include <TrackingTools/DetLayers/interface/DetLayer.h>
0193 #include <DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h>
0194
0195 template <class T>
0196 void TrackProducerBase<T>::setSecondHitPattern(Trajectory* traj,
0197                                                T& track,
0198                                                const Propagator* prop,
0199                                                const MeasurementTrackerEvent* measTk,
0200                                                const TrackerTopology* ttopo) {
0201   using namespace std;
0202   /// have to clone the propagator in order to change its propagation direction.

0203   std::unique_ptr<Propagator> localProp(prop->clone());
0204
0205   //use negative sigma=-3.0 in order to use a more conservative definition of isInside() for Bounds classes.

0206   Chi2MeasurementEstimator estimator(30., -3.0, 0.5, 2.0, 0.5, 1.e12);  // same as defauts....

0207
0208   // WARNING: At the moment the trajectories has the measurements with reversed sorting after the track smoothing.

0209   // Therefore the lastMeasurement is the inner one (for LHC-like tracks)

0210   if (!traj->firstMeasurement().updatedState().isValid() || !traj->lastMeasurement().updatedState().isValid())
0211     return;
0212
0213   const FreeTrajectoryState* outerState = traj->firstMeasurement().updatedState().freeState();
0214   const FreeTrajectoryState* innerState = traj->lastMeasurement().updatedState().freeState();
0215   TrajectoryStateOnSurface const& outerTSOS = traj->firstMeasurement().updatedState();
0216   TrajectoryStateOnSurface const& innerTSOS = traj->lastMeasurement().updatedState();
0217   const DetLayer* outerLayer = traj->firstMeasurement().layer();
0218   const DetLayer* innerLayer = traj->lastMeasurement().layer();
0219
0220   if (!outerLayer || !innerLayer) {
0221     //means  that the trajectory was fit/smoothed in a special case: not setting those pointers

0222     edm::LogError("TrackProducer") << "the trajectory was fit/smoothed in a special case: not setting those pointers.\n"
0223                                    << " Filling the secondary hit patterns was requested. So I will bail out.";
0224     throw cms::Exception("TrackProducerBase")
0225         << "the trajectory was fit/smoothed in a special case: not setting those pointers.\n"
0226         << " Filling the secondary hit patterns was requested. So I will bail out.";
0227   }
0228
0229   //WARNING: we are assuming that the hits were originally sorted along momentum (and therefore oppositeToMomentum after smoothing)

0230   PropagationDirection dirForInnerLayers = oppositeToMomentum;
0231   PropagationDirection dirForOuterLayers = alongMomentum;
0232   auto outIn = traj->direction() != oppositeToMomentum;
0233   if (outIn) {
0234     dirForInnerLayers = alongMomentum;
0235     dirForOuterLayers = oppositeToMomentum;
0236     // std::cout << "Iin setSecondHitPattern() logic. Trajectory direction (after smoothing) was not oppositeToMomentum. .. algo= " << track.algo() << std::endl;

0237   }
0238   // ----------- this previous block of code is not very safe. It should rely less on the sorting of the trajectory measurements -----

0239
0240   // Now all code looks as InOut in particular names

0241   // we will take care of OutIn only where it matters (MISSING_INNER vs _OUTER)

0242
0243   LogDebug("TrackProducer") << "calling inner compLayers()...";
0244   auto const& innerCompLayers = (*theSchool).compatibleLayers(*innerLayer, *innerState, dirForInnerLayers);
0245   LogDebug("TrackProducer") << "calling outer compLayers()...";
0246   auto const& outerCompLayers = (*theSchool).compatibleLayers(*outerLayer, *outerState, dirForOuterLayers);
0247
0248   LogDebug("TrackProducer") << "inner DetLayer  sub: " << innerLayer->subDetector() << "\n"
0249                             << "outer DetLayer  sub: " << outerLayer->subDetector() << "\n"
0250                             << "innerstate position rho: " << innerState->position().perp()
0251                             << " z: " << innerState->position().z() << "\n"
0252                             << "innerstate state pT: " << innerState->momentum().perp()
0253                             << " pz: " << innerState->momentum().z() << "\n"
0254                             << "outerstate position rho: " << outerState->position().perp()
0255                             << " z: " << outerState->position().z() << "\n"
0256                             << "outerstate state pT: " << outerState->momentum().perp()
0257                             << " pz: " << outerState->momentum().z() << "\n"
0258
0259                             << "innerLayers: " << innerCompLayers.size() << "\n"
0260                             << "outerLayers: " << outerCompLayers.size() << "\n";
0261
0262   auto loopOverLayer =
0263       [&](decltype(innerCompLayers) compLayers, TrajectoryStateOnSurface const& tsos) {
0264         for (auto it : compLayers) {
0265           if (it->basicComponents().empty()) {
0266             //this should never happen. but better protect for it

0267             edm::LogWarning("TrackProducer")
0268                 << "a detlayer with no components: I cannot figure out a DetId from this layer. please investigate.";
0269             continue;
0270           }
0271           auto const& detWithState = it->compatibleDets(tsos, *localProp, estimator);
0272           if (detWithState.empty())
0273             continue;
0274           DetId id = detWithState.front().first->geographicalId();
0275           MeasurementDetWithData const& measDet = measTk->idToDet(id);
0276           if (measDet.isActive() && !measDet.hasBadComponents(detWithState.front().second)) {
0277             InvalidTrackingRecHit tmpHit(*detWithState.front().first,
0278                                          outIn ? TrackingRecHit::missing_outer : TrackingRecHit::missing_inner);
0279             track.appendHitPattern(tmpHit, *ttopo);
0280             //cout << "WARNING: this hit is marked as lost because the detector was marked as active" << endl;

0281           } else {
0282             //cout << "WARNING: this hit is NOT marked as lost because the detector was marked as inactive" << endl;

0283             InvalidTrackingRecHit tmpHit(*detWithState.front().first,
0284                                          outIn ? TrackingRecHit::inactive_outer : TrackingRecHit::inactive_inner);
0285             track.appendHitPattern(tmpHit, *ttopo);
0286           }
0287         }  //end loop over layers

0288       };  // end lambda

0289
0290   localProp->setPropagationDirection(oppositeToMomentum);
0291   loopOverLayer(innerCompLayers, innerTSOS);
0292
0293   localProp->setPropagationDirection(alongMomentum);
0294   outIn = !outIn;  // if inOut should mark missing_outer

0295   loopOverLayer(outerCompLayers, outerTSOS);
0296 }