TrackProducer/plugins/TwoBodyDecayConstraintProducer.cc

0001 /**
0002  *  \class TwoBodyDecayConstraintProducer TwoBodyDecayConstraintProducer.cc RecoTracker/ConstraintProducerTest/src/TwoBodyDecayConstraintProducer.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/TwoBodyDecayFitter.h"
0029 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayVirtualMeasurement.h"
0030 #include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectoryState.h"
0031
0032 // // debug
0033 // #include <map>
0034 // #include "TH1F.h"
0035 // #include "TFile.h"
0036 // #include "TLorentzVector.h"
0037 // #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
0038
0039 class TwoBodyDecayConstraintProducer : public edm::global::EDProducer<> {
0040 public:
0041   explicit TwoBodyDecayConstraintProducer(const edm::ParameterSet&);
0042   ~TwoBodyDecayConstraintProducer() override = default;
0043
0044 private:
0045   void produce(edm::StreamID streamid, edm::Event&, const edm::EventSetup&) const override;
0046
0047   std::pair<bool, TrajectoryStateOnSurface> innermostState(const reco::TransientTrack& ttrack) const;
0048   bool match(const TrajectoryStateOnSurface& newTsos, const TrajectoryStateOnSurface& oldTsos) const;
0049
0050   const edm::InputTag srcTag_;
0051   const edm::InputTag bsSrcTag_;
0052
0053   const TwoBodyDecayFitter tbdFitter_;
0054
0055   const double primaryMass_;
0056   const double primaryWidth_;
0057   const double secondaryMass_;
0058
0059   const double sigmaPositionCutValue_;
0060   const double chi2CutValue_;
0061   const double errorRescaleValue_;
0062
0063   edm::EDGetTokenT<reco::TrackCollection> trackCollToken_;
0064   edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0065
0066   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0067   const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> trackingGeometryToken_;
0068   //   // debug
0069   //   std::map<std::string, TH1F*> histos_;
0070 };
0071
0072 TwoBodyDecayConstraintProducer::TwoBodyDecayConstraintProducer(const edm::ParameterSet& iConfig)
0073     : srcTag_(iConfig.getParameter<edm::InputTag>("src")),
0074       bsSrcTag_(iConfig.getParameter<edm::InputTag>("beamSpot")),
0075       tbdFitter_(iConfig),
0076       primaryMass_(iConfig.getParameter<double>("primaryMass")),
0077       primaryWidth_(iConfig.getParameter<double>("primaryWidth")),
0078       secondaryMass_(iConfig.getParameter<double>("secondaryMass")),
0079       sigmaPositionCutValue_(iConfig.getParameter<double>("sigmaPositionCut")),
0080       chi2CutValue_(iConfig.getParameter<double>("chi2Cut")),
0081       errorRescaleValue_(iConfig.getParameter<double>("rescaleError")),
0082       magFieldToken_(esConsumes()),
0083       trackingGeometryToken_(esConsumes()) {
0084   trackCollToken_ = consumes<reco::TrackCollection>(edm::InputTag(srcTag_));
0085   bsToken_ = consumes<reco::BeamSpot>(edm::InputTag(bsSrcTag_));
0086
0087   produces<std::vector<TrackParamConstraint> >();
0088   produces<TrackParamConstraintAssociationCollection>();
0089
0090   //   //debug
0091   //   histos_["deltaEta1"] = new TH1F( "deltaEta1", "deltaEta1", 200, -1., 1. );
0092   //   histos_["deltaP1"] = new TH1F( "deltaP1", "deltaP1", 200, -50., 50. );
0093
0094   //   histos_["deltaEta2"] = new TH1F( "deltaEta2", "deltaEta2", 200, -1., 1. );
0095   //   histos_["deltaP2"] = new TH1F( "deltaP2", "deltaP2", 200, -50., 50. );
0096 }
0097
0098 void TwoBodyDecayConstraintProducer::produce(edm::StreamID streamid,
0099                                              edm::Event& iEvent,
0100                                              const edm::EventSetup& iSetup) const {
0101   using namespace edm;
0102
0103   Handle<reco::TrackCollection> trackColl;
0104   iEvent.getByToken(trackCollToken_, trackColl);
0105
0106   Handle<reco::BeamSpot> beamSpot;
0107   iEvent.getByToken(bsToken_, beamSpot);
0108
0109   auto const* magField = &iSetup.getData(magFieldToken_);
0110
0111   auto trackingGeometry = iSetup.getHandle(trackingGeometryToken_);
0112
0113   edm::RefProd<std::vector<TrackParamConstraint> > rPairs =
0114       iEvent.getRefBeforePut<std::vector<TrackParamConstraint> >();
0115   std::unique_ptr<std::vector<TrackParamConstraint> > pairs(new std::vector<TrackParamConstraint>);
0116   std::unique_ptr<TrackParamConstraintAssociationCollection> output(
0117       new TrackParamConstraintAssociationCollection(trackColl, rPairs));
0118
0119   if (trackColl->size() == 2) {
0120     /// Construct virtual measurement (for TBD)
0121     TwoBodyDecayVirtualMeasurement vm(primaryMass_, primaryWidth_, secondaryMass_, *beamSpot.product());
0122
0123     /// Get transient tracks from track collection
0124     std::vector<reco::TransientTrack> ttracks(2);
0125     ttracks[0] = reco::TransientTrack(reco::TrackRef(trackColl, 0), magField);
0126     ttracks[0].setTrackingGeometry(trackingGeometry);
0127     ttracks[1] = reco::TransientTrack(reco::TrackRef(trackColl, 1), magField);
0128     ttracks[1].setTrackingGeometry(trackingGeometry);
0129
0130     /// Fit the TBD
0131     TwoBodyDecay tbd = tbdFitter_.estimate(ttracks, vm);
0132
0133     if (!tbd.isValid() or (tbd.chi2() > chi2CutValue_))
0134       return;
0135
0136     /// Get the innermost trajectory states
0137     std::pair<bool, TrajectoryStateOnSurface> oldInnermostState1 = innermostState(ttracks[0]);
0138     std::pair<bool, TrajectoryStateOnSurface> oldInnermostState2 = innermostState(ttracks[1]);
0139     if (!oldInnermostState1.second.isValid() || !oldInnermostState2.second.isValid())
0140       return;
0141
0142     /// Construct the TBD trajectory states
0143     TwoBodyDecayTrajectoryState::TsosContainer trackTsos(oldInnermostState1.second, oldInnermostState2.second);
0144     TwoBodyDecayTrajectoryState tbdTrajState(trackTsos, tbd, secondaryMass_, magField, true);
0145     if (!tbdTrajState.isValid())
0146       return;
0147
0148     /// Match the old and the new estimates for the trajectory state
0149     bool match1 = match(tbdTrajState.trajectoryStates(true).first, oldInnermostState1.second);
0150     bool match2 = match(tbdTrajState.trajectoryStates(true).second, oldInnermostState2.second);
0151     if (!match1 || !match2)
0152       return;
0153
0154     // re-scale error of constraintTsos
0155     tbdTrajState.rescaleError(errorRescaleValue_);
0156
0157     // produce constraint for first track
0158     pairs->push_back(tbdTrajState.trajectoryStates(true).first);
0159     output->insert(reco::TrackRef(trackColl, 0), edm::Ref<std::vector<TrackParamConstraint> >(rPairs, 0));
0160
0161     // produce constraint for second track
0162     pairs->push_back(tbdTrajState.trajectoryStates(true).second);
0163     output->insert(reco::TrackRef(trackColl, 1), edm::Ref<std::vector<TrackParamConstraint> >(rPairs, 1));
0164
0165     //     // debug
0166     //     if ( tbd.isValid() ) {
0167     //       TwoBodyDecayModel model;
0168     //       const std::pair< AlgebraicVector, AlgebraicVector > fitMomenta = model.cartesianSecondaryMomenta( tbd );
0169
0170     //       TLorentzVector recoMomentum1( ttracks[0].track().px(), ttracks[0].track().py(), ttracks[0].track().pz(),
0171     //                  sqrt((ttracks[0].track().p()*ttracks[0].track().p())+0.105658*0.105658) );
0172     //       TLorentzVector fitMomentum1( fitMomenta.first[0], fitMomenta.first[1], fitMomenta.first[2],
0173     //                 sqrt( fitMomenta.first.normsq()+0.105658*0.105658) );
0174     //       histos_["deltaP1"]->Fill( recoMomentum1.P() - fitMomentum1.P() );
0175     //       histos_["deltaEta1"]->Fill( recoMomentum1.Eta() - fitMomentum1.Eta() );
0176
0177     //       TLorentzVector recoMomentum2( ttracks[1].track().px(), ttracks[1].track().py(), ttracks[1].track().pz(),
0178     //                  sqrt((ttracks[1].track().p()*ttracks[1].track().p())+0.105658*0.105658) );
0179     //       TLorentzVector fitMomentum2( fitMomenta.second[0], fitMomenta.second[1], fitMomenta.second[2],
0180     //                 sqrt( fitMomenta.second.normsq()+0.105658*0.105658) );
0181     //       histos_["deltaP2"]->Fill( recoMomentum2.P() - fitMomentum2.P() );
0182     //       histos_["deltaEta2"]->Fill( recoMomentum2.Eta() - fitMomentum2.Eta() );
0183     //     }
0184   }
0185
0186   iEvent.put(std::move(pairs));
0187   iEvent.put(std::move(output));
0188 }
0189
0190 std::pair<bool, TrajectoryStateOnSurface> TwoBodyDecayConstraintProducer::innermostState(
0191     const reco::TransientTrack& ttrack) const {
0192   double outerR = ttrack.outermostMeasurementState().globalPosition().perp();
0193   double innerR = ttrack.innermostMeasurementState().globalPosition().perp();
0194   bool insideOut = (outerR > innerR);
0195   return std::make_pair(insideOut, insideOut ? ttrack.innermostMeasurementState() : ttrack.outermostMeasurementState());
0196 }
0197
0198 bool TwoBodyDecayConstraintProducer::match(const TrajectoryStateOnSurface& newTsos,
0199                                            const TrajectoryStateOnSurface& oldTsos) const {
0200   LocalPoint lp1 = newTsos.localPosition();
0201   LocalPoint lp2 = oldTsos.localPosition();
0202
0203   double deltaX = lp1.x() - lp2.x();
0204   double deltaY = lp1.y() - lp2.y();
0205
0206   return ((fabs(deltaX) < sigmaPositionCutValue_) && (fabs(deltaY) < sigmaPositionCutValue_));
0207 }
0208
0209 //define this as a plug-in
0210 DEFINE_FWK_MODULE(TwoBodyDecayConstraintProducer);