File indexing completed on 2024-04-06 12:15:50
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
0009 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0010 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0011
0012 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0013 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0014
0015 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0016 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0017 #include "DataFormats/VertexReco/interface/Vertex.h"
0018 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0019 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0020 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0021 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0022 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0023
0024 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0025
0026 #include "DataFormats/Math/interface/deltaPhi.h"
0027
0028 #include "HLTmmkFilter.h"
0029
0030 using namespace edm;
0031 using namespace reco;
0032 using namespace std;
0033 using namespace trigger;
0034
0035
0036 HLTmmkFilter::HLTmmkFilter(const edm::ParameterSet& iConfig)
0037 : HLTFilter(iConfig),
0038 transientTrackRecordToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0039 idealMagneticFieldRecordToken_(esConsumes()),
0040 muCandTag_(iConfig.getParameter<edm::InputTag>("MuCand")),
0041 muCandToken_(consumes<reco::RecoChargedCandidateCollection>(muCandTag_)),
0042 trkCandTag_(iConfig.getParameter<edm::InputTag>("TrackCand")),
0043 trkCandToken_(consumes<reco::RecoChargedCandidateCollection>(trkCandTag_)),
0044 thirdTrackMass_(iConfig.getParameter<double>("ThirdTrackMass")),
0045 maxEta_(iConfig.getParameter<double>("MaxEta")),
0046 minPt_(iConfig.getParameter<double>("MinPt")),
0047 minInvMass_(iConfig.getParameter<double>("MinInvMass")),
0048 maxInvMass_(iConfig.getParameter<double>("MaxInvMass")),
0049 maxNormalisedChi2_(iConfig.getParameter<double>("MaxNormalisedChi2")),
0050 minLxySignificance_(iConfig.getParameter<double>("MinLxySignificance")),
0051 minCosinePointingAngle_(iConfig.getParameter<double>("MinCosinePointingAngle")),
0052 minD0Significance_(iConfig.getParameter<double>("MinD0Significance")),
0053 fastAccept_(iConfig.getParameter<bool>("FastAccept")),
0054 beamSpotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0055 beamSpotToken_(consumes<reco::BeamSpot>(beamSpotTag_)) {
0056 produces<VertexCollection>();
0057 produces<CandidateCollection>();
0058 }
0059
0060
0061 HLTmmkFilter::~HLTmmkFilter() = default;
0062
0063 void HLTmmkFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0064 edm::ParameterSetDescription desc;
0065 makeHLTFilterDescription(desc);
0066 desc.add<edm::InputTag>("MuCand", edm::InputTag("hltMuTracks"));
0067 desc.add<edm::InputTag>("TrackCand", edm::InputTag("hltMumukAllConeTracks"));
0068 desc.add<double>("ThirdTrackMass", 0.106);
0069 desc.add<double>("MaxEta", 2.5);
0070 desc.add<double>("MinPt", 3.0);
0071 desc.add<double>("MinInvMass", 1.2);
0072 desc.add<double>("MaxInvMass", 2.2);
0073 desc.add<double>("MaxNormalisedChi2", 10.0);
0074 desc.add<double>("MinLxySignificance", 3.0);
0075 desc.add<double>("MinCosinePointingAngle", 0.9);
0076 desc.add<double>("MinD0Significance", 0.0);
0077 desc.add<bool>("FastAccept", false);
0078 desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOfflineBeamSpot"));
0079 descriptions.add("hltmmkFilter", desc);
0080 }
0081
0082
0083 void HLTmmkFilter::beginJob() {}
0084
0085
0086 void HLTmmkFilter::endJob() {}
0087
0088
0089 bool HLTmmkFilter::hltFilter(edm::Event& iEvent,
0090 const edm::EventSetup& iSetup,
0091 trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0092 const double MuMass(0.106);
0093 const double MuMass2(MuMass * MuMass);
0094
0095 const double thirdTrackMass2(thirdTrackMass_ * thirdTrackMass_);
0096
0097 unique_ptr<CandidateCollection> output(new CandidateCollection());
0098 unique_ptr<VertexCollection> vertexCollection(new VertexCollection());
0099
0100
0101 auto const& theB = iSetup.getHandle(transientTrackRecordToken_);
0102
0103
0104 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0105 iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
0106 const reco::BeamSpot& vertexBeamSpot = *recoBeamSpotHandle;
0107
0108 auto const& bFieldHandle = iSetup.getHandle(idealMagneticFieldRecordToken_);
0109 const MagneticField* magField = bFieldHandle.product();
0110
0111 TSCBLBuilderNoMaterial blsBuilder;
0112
0113
0114 RecoChargedCandidateRef refMu1;
0115 RecoChargedCandidateRef refMu2;
0116 RecoChargedCandidateRef refTrk;
0117
0118
0119 Handle<RecoChargedCandidateCollection> mucands;
0120 iEvent.getByToken(muCandToken_, mucands);
0121
0122
0123 Handle<RecoChargedCandidateCollection> trkcands;
0124 iEvent.getByToken(trkCandToken_, trkcands);
0125
0126 if (saveTags()) {
0127 filterproduct.addCollectionTag(muCandTag_);
0128 filterproduct.addCollectionTag(trkCandTag_);
0129 }
0130
0131 double e1, e2, e3;
0132 Particle::LorentzVector p, p1, p2, p3;
0133
0134
0135 vector<TrackRef> trkMuCands;
0136
0137
0138 vector<bool> isUsedCand(trkcands->size(), false);
0139
0140 int counter = 0;
0141
0142
0143 for (auto mucand1 = mucands->begin(), endCand1 = mucands->end(); mucand1 != endCand1; ++mucand1) {
0144 if (mucands->size() < 2)
0145 break;
0146 if (trkcands->empty())
0147 break;
0148
0149 TrackRef trk1 = mucand1->get<TrackRef>();
0150 LogDebug("HLTDisplacedMumukFilter") << " 1st muon: q*pt= " << trk1->charge() * trk1->pt()
0151 << ", eta= " << trk1->eta() << ", hits= " << trk1->numberOfValidHits();
0152
0153
0154 if (fabs(trk1->eta()) > maxEta_)
0155 continue;
0156
0157
0158 if (trk1->pt() < minPt_)
0159 continue;
0160
0161 auto mucand2 = mucand1;
0162 ++mucand2;
0163
0164 for (auto endCand2 = mucands->end(); mucand2 != endCand2; ++mucand2) {
0165 TrackRef trk2 = mucand2->get<TrackRef>();
0166
0167 LogDebug("HLTDisplacedMumukFilter") << " 2nd muon: q*pt= " << trk2->charge() * trk2->pt()
0168 << ", eta= " << trk2->eta() << ", hits= " << trk2->numberOfValidHits();
0169
0170
0171 if (fabs(trk2->eta()) > maxEta_)
0172 continue;
0173
0174
0175 if (trk2->pt() < minPt_)
0176 continue;
0177
0178 RecoChargedCandidateCollection::const_iterator trkcand, endCandTrk;
0179
0180 std::vector<bool>::iterator isUsedIter, endIsUsedCand;
0181
0182
0183 for (trkcand = trkcands->begin(),
0184 endCandTrk = trkcands->end(),
0185 isUsedIter = isUsedCand.begin(),
0186 endIsUsedCand = isUsedCand.end();
0187 trkcand != endCandTrk && isUsedIter != endIsUsedCand;
0188 ++trkcand, ++isUsedIter) {
0189 TrackRef trk3 = trkcand->get<TrackRef>();
0190
0191
0192 if (overlap(*mucand1, *trkcand)) {
0193 trkMuCands.push_back(trk3);
0194 *isUsedIter = true;
0195 continue;
0196 } else if (overlap(*mucand2, *trkcand)) {
0197 trkMuCands.push_back(trk3);
0198 *isUsedIter = true;
0199 continue;
0200 }
0201
0202 if (trkMuCands.size() == 2)
0203 break;
0204 }
0205
0206
0207
0208
0209
0210 for (trkcand = trkcands->begin(),
0211 endCandTrk = trkcands->end(),
0212 isUsedIter = isUsedCand.begin(),
0213 endIsUsedCand = isUsedCand.end();
0214 trkcand != endCandTrk && isUsedIter != endIsUsedCand;
0215 ++trkcand, ++isUsedIter) {
0216 TrackRef trk3 = trkcand->get<TrackRef>();
0217
0218 LogDebug("HLTDisplacedMumukFilter") << " 3rd track: q*pt= " << trk3->charge() * trk3->pt()
0219 << ", eta= " << trk3->eta() << ", hits= " << trk3->numberOfValidHits();
0220
0221
0222 bool skip = false;
0223 for (auto& trkMuCand : trkMuCands)
0224 if (trk3 == trkMuCand)
0225 skip = true;
0226 if (skip)
0227 continue;
0228
0229
0230 if (*isUsedIter)
0231 continue;
0232
0233
0234 if (fabs(trk3->eta()) > maxEta_)
0235 continue;
0236
0237
0238 if (trk3->pt() < minPt_)
0239 continue;
0240
0241
0242 e1 = sqrt(trk1->momentum().Mag2() + MuMass2);
0243 e2 = sqrt(trk2->momentum().Mag2() + MuMass2);
0244 e3 = sqrt(trk3->momentum().Mag2() + thirdTrackMass2);
0245
0246 p1 = Particle::LorentzVector(trk1->px(), trk1->py(), trk1->pz(), e1);
0247 p2 = Particle::LorentzVector(trk2->px(), trk2->py(), trk2->pz(), e2);
0248 p3 = Particle::LorentzVector(trk3->px(), trk3->py(), trk3->pz(), e3);
0249
0250 p = p1 + p2 + p3;
0251
0252
0253 double invmass = abs(p.mass());
0254
0255 LogDebug("HLTDisplacedMumukFilter") << " Invmass= " << invmass;
0256
0257 if (invmass > maxInvMass_ || invmass < minInvMass_)
0258 continue;
0259
0260
0261 vector<TransientTrack> t_tks;
0262 t_tks.push_back((*theB).build(&trk1));
0263 t_tks.push_back((*theB).build(&trk2));
0264 t_tks.push_back((*theB).build(&trk3));
0265
0266 if (t_tks.size() != 3)
0267 continue;
0268
0269 FreeTrajectoryState InitialFTS = initialFreeState(*trk3, magField);
0270 TrajectoryStateClosestToBeamLine tscb(blsBuilder(InitialFTS, *recoBeamSpotHandle));
0271 double d0sig = tscb.transverseImpactParameter().significance();
0272
0273 if (d0sig < minD0Significance_)
0274 continue;
0275
0276 KalmanVertexFitter kvf;
0277 TransientVertex tv = kvf.vertex(t_tks);
0278
0279 if (!tv.isValid())
0280 continue;
0281
0282 Vertex vertex = tv;
0283
0284
0285 GlobalPoint secondaryVertex = tv.position();
0286 GlobalError err = tv.positionError();
0287
0288
0289 GlobalPoint displacementFromBeamspot(-1 * ((vertexBeamSpot.x0() - secondaryVertex.x()) +
0290 (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
0291 -1 * ((vertexBeamSpot.y0() - secondaryVertex.y()) +
0292 (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
0293 0);
0294
0295 float lxy = displacementFromBeamspot.perp();
0296 float lxyerr = sqrt(err.rerr(displacementFromBeamspot));
0297
0298
0299 float normChi2 = tv.normalisedChiSquared();
0300
0301
0302 Vertex::Point vperp(displacementFromBeamspot.x(), displacementFromBeamspot.y(), 0.);
0303 math::XYZVector pperp(p.x(), p.y(), 0.);
0304
0305 float cosAlpha = vperp.Dot(pperp) / (vperp.R() * pperp.R());
0306
0307 LogDebug("HLTDisplacedMumukFilter")
0308 << " vertex fit normalised chi2: " << normChi2 << ", Lxy significance: " << lxy / lxyerr
0309 << ", cosine pointing angle: " << cosAlpha;
0310
0311
0312 vertexCollection->push_back(vertex);
0313
0314 if (normChi2 > maxNormalisedChi2_)
0315 continue;
0316 if (lxy / lxyerr < minLxySignificance_)
0317 continue;
0318 if (cosAlpha < minCosinePointingAngle_)
0319 continue;
0320
0321 LogDebug("HLTDisplacedMumukFilter") << " Event passed!";
0322
0323
0324 ++counter;
0325
0326
0327 bool i1done = false;
0328 bool i2done = false;
0329 bool i3done = false;
0330 vector<RecoChargedCandidateRef> vref;
0331 filterproduct.getObjects(TriggerMuon, vref);
0332 for (auto& i : vref) {
0333 RecoChargedCandidateRef candref = RecoChargedCandidateRef(i);
0334 TrackRef trktmp = candref->get<TrackRef>();
0335 if (trktmp == trk1) {
0336 i1done = true;
0337 } else if (trktmp == trk2) {
0338 i2done = true;
0339 } else if (trktmp == trk3) {
0340 i3done = true;
0341 }
0342 if (i1done && i2done && i3done)
0343 break;
0344 }
0345
0346 if (!i1done) {
0347 refMu1 = RecoChargedCandidateRef(
0348 Ref<RecoChargedCandidateCollection>(mucands, distance(mucands->begin(), mucand1)));
0349 filterproduct.addObject(TriggerMuon, refMu1);
0350 }
0351 if (!i2done) {
0352 refMu2 = RecoChargedCandidateRef(
0353 Ref<RecoChargedCandidateCollection>(mucands, distance(mucands->begin(), mucand2)));
0354 filterproduct.addObject(TriggerMuon, refMu2);
0355 }
0356 if (!i3done) {
0357 refTrk = RecoChargedCandidateRef(
0358 Ref<RecoChargedCandidateCollection>(trkcands, distance(trkcands->begin(), trkcand)));
0359 filterproduct.addObject(TriggerTrack, refTrk);
0360 }
0361
0362 if (fastAccept_)
0363 break;
0364 }
0365
0366 trkMuCands.clear();
0367 }
0368 }
0369
0370
0371 const bool accept(counter >= 1);
0372
0373 LogDebug("HLTDisplacedMumukFilter") << " >>>>> Result of HLTDisplacedMumukFilter is " << accept
0374 << ", number of muon pairs passing thresholds= " << counter;
0375
0376 iEvent.put(std::move(vertexCollection));
0377
0378 return accept;
0379 }
0380
0381
0382 FreeTrajectoryState HLTmmkFilter::initialFreeState(const reco::Track& tk, const MagneticField* field) {
0383 Basic3DVector<float> pos(tk.vertex());
0384 GlobalPoint gpos(pos);
0385 Basic3DVector<float> mom(tk.momentum());
0386 GlobalVector gmom(mom);
0387 GlobalTrajectoryParameters par(gpos, gmom, tk.charge(), field);
0388 CurvilinearTrajectoryError err(tk.covariance());
0389 return FreeTrajectoryState(par, err);
0390 }
0391
0392 int HLTmmkFilter::overlap(const reco::Candidate& a, const reco::Candidate& b) {
0393 double eps(1.44e-4);
0394
0395 double dpt = a.pt() - b.pt();
0396 dpt *= dpt;
0397
0398 double dphi = deltaPhi(a.phi(), b.phi());
0399 dphi *= dphi;
0400
0401 double deta = a.eta() - b.eta();
0402 deta *= deta;
0403
0404 if ((dphi + deta) < eps) {
0405 return 1;
0406 }
0407
0408 return 0;
0409 }