TrackProducer/interface/TrackProducerBase.icc

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

0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include <FWCore/Utilities/interface/ESInputTag.h>
0008
0009 #include "MagneticField/Engine/interface/MagneticField.h"
0010 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0011 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0012 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0013 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0014 #include "TrackingTools/TrackFitters/interface/TrajectoryFitterRecord.h"
0015 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0016
0017 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0018 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0019
0020 #include "DataFormats/Common/interface/AssociationMap.h"
0021
0022 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0023 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
0024 // #include "TrackingTools/DetLayers/interface/NavigationSetter.h"

0025
0026 #include <TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h>
0027 #include <TrackingTools/DetLayers/interface/GeometricSearchDet.h>
0028 #include <RecoTracker/MeasurementDet/interface/MeasurementTracker.h>
0029 #include <RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h>
0030 #include <TrackingTools/MeasurementDet/interface/MeasurementDet.h>
0031 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0032
0033 #include "DataFormats/DetId/interface/DetId.h"
0034
0035 //destructor

0036 template <class T>
0037 TrackProducerBase<T>::~TrackProducerBase() noexcept(false) {}
0038
0039 // member functions

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

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

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

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

0056   bool useSimpleMF = false;
0057   if (conf_.exists("useSimpleMF"))
0058     useSimpleMF = conf_.getParameter<bool>("useSimpleMF");
0059   std::string mfName = "";
0060   if (useSimpleMF) {
0061     mfName = conf_.getParameter<std::string>("SimpleMagneticField");
0062   }
0063   mfSrc_ = cc.esConsumes(edm::ESInputTag("", mfName));
0064
0065   std::string fitterName = conf_.getParameter<std::string>("Fitter");
0066   fitterSrc_ = cc.esConsumes(edm::ESInputTag("", fitterName));
0067
0068   std::string propagatorName = conf_.getParameter<std::string>("Propagator");
0069   propagatorSrc_ = cc.esConsumes(edm::ESInputTag("", propagatorName));
0070
0071   std::string builderName = conf_.getParameter<std::string>("TTRHBuilder");
0072   builderSrc_ = cc.esConsumes(edm::ESInputTag("", builderName));
0073
0074   //

0075   // get also the measurementTracker and the NavigationSchool

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

0077   //

0078
0079   std::string theNavigationSchool = "";
0080   if (conf_.exists("NavigationSchool"))
0081     theNavigationSchool = conf_.getParameter<std::string>("NavigationSchool");
0082   else
0083     edm::LogWarning("TrackProducerBase")
0084         << " NavigationSchool parameter not set. secondary hit pattern will not be filled.";
0085   if (!theNavigationSchool.empty()) {
0086     useSchool_ = true;
0087     schoolSrc_ = cc.esConsumes(edm::ESInputTag("", theNavigationSchool));
0088     std::string measTkName = conf_.getParameter<std::string>("MeasurementTracker");
0089     measTkSrc_ = cc.esConsumes(edm::ESInputTag("", measTkName));
0090   }
0091 }
0092
0093 // ------------ method called to produce the data  ------------

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

0104   //get geometry

0105   //

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

0110   //get magnetic field

0111   //

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

0117   // get the fitter from the ES

0118   //

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

0123   // get also the propagator

0124   //

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

0130   // get the builder

0131   //

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

0137   // get also the measurementTracker and the NavigationSchool

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

0139   //

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

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

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

0158   //get the TrackCandidateCollection from the event

0159   //

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

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

0178   //get the TrackCollection from the event

0179   //

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

0289       };  // end lambda

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

0296   loopOverLayer(outerCompLayers, outerTSOS);
0297 }