File indexing completed on 2021-03-25 23:59:49
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
0085 Handle<RecoChargedCandidateCollection> mucands;
0086 iEvent.getByToken(muCandToken_, mucands);
0087
0088
0089 auto const& theB = iSetup.getHandle(transientTrackRecordToken_);
0090
0091
0092 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0093 iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
0094
0095
0096 auto const& bFieldHandle = iSetup.getHandle(idealMagneticFieldRecordToken_);
0097 const MagneticField* magField = bFieldHandle.product();
0098 TSCBLBuilderNoMaterial blsBuilder;
0099
0100
0101 Handle<RecoChargedCandidateCollection> trkcands;
0102 iEvent.getByToken(trkCandToken_, trkcands);
0103
0104 unique_ptr<VertexCollection> vertexCollection(new VertexCollection());
0105
0106
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
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
0138 if (!checkPreviousCand(trk1, vPrevCands))
0139 continue;
0140
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
0157 if (!checkPreviousCand(trk2, vPrevCands))
0158 continue;
0159
0160 if (fabs(trk2->eta()) > maxEta_)
0161 continue;
0162 if (trk2->pt() < minPt_)
0163 continue;
0164
0165
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
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
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
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
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
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
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 }