TrackProducer/plugins/TwoBodyDecayMomConstraintProducer.cc

0001 /**
0002  *  \class TwoBodyDecayMomConstraintProducer TwoBodyDecayMomConstraintProducer.cc RecoTracker/ConstraintProducerTest/src/TwoBodyDecayMomConstraintProducer.cc
0003  *
0004  *  Description: Produces track parameter constraints for refitting tracks, according to information TwoBodyDecay kinematic fit.
0005  *
0006  *  Original Author:  Edmund Widl
0007  *
0008  */
0009
0010 #include <memory>
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/global/EDProducer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/Utilities/interface/InputTag.h"
0018
0019 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0020 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0021
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0024
0025 #include "TrackingTools/PatternTools/interface/TrackConstraintAssociation.h"
0026
0027 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecay.h"
0028 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
0029 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayFitter.h"
0030 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayVirtualMeasurement.h"
0031 #include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectoryState.h"
0032
0033 // // debug
0034 // #include <map>
0035 // #include "TH1F.h"
0036 // #include "TFile.h"
0037 // #include "TLorentzVector.h"
0038 // #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
0039
0040 class TwoBodyDecayMomConstraintProducer : public edm::global::EDProducer<> {
0041 public:
0042   explicit TwoBodyDecayMomConstraintProducer(const edm::ParameterSet&);
0043   ~TwoBodyDecayMomConstraintProducer() override = default;
0044
0045 private:
0046   void produce(edm::StreamID streamid, edm::Event&, const edm::EventSetup&) const override;
0047
0048   std::pair<double, double> momentaAtVertex(const TwoBodyDecay& tbd) const;
0049
0050   std::pair<double, double> momentaAtInnermostSurface(const TwoBodyDecay& tbd,
0051                                                       const std::vector<reco::TransientTrack>& ttracks,
0052                                                       const MagneticField* magField) const;
0053
0054   std::pair<bool, TrajectoryStateOnSurface> innermostState(const reco::TransientTrack& ttrack) const;
0055   bool match(const TrajectoryStateOnSurface& newTsos, const TrajectoryStateOnSurface& oldTsos) const;
0056
0057   const edm::InputTag srcTag_;
0058   const edm::InputTag bsSrcTag_;
0059
0060   const TwoBodyDecayFitter tbdFitter_;
0061
0062   const double primaryMass_;
0063   const double primaryWidth_;
0064   const double secondaryMass_;
0065
0066   const double sigmaPositionCutValue_;
0067   const double chi2CutValue_;
0068   const double fixedMomentumError_;
0069
0070   enum MomentumForRefitting { atVertex, atInnermostSurface };
0071   const MomentumForRefitting momentumForRefitting_;
0072
0073   static MomentumForRefitting momentumForRefittingFromString(std::string momentumForRefittingString);
0074
0075   edm::EDGetTokenT<reco::TrackCollection> trackCollToken_;
0076   edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0077
0078   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0079   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> trackingGeometryToken_;
0080   //   // debug
0081   //   std::map<std::string, TH1F*> histos_;
0082 };
0083
0084 TwoBodyDecayMomConstraintProducer::TwoBodyDecayMomConstraintProducer(const edm::ParameterSet& iConfig)
0085     : srcTag_(iConfig.getParameter<edm::InputTag>("src")),
0086       bsSrcTag_(iConfig.getParameter<edm::InputTag>("beamSpot")),
0087       tbdFitter_(iConfig),
0088       primaryMass_(iConfig.getParameter<double>("primaryMass")),
0089       primaryWidth_(iConfig.getParameter<double>("primaryWidth")),
0090       secondaryMass_(iConfig.getParameter<double>("secondaryMass")),
0091       sigmaPositionCutValue_(iConfig.getParameter<double>("sigmaPositionCut")),
0092       chi2CutValue_(iConfig.getParameter<double>("chi2Cut")),
0093       fixedMomentumError_(iConfig.getParameter<double>("fixedMomentumError")),
0094       momentumForRefitting_(momentumForRefittingFromString(iConfig.getParameter<std::string>("momentumForRefitting"))),
0095       magFieldToken_(esConsumes()),
0096       trackingGeometryToken_(esConsumes()) {
0097   trackCollToken_ = consumes<reco::TrackCollection>(edm::InputTag(srcTag_));
0098   bsToken_ = consumes<reco::BeamSpot>(edm::InputTag(bsSrcTag_));
0099
0100   produces<std::vector<MomentumConstraint> >();
0101   produces<TrackMomConstraintAssociationCollection>();
0102
0103   //   //debug
0104   //   histos_["deltaEta1"] = new TH1F( "deltaEta1", "deltaEta1", 200, -1., 1. );
0105   //   histos_["deltaP1"] = new TH1F( "deltaP1", "deltaP1", 200, -50., 50. );
0106
0107   //   histos_["deltaEta2"] = new TH1F( "deltaEta2", "deltaEta2", 200, -1., 1. );
0108   //   histos_["deltaP2"] = new TH1F( "deltaP2", "deltaP2", 200, -50., 50. );
0109 }
0110
0111 void TwoBodyDecayMomConstraintProducer::produce(edm::StreamID streamid,
0112                                                 edm::Event& iEvent,
0113                                                 const edm::EventSetup& iSetup) const {
0114   using namespace edm;
0115
0116   Handle<reco::TrackCollection> trackColl;
0117   iEvent.getByToken(trackCollToken_, trackColl);
0118
0119   Handle<reco::BeamSpot> beamSpot;
0120   iEvent.getByToken(bsToken_, beamSpot);
0121
0122   auto const* magField = &iSetup.getData(magFieldToken_);
0123
0124   auto trackingGeometry = iSetup.getHandle(trackingGeometryToken_);
0125
0126   edm::RefProd<std::vector<MomentumConstraint> > rPairs = iEvent.getRefBeforePut<std::vector<MomentumConstraint> >();
0127   std::unique_ptr<std::vector<MomentumConstraint> > pairs(new std::vector<MomentumConstraint>);
0128   std::unique_ptr<TrackMomConstraintAssociationCollection> output(
0129       new TrackMomConstraintAssociationCollection(trackColl, rPairs));
0130
0131   if (trackColl->size() == 2) {
0132     /// Construct virtual measurement (for TBD)
0133     TwoBodyDecayVirtualMeasurement vm(primaryMass_, primaryWidth_, secondaryMass_, *beamSpot.product());
0134
0135     /// Get transient tracks from track collection
0136     std::vector<reco::TransientTrack> ttracks(2);
0137     ttracks[0] = reco::TransientTrack(reco::TrackRef(trackColl, 0), magField);
0138     ttracks[0].setTrackingGeometry(trackingGeometry);
0139     ttracks[1] = reco::TransientTrack(reco::TrackRef(trackColl, 1), magField);
0140     ttracks[1].setTrackingGeometry(trackingGeometry);
0141
0142     /// Fit the TBD
0143     TwoBodyDecay tbd = tbdFitter_.estimate(ttracks, vm);
0144
0145     if (!tbd.isValid() or (tbd.chi2() > chi2CutValue_))
0146       return;
0147
0148     std::pair<double, double> fitMomenta;
0149     if (momentumForRefitting_ == atVertex) {
0150       fitMomenta = momentaAtVertex(tbd);
0151     } else if (momentumForRefitting_ == atInnermostSurface) {
0152       fitMomenta = momentaAtInnermostSurface(tbd, ttracks, magField);
0153     }  // no other possibility!!!
0154
0155     if ((fitMomenta.first > 0.) and (fitMomenta.second > 0.)) {
0156       // produce constraint for first track
0157       MomentumConstraint constraint1(fitMomenta.first, fixedMomentumError_);
0158       pairs->push_back(constraint1);
0159       output->insert(reco::TrackRef(trackColl, 0), edm::Ref<std::vector<MomentumConstraint> >(rPairs, 0));
0160
0161       // produce constraint for second track
0162       MomentumConstraint constraint2(fitMomenta.second, fixedMomentumError_);
0163       pairs->push_back(constraint2);
0164       output->insert(reco::TrackRef(trackColl, 1), edm::Ref<std::vector<MomentumConstraint> >(rPairs, 1));
0165     }
0166
0167     //     // debug
0168     //     if ( tbd.isValid() and ( fitMomenta.first > 0. ) and ( fitMomenta.second > 0. ) ) {
0169     //       TwoBodyDecayModel model;
0170     //       const std::pair< AlgebraicVector, AlgebraicVector > fitMomenta = model.cartesianSecondaryMomenta( tbd );
0171
0172     //       TLorentzVector recoMomentum1( ttracks[0].track().px(), ttracks[0].track().py(), ttracks[0].track().pz(),
0173     //                  sqrt((ttracks[0].track().p()*ttracks[0].track().p())+0.105658*0.105658) );
0174     //       TLorentzVector fitMomentum1( fitMomenta.first[0], fitMomenta.first[1], fitMomenta.first[2],
0175     //                 sqrt( fitMomenta.first.normsq()+0.105658*0.105658) );
0176     //       histos_["deltaP1"]->Fill( recoMomentum1.P() - fitMomentum1.P() );
0177     //       histos_["deltaEta1"]->Fill( recoMomentum1.Eta() - fitMomentum1.Eta() );
0178
0179     //       TLorentzVector recoMomentum2( ttracks[1].track().px(), ttracks[1].track().py(), ttracks[1].track().pz(),
0180     //                  sqrt((ttracks[1].track().p()*ttracks[1].track().p())+0.105658*0.105658) );
0181     //       TLorentzVector fitMomentum2( fitMomenta.second[0], fitMomenta.second[1], fitMomenta.second[2],
0182     //                 sqrt( fitMomenta.second.normsq()+0.105658*0.105658) );
0183     //       histos_["deltaP2"]->Fill( recoMomentum2.P() - fitMomentum2.P() );
0184     //       histos_["deltaEta2"]->Fill( recoMomentum2.Eta() - fitMomentum2.Eta() );
0185     //     }
0186   }
0187
0188   iEvent.put(std::move(pairs));
0189   iEvent.put(std::move(output));
0190 }
0191
0192 std::pair<double, double> TwoBodyDecayMomConstraintProducer::momentaAtVertex(const TwoBodyDecay& tbd) const {
0193   // construct global trajectory parameters at the starting point
0194   TwoBodyDecayModel tbdDecayModel(tbd[TwoBodyDecayParameters::mass], secondaryMass_);
0195   std::pair<AlgebraicVector, AlgebraicVector> secondaryMomenta = tbdDecayModel.cartesianSecondaryMomenta(tbd);
0196
0197   return std::make_pair(secondaryMomenta.first.norm(), secondaryMomenta.second.norm());
0198 }
0199
0200 std::pair<double, double> TwoBodyDecayMomConstraintProducer::momentaAtInnermostSurface(
0201     const TwoBodyDecay& tbd, const std::vector<reco::TransientTrack>& ttracks, const MagneticField* magField) const {
0202   std::pair<double, double> result = std::make_pair(-1., -1.);
0203
0204   /// Get the innermost trajectory states
0205   std::pair<bool, TrajectoryStateOnSurface> oldInnermostState1 = innermostState(ttracks[0]);
0206   std::pair<bool, TrajectoryStateOnSurface> oldInnermostState2 = innermostState(ttracks[1]);
0207   if (!oldInnermostState1.second.isValid() || !oldInnermostState2.second.isValid())
0208     return result;
0209
0210   /// Construct the TBD trajectory states
0211   TwoBodyDecayTrajectoryState::TsosContainer trackTsos(oldInnermostState1.second, oldInnermostState2.second);
0212   TwoBodyDecayTrajectoryState tbdTrajState(trackTsos, tbd, secondaryMass_, magField, true);
0213   if (!tbdTrajState.isValid())
0214     return result;
0215
0216   /// Match the old and the new estimates for the trajectory state
0217   bool match1 = match(tbdTrajState.trajectoryStates(true).first, oldInnermostState1.second);
0218   bool match2 = match(tbdTrajState.trajectoryStates(true).second, oldInnermostState2.second);
0219   if (!match1 || !match2)
0220     return result;
0221
0222   result = std::make_pair(fabs(1. / tbdTrajState.trajectoryStates(true).first.localParameters().qbp()),
0223                           fabs(1. / tbdTrajState.trajectoryStates(true).second.localParameters().qbp()));
0224   return result;
0225 }
0226
0227 std::pair<bool, TrajectoryStateOnSurface> TwoBodyDecayMomConstraintProducer::innermostState(
0228     const reco::TransientTrack& ttrack) const {
0229   double outerR = ttrack.outermostMeasurementState().globalPosition().perp();
0230   double innerR = ttrack.innermostMeasurementState().globalPosition().perp();
0231   bool insideOut = (outerR > innerR);
0232   return std::make_pair(insideOut, insideOut ? ttrack.innermostMeasurementState() : ttrack.outermostMeasurementState());
0233 }
0234
0235 bool TwoBodyDecayMomConstraintProducer::match(const TrajectoryStateOnSurface& newTsos,
0236                                               const TrajectoryStateOnSurface& oldTsos) const {
0237   LocalPoint lp1 = newTsos.localPosition();
0238   LocalPoint lp2 = oldTsos.localPosition();
0239
0240   double deltaX = lp1.x() - lp2.x();
0241   double deltaY = lp1.y() - lp2.y();
0242
0243   return ((fabs(deltaX) < sigmaPositionCutValue_) && (fabs(deltaY) < sigmaPositionCutValue_));
0244 }
0245
0246 TwoBodyDecayMomConstraintProducer::MomentumForRefitting
0247 TwoBodyDecayMomConstraintProducer::momentumForRefittingFromString(std::string strMomentumForRefitting) {
0248   if (strMomentumForRefitting == "atVertex") {
0249     return atVertex;
0250   } else if (strMomentumForRefitting == "atInnermostSurface") {
0251     return atInnermostSurface;
0252   } else {
0253     throw cms::Exception("TwoBodyDecayMomConstraintProducer") << "value of config  variable 'momentumForRefitting': "
0254                                                               << "has to be 'atVertex' or 'atInnermostSurface'";
0255   }
0256 }
0257
0258 //define this as a plug-in
0259 DEFINE_FWK_MODULE(TwoBodyDecayMomConstraintProducer);