ReferenceTrajectories/plugins/TwoBodyDecayTrajectoryFactory.cc

0001 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0002 #include "DataFormats/GeometrySurface/interface/Surface.h"
0003 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0004 #include "DataFormats/Math/interface/Vector.h"
0005 #include "DataFormats/Math/interface/Error.h"
0006 #include "RecoVertex/VertexTools/interface/PerigeeLinearizedTrackState.h"
0007 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryPlugin.h"
0008
0009 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0010 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0011 #include "Geometry/CommonDetUnit/interface/TrackerGeomDet.h"
0012
0013 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0014
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include "Alignment/ReferenceTrajectories/interface/TrajectoryFactoryBase.h"
0019 #include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectory.h"
0020
0021 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayFitter.h"
0022
0023 /**
0024    by Edmund Widl, see CMS NOTE-2007/032.
0025    extension for BreakPoints or BrokenLines by Claus Kleinwort
0026  */
0027
0028 class TwoBodyDecayTrajectoryFactory : public TrajectoryFactoryBase {
0029 public:
0030   typedef TwoBodyDecayVirtualMeasurement VirtualMeasurement;
0031   typedef TwoBodyDecayTrajectoryState::TsosContainer TsosContainer;
0032   typedef TwoBodyDecayTrajectory::ConstRecHitCollection ConstRecHitCollection;
0033
0034   TwoBodyDecayTrajectoryFactory(const edm::ParameterSet &config, edm::ConsumesCollector &iC);
0035   ~TwoBodyDecayTrajectoryFactory() override;
0036   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_MagFieldToken;
0037   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_globTackingToken;
0038
0039   /// Produce the trajectories.
0040   const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
0041                                                    const ConstTrajTrackPairCollection &tracks,
0042                                                    const reco::BeamSpot &beamSpot) const override;
0043
0044   const ReferenceTrajectoryCollection trajectories(const edm::EventSetup &setup,
0045                                                    const ConstTrajTrackPairCollection &tracks,
0046                                                    const ExternalPredictionCollection &external,
0047                                                    const reco::BeamSpot &beamSpot) const override;
0048
0049   TwoBodyDecayTrajectoryFactory *clone() const override { return new TwoBodyDecayTrajectoryFactory(*this); }
0050
0051 protected:
0052   const ReferenceTrajectoryCollection constructTrajectories(const ConstTrajTrackPairCollection &tracks,
0053                                                             const TwoBodyDecay &tbd,
0054                                                             const MagneticField *magField,
0055                                                             const reco::BeamSpot &beamSpot,
0056                                                             bool setParameterErrors) const;
0057
0058   bool match(const TrajectoryStateOnSurface &state, const TransientTrackingRecHit::ConstRecHitPointer &recHit) const;
0059
0060   TwoBodyDecayFitter theFitter;
0061
0062   double thePrimaryMass;
0063   double thePrimaryWidth;
0064   double theSecondaryMass;
0065
0066   double theNSigmaCutValue;
0067   double theChi2CutValue;
0068   bool theUseRefittedStateFlag;
0069   bool theConstructTsosWithErrorsFlag;
0070 };
0071
0072 ////////////////////////////////////////////////////////////////////////////////
0073 ////////////////////////////////////////////////////////////////////////////////
0074 ////////////////////////////////////////////////////////////////////////////////
0075
0076 TwoBodyDecayTrajectoryFactory::TwoBodyDecayTrajectoryFactory(const edm::ParameterSet &config,
0077                                                              edm::ConsumesCollector &iC)
0078     : TrajectoryFactoryBase(config, 2, iC),
0079       m_MagFieldToken(iC.esConsumes()),
0080       m_globTackingToken(iC.esConsumes()),
0081       theFitter(config) {
0082   const edm::ParameterSet ppc = config.getParameter<edm::ParameterSet>("ParticleProperties");
0083   thePrimaryMass = ppc.getParameter<double>("PrimaryMass");
0084   thePrimaryWidth = ppc.getParameter<double>("PrimaryWidth");
0085   theSecondaryMass = ppc.getParameter<double>("SecondaryMass");
0086
0087   theNSigmaCutValue = config.getParameter<double>("NSigmaCut");
0088   theChi2CutValue = config.getParameter<double>("Chi2Cut");
0089   theUseRefittedStateFlag = config.getParameter<bool>("UseRefittedState");
0090   theConstructTsosWithErrorsFlag = config.getParameter<bool>("ConstructTsosWithErrors");
0091 }
0092
0093 TwoBodyDecayTrajectoryFactory::~TwoBodyDecayTrajectoryFactory() {}
0094
0095 const TrajectoryFactoryBase::ReferenceTrajectoryCollection TwoBodyDecayTrajectoryFactory::trajectories(
0096     const edm::EventSetup &setup, const ConstTrajTrackPairCollection &tracks, const reco::BeamSpot &beamSpot) const {
0097   ReferenceTrajectoryCollection trajectories;
0098
0099   const MagneticField *magneticField = &setup.getData(m_MagFieldToken);
0100   const GlobalTrackingGeometry *trackingGeometry = &setup.getData(m_globTackingToken);
0101   if (tracks.size() == 2) {
0102     // produce transient tracks from persistent tracks
0103     std::vector<reco::TransientTrack> transientTracks(2);
0104
0105     transientTracks[0] = reco::TransientTrack(*tracks[0].second, magneticField);
0106     transientTracks[0].setTrackingGeometry(trackingGeometry);
0107
0108     transientTracks[1] = reco::TransientTrack(*tracks[1].second, magneticField);
0109     transientTracks[1].setTrackingGeometry(trackingGeometry);
0110
0111     // estimate the decay parameters
0112     VirtualMeasurement vm(thePrimaryMass, thePrimaryWidth, theSecondaryMass, beamSpot);
0113     TwoBodyDecay tbd = theFitter.estimate(transientTracks, vm);
0114
0115     if (!tbd.isValid() || (tbd.chi2() > theChi2CutValue)) {
0116       trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
0117       return trajectories;
0118     }
0119
0120     return constructTrajectories(tracks, tbd, magneticField, beamSpot, false);
0121   } else {
0122     edm::LogInfo("ReferenceTrajectories") << "@SUB=TwoBodyDecayTrajectoryFactory::trajectories"
0123                                           << "Need 2 tracks, got " << tracks.size() << ".\n";
0124   }
0125
0126   return trajectories;
0127 }
0128
0129 const TrajectoryFactoryBase::ReferenceTrajectoryCollection TwoBodyDecayTrajectoryFactory::trajectories(
0130     const edm::EventSetup &setup,
0131     const ConstTrajTrackPairCollection &tracks,
0132     const ExternalPredictionCollection &external,
0133     const reco::BeamSpot &beamSpot) const {
0134   ReferenceTrajectoryCollection trajectories;
0135
0136   const MagneticField *magneticField = &setup.getData(m_MagFieldToken);
0137
0138   const GlobalTrackingGeometry *trackingGeometry = &setup.getData(m_globTackingToken);
0139
0140   if (tracks.size() == 2 && external.size() == 2) {
0141     if (external[0].isValid() && external[1].isValid())  // Include external estimates
0142     {
0143       // produce transient tracks from persistent tracks
0144       std::vector<reco::TransientTrack> transientTracks(2);
0145
0146       transientTracks[0] = reco::TransientTrack(*tracks[0].second, magneticField);
0147       transientTracks[0].setTrackingGeometry(trackingGeometry);
0148
0149       transientTracks[1] = reco::TransientTrack(*tracks[1].second, magneticField);
0150       transientTracks[1].setTrackingGeometry(trackingGeometry);
0151
0152       // estimate the decay parameters. the transient tracks are not really associated to the
0153       // the external tsos, but this is o.k., because the only information retrieved from them
0154       // is the magnetic field.
0155       VirtualMeasurement vm(thePrimaryMass, thePrimaryWidth, theSecondaryMass, beamSpot);
0156       TwoBodyDecay tbd = theFitter.estimate(transientTracks, external, vm);
0157
0158       if (!tbd.isValid() || (tbd.chi2() > theChi2CutValue)) {
0159         trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
0160         return trajectories;
0161       }
0162
0163       return constructTrajectories(tracks, tbd, magneticField, beamSpot, true);
0164     } else {
0165       // Return without external estimate
0166       trajectories = this->trajectories(setup, tracks, beamSpot);
0167     }
0168   } else {
0169     edm::LogInfo("ReferenceTrajectories") << "@SUB=TwoBodyDecayTrajectoryFactory::trajectories"
0170                                           << "Need 2 tracks, got " << tracks.size() << ".\n";
0171   }
0172
0173   return trajectories;
0174 }
0175
0176 const TwoBodyDecayTrajectoryFactory::ReferenceTrajectoryCollection TwoBodyDecayTrajectoryFactory::constructTrajectories(
0177     const ConstTrajTrackPairCollection &tracks,
0178     const TwoBodyDecay &tbd,
0179     const MagneticField *magField,
0180     const reco::BeamSpot &beamSpot,
0181     bool setParameterErrors) const {
0182   ReferenceTrajectoryCollection trajectories;
0183
0184   // get innermost valid trajectory state and hits from the tracks
0185   TrajectoryInput input1 = this->innermostStateAndRecHits(tracks[0]);
0186   TrajectoryInput input2 = this->innermostStateAndRecHits(tracks[1]);
0187
0188   if (!(input1.first.isValid() && input2.first.isValid()))
0189     return trajectories;
0190
0191   // produce TwoBodyDecayTrajectoryState (input for TwoBodyDecayTrajectory)
0192   TsosContainer tsos(input1.first, input2.first);
0193   ConstRecHitCollection recHits(input1.second, input2.second);
0194   TwoBodyDecayTrajectoryState trajectoryState(tsos, tbd, theSecondaryMass, magField);
0195
0196   if (!trajectoryState.isValid()) {
0197     trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
0198     return trajectories;
0199   }
0200
0201   // always use the refitted trajectory state for matching
0202   // FIXME FIXME CLONE
0203   //TrackingRecHit::ConstRecHitPointer updatedRecHit1( recHits.first.front()->clone( tsos.first ) );
0204   // TrackingRecHit::ConstRecHitPointer updatedRecHit2( recHits.second.front()->clone( tsos.second ) );
0205
0206   TrackingRecHit::ConstRecHitPointer updatedRecHit1(recHits.first.front());
0207   TrackingRecHit::ConstRecHitPointer updatedRecHit2(recHits.second.front());
0208
0209   bool valid1 = match(trajectoryState.trajectoryStates(true).first, updatedRecHit1);
0210
0211   bool valid2 = match(trajectoryState.trajectoryStates(true).second, updatedRecHit2);
0212
0213   if (!valid1 || !valid2) {
0214     trajectories.push_back(ReferenceTrajectoryPtr(new TwoBodyDecayTrajectory()));
0215     return trajectories;
0216   }
0217
0218   ReferenceTrajectoryBase::Config config(materialEffects(), propagationDirection());
0219   config.useBeamSpot = useBeamSpot_;
0220   config.includeAPEs = includeAPEs_;
0221   config.allowZeroMaterial = allowZeroMaterial_;
0222   config.useRefittedState = theUseRefittedStateFlag;
0223   config.constructTsosWithErrors = theConstructTsosWithErrorsFlag;
0224   // set the flag for reversing the RecHits to false, since they are already in the correct order.
0225   config.hitsAreReverse = false;
0226   TwoBodyDecayTrajectory *result = new TwoBodyDecayTrajectory(trajectoryState, recHits, magField, beamSpot, config);
0227   if (setParameterErrors && tbd.hasError())
0228     result->setParameterErrors(tbd.covariance());
0229   trajectories.push_back(ReferenceTrajectoryPtr(result));
0230   return trajectories;
0231 }
0232
0233 bool TwoBodyDecayTrajectoryFactory::match(const TrajectoryStateOnSurface &state,
0234                                           const TransientTrackingRecHit::ConstRecHitPointer &recHit) const {
0235   LocalPoint lp1 = state.localPosition();
0236   LocalPoint lp2 = recHit->localPosition();
0237
0238   double deltaX = lp1.x() - lp2.x();
0239   double deltaY = lp1.y() - lp2.y();
0240
0241   LocalError le = recHit->localPositionError();
0242
0243   double varX = le.xx();
0244   double varY = le.yy();
0245
0246   return ((fabs(deltaX) / sqrt(varX) < theNSigmaCutValue) && (fabs(deltaY) / sqrt(varY) < theNSigmaCutValue));
0247 }
0248
0249 DEFINE_EDM_PLUGIN(TrajectoryFactoryPlugin, TwoBodyDecayTrajectoryFactory, "TwoBodyDecayTrajectoryFactory");