btau/plugins/HLTmmkFilter.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
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   // get the transient track builder
0101   auto const& theB = iSetup.getHandle(transientTrackRecordToken_);
0102
0103   // get the beamspot position
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   // Ref to Candidate object to be recorded in filter object
0114   RecoChargedCandidateRef refMu1;
0115   RecoChargedCandidateRef refMu2;
0116   RecoChargedCandidateRef refTrk;
0117
0118   // get hold of muon trks
0119   Handle<RecoChargedCandidateCollection> mucands;
0120   iEvent.getByToken(muCandToken_, mucands);
0121
0122   // get track candidates around displaced muons
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   //TrackRefs to mu cands in trkcand
0135   vector<TrackRef> trkMuCands;
0136
0137   //Already used mu tracks to avoid double counting candidates
0138   vector<bool> isUsedCand(trkcands->size(), false);
0139
0140   int counter = 0;
0141
0142   //run algorithm
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     // eta cut
0154     if (fabs(trk1->eta()) > maxEta_)
0155       continue;
0156
0157     // Pt threshold cut
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       // eta cut
0171       if (fabs(trk2->eta()) > maxEta_)
0172         continue;
0173
0174       // Pt threshold cut
0175       if (trk2->pt() < minPt_)
0176         continue;
0177
0178       RecoChargedCandidateCollection::const_iterator trkcand, endCandTrk;
0179
0180       std::vector<bool>::iterator isUsedIter, endIsUsedCand;
0181
0182       //get overlapping muon candidates
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         //check for overlapping muon tracks and save TrackRefs
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       //Not all overlapping candidates found, skip event
0207       //if (trkMuCands.size()!=2) continue;
0208
0209       //combine muons with all tracks
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         //skip overlapping muon candidates
0222         bool skip = false;
0223         for (auto& trkMuCand : trkMuCands)
0224           if (trk3 == trkMuCand)
0225             skip = true;
0226         if (skip)
0227           continue;
0228
0229         //skip already used tracks
0230         if (*isUsedIter)
0231           continue;
0232
0233         // eta cut
0234         if (fabs(trk3->eta()) > maxEta_)
0235           continue;
0236
0237         // Pt threshold cut
0238         if (trk3->pt() < minPt_)
0239           continue;
0240
0241         // Combined system
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         //invariant mass cut
0253         double invmass = abs(p.mass());
0254
0255         LogDebug("HLTDisplacedMumukFilter") << " Invmass= " << invmass;
0256
0257         if (invmass > maxInvMass_ || invmass < minInvMass_)
0258           continue;
0259
0260         // do the vertex fit
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         // get vertex position and error to calculate the decay length significance
0285         GlobalPoint secondaryVertex = tv.position();
0286         GlobalError err = tv.positionError();
0287
0288         //calculate decay length  significance w.r.t. the beamspot
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         // get normalizes chi2
0299         float normChi2 = tv.normalisedChiSquared();
0300
0301         //calculate the angle between the decay length and the mumu momentum
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         // put vertex in the event
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         //Add event
0324         ++counter;
0325
0326         //Get refs
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   // filter decision
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 }