btau/plugins/HLTmmkkFilter.cc

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