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