btau/plugins/HLTmumutktkVtxProducer.cc

0001 #include <algorithm>
0002 #include <cmath>
0003 #include <vector>
0004
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0011 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0012 #include "DataFormats/VertexReco/interface/Vertex.h"
0013 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0014 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0015 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0016 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0017 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0018 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0019 #include "HLTmumutktkVtxProducer.h"
0020
0021 using namespace edm;
0022 using namespace reco;
0023 using namespace std;
0024 using namespace trigger;
0025
0026 // ----------------------------------------------------------------------
0027 HLTmumutktkVtxProducer::HLTmumutktkVtxProducer(const edm::ParameterSet& iConfig)
0028     : transientTrackRecordToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0029       muCandTag_(iConfig.getParameter<edm::InputTag>("MuCand")),
0030       muCandToken_(consumes<reco::RecoChargedCandidateCollection>(muCandTag_)),
0031       trkCandTag_(iConfig.getParameter<edm::InputTag>("TrackCand")),
0032       trkCandToken_(consumes<reco::RecoChargedCandidateCollection>(trkCandTag_)),
0033       previousCandTag_(iConfig.getParameter<edm::InputTag>("PreviousCandTag")),
0034       previousCandToken_(consumes<trigger::TriggerFilterObjectWithRefs>(previousCandTag_)),
0035       mfName_(iConfig.getParameter<std::string>("SimpleMagneticField")),
0036       idealMagneticFieldRecordToken_(esConsumes(edm::ESInputTag("", mfName_))),
0037       thirdTrackMass_(iConfig.getParameter<double>("ThirdTrackMass")),
0038       fourthTrackMass_(iConfig.getParameter<double>("FourthTrackMass")),
0039       maxEta_(iConfig.getParameter<double>("MaxEta")),
0040       minPt_(iConfig.getParameter<double>("MinPt")),
0041       minInvMass_(iConfig.getParameter<double>("MinInvMass")),
0042       maxInvMass_(iConfig.getParameter<double>("MaxInvMass")),
0043       minTrkTrkMass_(iConfig.getParameter<double>("MinTrkTrkMass")),
0044       maxTrkTrkMass_(iConfig.getParameter<double>("MaxTrkTrkMass")),
0045       minD0Significance_(iConfig.getParameter<double>("MinD0Significance")),
0046       oppositeSign_(iConfig.getParameter<bool>("OppositeSign")),
0047       overlapDR_(iConfig.getParameter<double>("OverlapDR")),
0048       beamSpotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0049       beamSpotToken_(consumes<reco::BeamSpot>(beamSpotTag_)) {
0050   produces<VertexCollection>();
0051 }
0052
0053 // ----------------------------------------------------------------------
0054 HLTmumutktkVtxProducer::~HLTmumutktkVtxProducer() = default;
0055
0056 void HLTmumutktkVtxProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0057   edm::ParameterSetDescription desc;
0058   desc.add<edm::InputTag>("MuCand", edm::InputTag("hltMuTracks"));
0059   desc.add<edm::InputTag>("TrackCand", edm::InputTag("hltMumukAllConeTracks"));
0060   desc.add<edm::InputTag>("PreviousCandTag", edm::InputTag("hltDisplacedmumuFilterDoubleMu4Jpsi"));
0061   desc.add<std::string>("SimpleMagneticField", "");
0062   desc.add<double>("ThirdTrackMass", 0.493677);
0063   desc.add<double>("FourthTrackMass", 0.493677);
0064   desc.add<double>("MaxEta", 2.5);
0065   desc.add<double>("MinPt", 0.0);
0066   desc.add<double>("MinInvMass", 0.0);
0067   desc.add<double>("MaxInvMass", 99999.);
0068   desc.add<double>("MinTrkTrkMass", 0.0);
0069   desc.add<double>("MaxTrkTrkMass", 99999.);
0070   desc.add<double>("MinD0Significance", 0.0);
0071   desc.add<bool>("OppositeSign", false);
0072   desc.add<double>("OverlapDR", 0.001);
0073   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOfflineBeamSpot"));
0074   descriptions.add("HLTmumutktkVtxProducer", desc);
0075 }
0076
0077 // ----------------------------------------------------------------------
0078 void HLTmumutktkVtxProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0079   const double MuMass(0.106);
0080   const double MuMass2(MuMass * MuMass);
0081   const double thirdTrackMass2(thirdTrackMass_ * thirdTrackMass_);
0082   const double fourthTrackMass2(fourthTrackMass_ * fourthTrackMass_);
0083
0084   // get hold of muon trks
0085   Handle<RecoChargedCandidateCollection> mucands;
0086   iEvent.getByToken(muCandToken_, mucands);
0087
0088   // get the transient track builder
0089   auto const& theB = iSetup.getHandle(transientTrackRecordToken_);
0090
0091   // get the beamspot position
0092   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0093   iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
0094
0095   // get the b field
0096   auto const& bFieldHandle = iSetup.getHandle(idealMagneticFieldRecordToken_);
0097   const MagneticField* magField = bFieldHandle.product();
0098   TSCBLBuilderNoMaterial blsBuilder;
0099
0100   // get track candidates around displaced muons
0101   Handle<RecoChargedCandidateCollection> trkcands;
0102   iEvent.getByToken(trkCandToken_, trkcands);
0103
0104   unique_ptr<VertexCollection> vertexCollection(new VertexCollection());
0105
0106   // Ref to Candidate object to be recorded in filter object
0107   RecoChargedCandidateRef refMu1;
0108   RecoChargedCandidateRef refMu2;
0109   RecoChargedCandidateRef refTrk1;
0110   RecoChargedCandidateRef refTrk2;
0111
0112   double e1, e2, e3_m3, e3_m4, e4_m3, e4_m4;
0113   Particle::LorentzVector p, pBar, p1, p2, p3_m3, p3_m4, p4_m3, p4_m4, p_m3m4, p_m4m3;
0114
0115   if (mucands->size() < 2)
0116     return;
0117   if (trkcands->size() < 2)
0118     return;
0119
0120   RecoChargedCandidateCollection::const_iterator mucand1;
0121   RecoChargedCandidateCollection::const_iterator mucand2;
0122   RecoChargedCandidateCollection::const_iterator trkcand1;
0123   RecoChargedCandidateCollection::const_iterator trkcand2;
0124
0125   // get the objects passing the previous filter
0126   Handle<TriggerFilterObjectWithRefs> previousCands;
0127   iEvent.getByToken(previousCandToken_, previousCands);
0128
0129   vector<RecoChargedCandidateRef> vPrevCands;
0130   previousCands->getObjects(TriggerMuon, vPrevCands);
0131
0132   for (mucand1 = mucands->begin(); mucand1 != mucands->end(); ++mucand1) {
0133     TrackRef trk1 = mucand1->get<TrackRef>();
0134     LogDebug("HLTmumutktkVtxProducer") << " 1st muon: q*pt= " << trk1->charge() * trk1->pt() << ", eta= " << trk1->eta()
0135                                        << ", hits= " << trk1->numberOfValidHits();
0136
0137     //first check if this muon passed the previous filter
0138     if (!checkPreviousCand(trk1, vPrevCands))
0139       continue;
0140     // eta and pt cut
0141     if (fabs(trk1->eta()) > maxEta_)
0142       continue;
0143     if (trk1->pt() < minPt_)
0144       continue;
0145
0146     mucand2 = mucand1;
0147     ++mucand2;
0148     for (; mucand2 != mucands->end(); mucand2++) {
0149       TrackRef trk2 = mucand2->get<TrackRef>();
0150       if (overlap(trk1, trk2))
0151         continue;
0152
0153       LogDebug("HLTDisplacedMumukFilter") << " 2nd muon: q*pt= " << trk2->charge() * trk2->pt()
0154                                           << ", eta= " << trk2->eta() << ", hits= " << trk2->numberOfValidHits();
0155
0156       //first check if this muon passed the previous filter
0157       if (!checkPreviousCand(trk2, vPrevCands))
0158         continue;
0159       // eta and pt cut
0160       if (fabs(trk2->eta()) > maxEta_)
0161         continue;
0162       if (trk2->pt() < minPt_)
0163         continue;
0164
0165       //loop on track collection - trk1
0166       for (trkcand1 = trkcands->begin(); trkcand1 != trkcands->end(); ++trkcand1) {
0167         TrackRef trk3 = trkcand1->get<TrackRef>();
0168
0169         if (overlap(trk1, trk3))
0170           continue;
0171         if (overlap(trk2, trk3))
0172           continue;
0173
0174         LogDebug("HLTDisplacedMumukFilter") << " 3rd track: q*pt= " << trk3->charge() * trk3->pt()
0175                                             << ", eta= " << trk3->eta() << ", hits= " << trk3->numberOfValidHits();
0176
0177         // eta and pt cut
0178         if (fabs(trk3->eta()) > maxEta_)
0179           continue;
0180         if (trk3->pt() < minPt_)
0181           continue;
0182
0183         FreeTrajectoryState InitialFTS_Trk3 = initialFreeState(*trk3, magField);
0184         TrajectoryStateClosestToBeamLine tscb_Trk3(blsBuilder(InitialFTS_Trk3, *recoBeamSpotHandle));
0185         double d0sigTrk3 = tscb_Trk3.transverseImpactParameter().significance();
0186         if (d0sigTrk3 < minD0Significance_)
0187           continue;
0188
0189         //loop on track collection - trk2
0190         for (trkcand2 = trkcands->begin(); trkcand2 != trkcands->end(); ++trkcand2) {
0191           TrackRef trk4 = trkcand2->get<TrackRef>();
0192
0193           if (oppositeSign_) {
0194             if (trk3->charge() * trk4->charge() != -1)
0195               continue;
0196           }
0197           if (overlap(trk1, trk4))
0198             continue;
0199           if (overlap(trk2, trk4))
0200             continue;
0201           if (overlap(trk3, trk4))
0202             continue;
0203
0204           LogDebug("HLTDisplacedMumukFilter") << " 4th track: q*pt= " << trk4->charge() * trk4->pt()
0205                                               << ", eta= " << trk4->eta() << ", hits= " << trk4->numberOfValidHits();
0206
0207           // eta and pt cut
0208           if (fabs(trk4->eta()) > maxEta_)
0209             continue;
0210           if (trk4->pt() < minPt_)
0211             continue;
0212
0213           FreeTrajectoryState InitialFTS_Trk4 = initialFreeState(*trk4, magField);
0214           TrajectoryStateClosestToBeamLine tscb_Trk4(blsBuilder(InitialFTS_Trk4, *recoBeamSpotHandle));
0215           double d0sigTrk4 = tscb_Trk4.transverseImpactParameter().significance();
0216           if (d0sigTrk4 < minD0Significance_)
0217             continue;
0218
0219           // Combined system
0220           e1 = sqrt(trk1->momentum().Mag2() + MuMass2);
0221           e2 = sqrt(trk2->momentum().Mag2() + MuMass2);
0222           e3_m3 = sqrt(trk3->momentum().Mag2() + thirdTrackMass2);
0223           e3_m4 = sqrt(trk3->momentum().Mag2() + fourthTrackMass2);
0224           e4_m3 = sqrt(trk4->momentum().Mag2() + thirdTrackMass2);
0225           e4_m4 = sqrt(trk4->momentum().Mag2() + fourthTrackMass2);
0226
0227           p1 = Particle::LorentzVector(trk1->px(), trk1->py(), trk1->pz(), e1);
0228           p2 = Particle::LorentzVector(trk2->px(), trk2->py(), trk2->pz(), e2);
0229           p3_m3 = Particle::LorentzVector(trk3->px(), trk3->py(), trk3->pz(), e3_m3);
0230           p3_m4 = Particle::LorentzVector(trk3->px(), trk3->py(), trk3->pz(), e3_m4);
0231           p4_m3 = Particle::LorentzVector(trk4->px(), trk4->py(), trk4->pz(), e4_m3);
0232           p4_m4 = Particle::LorentzVector(trk4->px(), trk4->py(), trk4->pz(), e4_m4);
0233
0234           p = p1 + p2 + p3_m3 + p4_m4;
0235           pBar = p1 + p2 + p3_m4 + p4_m3;
0236           p_m3m4 = p3_m3 + p4_m4;
0237           p_m4m3 = p3_m4 + p4_m3;
0238
0239           //invariant mass cut
0240           if (!((p_m3m4.mass() > minTrkTrkMass_ && p_m3m4.mass() < maxTrkTrkMass_) ||
0241                 (p_m4m3.mass() > minTrkTrkMass_ && p_m4m3.mass() < maxTrkTrkMass_)))
0242             continue;
0243           if (!((p.mass() > minInvMass_ && p.mass() < maxInvMass_) ||
0244                 (pBar.mass() > minInvMass_ && pBar.mass() < maxInvMass_)))
0245             continue;
0246
0247           // do the vertex fit
0248           vector<TransientTrack> t_tks;
0249           t_tks.push_back((*theB).build(&trk1));
0250           t_tks.push_back((*theB).build(&trk2));
0251           t_tks.push_back((*theB).build(&trk3));
0252           t_tks.push_back((*theB).build(&trk4));
0253           if (t_tks.size() != 4)
0254             continue;
0255
0256           KalmanVertexFitter kvf;
0257           TransientVertex tv = kvf.vertex(t_tks);
0258           if (!tv.isValid())
0259             continue;
0260           Vertex vertex = tv;
0261
0262           vertexCollection->push_back(vertex);
0263         }
0264       }
0265     }
0266   }
0267   iEvent.put(std::move(vertexCollection));
0268 }
0269
0270 FreeTrajectoryState HLTmumutktkVtxProducer::initialFreeState(const reco::Track& tk, const MagneticField* field) {
0271   Basic3DVector<float> pos(tk.vertex());
0272   GlobalPoint gpos(pos);
0273   Basic3DVector<float> mom(tk.momentum());
0274   GlobalVector gmom(mom);
0275   GlobalTrajectoryParameters par(gpos, gmom, tk.charge(), field);
0276   CurvilinearTrajectoryError err(tk.covariance());
0277   return FreeTrajectoryState(par, err);
0278 }
0279
0280 bool HLTmumutktkVtxProducer::overlap(const TrackRef& trackref1, const TrackRef& trackref2) {
0281   if (deltaR(trackref1->eta(), trackref1->phi(), trackref2->eta(), trackref2->phi()) < overlapDR_)
0282     return true;
0283   return false;
0284 }
0285
0286 bool HLTmumutktkVtxProducer::checkPreviousCand(const TrackRef& trackref,
0287                                                const vector<RecoChargedCandidateRef>& refVect) const {
0288   bool ok = false;
0289   for (auto& i : refVect) {
0290     if (i->get<TrackRef>() == trackref) {
0291       ok = true;
0292       break;
0293     }
0294   }
0295   return ok;
0296 }