ReferenceTrajectories/plugins/DualTrajectoryFactory.cc

0001 #include "FWCore/Framework/interface/ESHandle.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryPlugin.h"
0008 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0009
0010 #include <algorithm>
0011
0012 #include "Alignment/ReferenceTrajectories/interface/DualReferenceTrajectory.h"
0013 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryBase.h"
0014
0015 /// A factory that produces instances of class ReferenceTrajectory from a given TrajTrackPairCollection.
0016
0017 class DualTrajectoryFactory : public TrajectoryFactoryBase {
0018 public:
0019   DualTrajectoryFactory(const edm::ParameterSet &config, edm::ConsumesCollector &iC);
0020   ~DualTrajectoryFactory() override;
0021   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagFieldToken;
0022
0023   /// Produce the reference trajectories.
0024   const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
0025                                                    const ConstTrajTrackPairCollection &tracks,
0026                                                    const reco::BeamSpot &beamSpot) const override;
0027
0028   const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
0029                                                    const ConstTrajTrackPairCollection &tracks,
0030                                                    const ExternalPredictionCollection &external,
0031                                                    const reco::BeamSpot &beamSpot) const override;
0032
0033   DualTrajectoryFactory *clone() const override { return new DualTrajectoryFactory(*this); }
0034
0035 protected:
0036   struct DualTrajectoryInput {
0037     TrajectoryStateOnSurface refTsos;
0038     TransientTrackingRecHit::ConstRecHitContainer fwdRecHits;
0039     TransientTrackingRecHit::ConstRecHitContainer bwdRecHits;
0040   };
0041
0042   const DualTrajectoryInput referenceStateAndRecHits(const ConstTrajTrackPair &track) const;
0043
0044   const TrajectoryStateOnSurface propagateExternal(const TrajectoryStateOnSurface &external,
0045                                                    const Surface &surface,
0046                                                    const MagneticField *magField) const;
0047
0048   double theMass;
0049 };
0050
0051 ////////////////////////////////////////////////////////////////////////////////
0052 ////////////////////////////////////////////////////////////////////////////////
0053 ////////////////////////////////////////////////////////////////////////////////
0054
0055 DualTrajectoryFactory::DualTrajectoryFactory(const edm::ParameterSet &config, edm::ConsumesCollector &iC)
0056     : TrajectoryFactoryBase(config, iC),
0057       m_MagFieldToken(iC.esConsumes()),
0058       theMass(config.getParameter<double>("ParticleMass")) {
0059   edm::LogInfo("Alignment") << "@SUB=DualTrajectoryFactory"
0060                             << "mass: " << theMass;
0061 }
0062
0063 DualTrajectoryFactory::~DualTrajectoryFactory(void) {}
0064
0065 const DualTrajectoryFactory::ReferenceTrajectoryCollection DualTrajectoryFactory::trajectories(
0066     const edm::EventSetup &setup, const ConstTrajTrackPairCollection &tracks, const reco::BeamSpot &beamSpot) const {
0067   ReferenceTrajectoryCollection trajectories;
0068
0069   const MagneticField *magneticField = &setup.getData(m_MagFieldToken);
0070
0071   if (magneticField->inTesla(GlobalPoint(0., 0., 0.)).mag2() < 1.e-6) {
0072     edm::LogWarning("Alignment") << "@SUB=DualTrajectoryFactory::trajectories"
0073                                  << "B-field in z is " << magneticField->inTesla(GlobalPoint(0., 0., 0.)).z()
0074                                  << ": You should probably use the DualBzeroTrajectoryFactory\n"
0075                                  << "or fix this code to switch automatically as the ReferenceTrajectoryFactory does.";
0076   }
0077
0078   ConstTrajTrackPairCollection::const_iterator itTracks = tracks.begin();
0079
0080   while (itTracks != tracks.end()) {
0081     const DualTrajectoryInput input = this->referenceStateAndRecHits(*itTracks);
0082     // Check input: If all hits were rejected, the TSOS is initialized as invalid.
0083     if (input.refTsos.isValid()) {
0084       ReferenceTrajectoryBase::Config config(materialEffects(), propagationDirection(), theMass);
0085       config.useBeamSpot = useBeamSpot_;
0086       config.includeAPEs = includeAPEs_;
0087       config.allowZeroMaterial = allowZeroMaterial_;
0088       ReferenceTrajectoryPtr ptr(new DualReferenceTrajectory(
0089           input.refTsos, input.fwdRecHits, input.bwdRecHits, magneticField, beamSpot, config));
0090       trajectories.push_back(ptr);
0091     }
0092
0093     ++itTracks;
0094   }
0095
0096   return trajectories;
0097 }
0098
0099 const DualTrajectoryFactory::ReferenceTrajectoryCollection DualTrajectoryFactory::trajectories(
0100     const edm::EventSetup &setup,
0101     const ConstTrajTrackPairCollection &tracks,
0102     const ExternalPredictionCollection &external,
0103     const reco::BeamSpot &beamSpot) const {
0104   ReferenceTrajectoryCollection trajectories;
0105
0106   if (tracks.size() != external.size()) {
0107     edm::LogInfo("ReferenceTrajectories")
0108         << "@SUB=DualTrajectoryFactory::trajectories"
0109         << "Inconsistent input:\n"
0110         << "\tnumber of tracks = " << tracks.size() << "\tnumber of external predictions = " << external.size();
0111     return trajectories;
0112   }
0113   const MagneticField *magneticField = &setup.getData(m_MagFieldToken);
0114
0115   if (magneticField->inTesla(GlobalPoint(0., 0., 0.)).mag2() < 1.e-6) {
0116     edm::LogWarning("Alignment") << "@SUB=DualTrajectoryFactory::trajectories"
0117                                  << "B-field in z is " << magneticField->inTesla(GlobalPoint(0., 0., 0.)).z()
0118                                  << ": You should probably use the DualBzeroTrajectoryFactory\n"
0119                                  << "or fix this code to switch automatically as the ReferenceTrajectoryFactory does.";
0120   }
0121
0122   ConstTrajTrackPairCollection::const_iterator itTracks = tracks.begin();
0123   ExternalPredictionCollection::const_iterator itExternal = external.begin();
0124
0125   while (itTracks != tracks.end()) {
0126     const DualTrajectoryInput input = referenceStateAndRecHits(*itTracks);
0127     // Check input: If all hits were rejected, the TSOS is initialized as invalid.
0128     if (input.refTsos.isValid()) {
0129       if ((*itExternal).isValid()) {
0130         TrajectoryStateOnSurface propExternal = propagateExternal(*itExternal, input.refTsos.surface(), magneticField);
0131
0132         if (!propExternal.isValid())
0133           continue;
0134
0135         ReferenceTrajectoryBase::Config config(materialEffects(), propagationDirection(), theMass);
0136         config.useBeamSpot = useBeamSpot_;
0137         config.includeAPEs = includeAPEs_;
0138         config.allowZeroMaterial = allowZeroMaterial_;
0139         ReferenceTrajectoryPtr ptr(new DualReferenceTrajectory(
0140             propExternal, input.fwdRecHits, input.bwdRecHits, magneticField, beamSpot, config));
0141
0142         AlgebraicSymMatrix externalParamErrors(asHepMatrix<5>(propExternal.localError().matrix()));
0143         ptr->setParameterErrors(externalParamErrors);
0144         trajectories.push_back(ptr);
0145       } else {
0146         ReferenceTrajectoryBase::Config config(materialEffects(), propagationDirection(), theMass);
0147         config.useBeamSpot = useBeamSpot_;
0148         config.includeAPEs = includeAPEs_;
0149         config.allowZeroMaterial = allowZeroMaterial_;
0150         ReferenceTrajectoryPtr ptr(new DualReferenceTrajectory(
0151             input.refTsos, input.fwdRecHits, input.bwdRecHits, magneticField, beamSpot, config));
0152         trajectories.push_back(ptr);
0153       }
0154     }
0155
0156     ++itTracks;
0157     ++itExternal;
0158   }
0159
0160   return trajectories;
0161 }
0162
0163 const DualTrajectoryFactory::DualTrajectoryInput DualTrajectoryFactory::referenceStateAndRecHits(
0164     const ConstTrajTrackPair &track) const {
0165   DualTrajectoryInput input;
0166
0167   // get the trajectory measurements in the correct order, i.e. reverse if needed
0168   Trajectory::DataContainer allTrajMeas = this->orderedTrajectoryMeasurements(*track.first);
0169   Trajectory::DataContainer usedTrajMeas;
0170   Trajectory::DataContainer::iterator itM;
0171   // get all relevant trajectory measurements
0172   for (itM = allTrajMeas.begin(); itM != allTrajMeas.end(); itM++) {
0173     if (useRecHit((*itM).recHit()))
0174       usedTrajMeas.push_back(*itM);
0175   }
0176
0177   unsigned int iMeas = 0;
0178   unsigned int nMeas = usedTrajMeas.size();
0179   unsigned int nRefStateMeas = nMeas / 2;
0180   // get the valid RecHits
0181   for (itM = usedTrajMeas.begin(); itM != usedTrajMeas.end(); itM++, iMeas++) {
0182     TransientTrackingRecHit::ConstRecHitPointer aRecHit = (*itM).recHit();
0183
0184     if (iMeas < nRefStateMeas) {
0185       input.bwdRecHits.push_back(aRecHit);
0186     } else if (iMeas > nRefStateMeas) {
0187       input.fwdRecHits.push_back(aRecHit);
0188     } else {  // iMeas == nRefStateMeas
0189       if ((*itM).updatedState().isValid()) {
0190         input.refTsos = (*itM).updatedState();
0191         input.bwdRecHits.push_back(aRecHit);
0192         input.fwdRecHits.push_back(aRecHit);
0193       } else {
0194         // if the tsos of the middle hit is not valid, try the next one ...
0195         nRefStateMeas++;
0196         input.bwdRecHits.push_back(aRecHit);
0197       }
0198     }
0199   }
0200
0201   // bring input.fwdRecHits into correct order
0202   std::reverse(input.bwdRecHits.begin(), input.bwdRecHits.end());
0203
0204   return input;
0205 }
0206
0207 const TrajectoryStateOnSurface DualTrajectoryFactory::propagateExternal(const TrajectoryStateOnSurface &external,
0208                                                                         const Surface &surface,
0209                                                                         const MagneticField *magField) const {
0210   AnalyticalPropagator propagator(magField, anyDirection);
0211   const std::pair<TrajectoryStateOnSurface, double> tsosWithPath = propagator.propagateWithPath(external, surface);
0212   return tsosWithPath.first;
0213 }
0214
0215 DEFINE_EDM_PLUGIN(TrajectoryFactoryPlugin, DualTrajectoryFactory, "DualTrajectoryFactory");