Muon/plugins/HLTMuonTrackMassFilter.cc

0001 #include "HLTMuonTrackMassFilter.h"
0002
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008
0009 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0010 #include "DataFormats/HLTReco/interface/TriggerRefsCollections.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/MessageLogger/interface/MessageDrop.h"
0014
0015 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0016 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0017 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021
0022 #include "DataFormats/Math/interface/deltaR.h"
0023
0024 #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
0025 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0026
0027 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0029 #include "FWCore/Utilities/interface/InputTag.h"
0030
0031 #include <vector>
0032 #include <memory>
0033 #include <sstream>
0034
0035 HLTMuonTrackMassFilter::HLTMuonTrackMassFilter(const edm::ParameterSet& iConfig)
0036     : HLTFilter(iConfig),
0037       idealMagneticFieldRecordToken_(esConsumes()),
0038       beamspotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0039       beamspotToken_(consumes<reco::BeamSpot>(beamspotTag_)),
0040       muonTag_(iConfig.getParameter<edm::InputTag>("CandTag")),
0041       muonToken_(consumes<reco::RecoChargedCandidateCollection>(muonTag_)),
0042       trackTag_(iConfig.getParameter<edm::InputTag>("TrackTag")),
0043       trackToken_(consumes<reco::RecoChargedCandidateCollection>(trackTag_)),
0044       prevCandTag_(iConfig.getParameter<edm::InputTag>("PreviousCandTag")),
0045       prevCandToken_(consumes<trigger::TriggerFilterObjectWithRefs>(prevCandTag_)),
0046       minMasses_(iConfig.getParameter<std::vector<double> >("MinMasses")),
0047       maxMasses_(iConfig.getParameter<std::vector<double> >("MaxMasses")),
0048       checkCharge_(iConfig.getParameter<bool>("checkCharge")),
0049       minTrackPt_(iConfig.getParameter<double>("MinTrackPt")),
0050       minTrackP_(iConfig.getParameter<double>("MinTrackP")),
0051       maxTrackEta_(iConfig.getParameter<double>("MaxTrackEta")),
0052       maxTrackDxy_(iConfig.getParameter<double>("MaxTrackDxy")),
0053       maxTrackDz_(iConfig.getParameter<double>("MaxTrackDz")),
0054       minTrackHits_(iConfig.getParameter<int>("MinTrackHits")),
0055       maxTrackNormChi2_(iConfig.getParameter<double>("MaxTrackNormChi2")),
0056       //   maxDzMuonTrack_(iConfig.getParameter<double>("MaxDzMuonTrack")),
0057       max_DCAMuonTrack_(iConfig.getParameter<double>("MaxDCAMuonTrack")),
0058       cutCowboys_(iConfig.getParameter<bool>("CutCowboys")) {
0059   //
0060   // verify mass windows
0061   //
0062   bool massesValid = minMasses_.size() == maxMasses_.size();
0063   if (massesValid) {
0064     for (unsigned int i = 0; i < minMasses_.size(); ++i) {
0065       if (minMasses_[i] < 0 || maxMasses_[i] < 0 || minMasses_[i] > maxMasses_[i])
0066         massesValid = false;
0067     }
0068   }
0069   if (!massesValid) {
0070     edm::LogError("HLTMuonTrackMassFilter") << "Inconsistency in definition of mass windows, "
0071                                             << "no event will pass the filter";
0072     minMasses_.clear();
0073     maxMasses_.clear();
0074   }
0075
0076   std::ostringstream stream;
0077   stream << "instantiated with parameters\n";
0078   stream << "  beamspot = " << beamspotTag_ << "\n";
0079   stream << "  muonCandidates = " << muonTag_ << "\n";
0080   stream << "  trackCandidates = " << trackTag_ << "\n";
0081   stream << "  previousCandidates = " << prevCandTag_ << "\n";
0082   stream << "  saveTags= " << saveTags() << "\n";
0083   stream << "  mass windows =";
0084   for (size_t i = 0; i < minMasses_.size(); ++i)
0085     stream << " (" << minMasses_[i] << "," << maxMasses_[i] << ")";
0086   stream << "\n";
0087   stream << "  checkCharge = " << checkCharge_ << "\n";
0088   stream << "  MinTrackPt = " << minTrackPt_ << "\n";
0089   stream << "  MinTrackP = " << minTrackP_ << "\n";
0090   stream << "  MaxTrackEta = " << maxTrackEta_ << "\n";
0091   stream << "  MaxTrackDxy = " << maxTrackDxy_ << "\n";
0092   stream << "  MaxTrackDz = " << maxTrackDz_ << "\n";
0093   stream << "  MinTrackHits = " << minTrackHits_ << "\n";
0094   stream << "  MaxTrackNormChi2 = " << maxTrackNormChi2_ << "\n";
0095   //   stream << "  MaxDzMuonTrack = " << maxDzMuonTrack_ << "\n";
0096   LogDebug("HLTMuonTrackMassFilter") << stream.str();
0097 }
0098
0099 void HLTMuonTrackMassFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0100   edm::ParameterSetDescription desc;
0101   makeHLTFilterDescription(desc);
0102   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOfflineBeamSpot"));
0103   desc.add<edm::InputTag>("CandTag", edm::InputTag("hltL3MuonCandidates"));
0104   //  desc.add<edm::InputTag>("TrackTag",edm::InputTag("hltMuTkMuJpsiTrackerMuonCands"));
0105   desc.add<edm::InputTag>("TrackTag", edm::InputTag(""));
0106   //  desc.add<edm::InputTag>("PreviousCandTag",edm::InputTag("hltMu0TkMuJpsiTrackMassFiltered"));
0107   desc.add<edm::InputTag>("PreviousCandTag", edm::InputTag(""));
0108   {
0109     std::vector<double> temp1;
0110     temp1.reserve(1);
0111     temp1.push_back(2.8);
0112     desc.add<std::vector<double> >("MinMasses", temp1);
0113   }
0114   {
0115     std::vector<double> temp1;
0116     temp1.reserve(1);
0117     temp1.push_back(3.4);
0118     desc.add<std::vector<double> >("MaxMasses", temp1);
0119   }
0120   desc.add<bool>("checkCharge", true);
0121   desc.add<double>("MinTrackPt", 0.0);
0122   desc.add<double>("MinTrackP", 3.0);
0123   desc.add<double>("MaxTrackEta", 999.0);
0124   desc.add<double>("MaxTrackDxy", 999.0);
0125   desc.add<double>("MaxTrackDz", 999.0);
0126   desc.add<int>("MinTrackHits", 5);
0127   desc.add<double>("MaxTrackNormChi2", 10000000000.0);
0128   desc.add<double>("MaxDCAMuonTrack", 99999.9);
0129   desc.add<bool>("CutCowboys", false);
0130   descriptions.add("hltMuonTrackMassFilter", desc);
0131 }
0132
0133 bool HLTMuonTrackMassFilter::hltFilter(edm::Event& iEvent,
0134                                        const edm::EventSetup& iSetup,
0135                                        trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0136   // The filter object
0137   if (saveTags()) {
0138     filterproduct.addCollectionTag(muonTag_);
0139     filterproduct.addCollectionTag(trackTag_);
0140   }
0141   //
0142   // Beamspot
0143   //
0144   edm::Handle<reco::BeamSpot> beamspotHandle;
0145   iEvent.getByToken(beamspotToken_, beamspotHandle);
0146   reco::BeamSpot::Point beamspot(beamspotHandle->position());
0147   // Needed for DCA calculation
0148   auto const& bFieldHandle = iSetup.getHandle(idealMagneticFieldRecordToken_);
0149   //
0150   // Muons
0151   //
0152   edm::Handle<reco::RecoChargedCandidateCollection> muonHandle;
0153   iEvent.getByToken(muonToken_, muonHandle);
0154   //
0155   // Tracks
0156   //
0157   edm::Handle<reco::RecoChargedCandidateCollection> trackHandle;
0158   iEvent.getByToken(trackToken_, trackHandle);
0159   //
0160   // Muons from previous filter
0161   //
0162   edm::Handle<trigger::TriggerFilterObjectWithRefs> prevCandHandle;
0163   iEvent.getByToken(prevCandToken_, prevCandHandle);
0164   std::vector<reco::RecoChargedCandidateRef> prevMuonRefs;
0165   prevCandHandle->getObjects(trigger::TriggerMuon, prevMuonRefs);
0166   std::vector<reco::RecoChargedCandidateRef> prevTrackRefs;
0167   prevCandHandle->getObjects(trigger::TriggerTrack, prevTrackRefs);
0168   bool checkPrevTracks = prevTrackRefs.size() == prevMuonRefs.size();
0169   //   LogDebug("HLTMuonTrackMassFilter") << "#previous track refs = " << prevTrackRefs.size();
0170   //
0171   // access to muons and selection according to configuration
0172   //   if the previous candidates are taken from a muon+track
0173   //   quarkonia filter we rely on the fact that only the muons
0174   //   are flagged as TriggerMuon
0175   //
0176   std::vector<reco::RecoChargedCandidateRef> selectedMuonRefs;
0177   selectedMuonRefs.reserve(muonHandle->size());
0178   //   std::vector<size_t> prevMuonIndices;
0179   //   std::ostringstream stream1;
0180   for (unsigned int i = 0; i < muonHandle->size(); ++i) {
0181     // Ref
0182     reco::RecoChargedCandidateRef muonRef(muonHandle, i);
0183     //     stream1 << "Checking muon with q / pt / p / eta = "
0184     //      << muonRef->charge() << " " << muonRef->pt() << " "
0185     //      << muonRef->p() << " " << muonRef->eta() << "\n";
0186     // passed previous filter?
0187     if (find(prevMuonRefs.begin(), prevMuonRefs.end(), muonRef) == prevMuonRefs.end())
0188       continue;
0189     //     prevMuonIndices.push_back(find(prevMuonRefs.begin(),prevMuonRefs.end(),muonRef)-
0190     //                prevMuonRefs.begin());
0191     // keep muon
0192     //     stream1 << "... accepted as #" << selectedMuonRefs.size() << "\n";
0193     selectedMuonRefs.push_back(muonRef);
0194   }
0195   //   LogDebug("HLTMuonTrackMassFilter") << stream1.str();
0196   //
0197   // access to tracks and selection according to configuration
0198   //
0199   std::vector<reco::RecoChargedCandidateRef> selectedTrackRefs;
0200   selectedTrackRefs.reserve(trackHandle->size());
0201   //   std::ostringstream stream2;
0202   for (unsigned int i = 0; i < trackHandle->size(); ++i) {
0203     // validity of REF
0204     reco::RecoChargedCandidateRef trackRef(trackHandle, i);
0205     const reco::RecoChargedCandidate& trackCand = *trackRef;
0206     //     stream2 << "Checking track with q / pt / p / eta = "
0207     //      << trackCand.charge() << " " << trackCand.pt() << " "
0208     //      << trackCand.p() << " " << trackCand.eta() << "\n";
0209     // cuts on the momentum
0210     if (trackCand.pt() < minTrackPt_ || trackCand.p() < minTrackP_ || fabs(trackCand.eta()) > maxTrackEta_)
0211       continue;
0212     if (trackCand.track().isNull())
0213       continue;
0214     // cuts on track quality
0215     const reco::Track& track = *trackCand.track();
0216     //     stream2 << "... with dxy / dz / #hits / chi2 = "
0217     //      << track.dxy(beamspot) << " "
0218     //      << track.dz(beamspot) << " "
0219     //      << track.numberOfValidHits() << " "
0220     //      << track.normalizedChi2();
0221     if (fabs(track.dxy(beamspot)) > maxTrackDxy_ || fabs(track.dz(beamspot)) > maxTrackDz_ ||
0222         track.numberOfValidHits() < minTrackHits_ || track.normalizedChi2() > maxTrackNormChi2_)
0223       continue;
0224     // keep track
0225     //     stream2 << "... accepted as #" << selectedTrackRefs.size() << "\n";
0226     selectedTrackRefs.push_back(trackRef);
0227   }
0228   //   LogDebug("HLTMuonTrackMassFilter") << stream2.str();
0229   //
0230   // combinations
0231   //
0232   //   unsigned int nDz(0);
0233   unsigned int nQ(0);
0234   unsigned int nCowboy(0);
0235   unsigned int nComb(0);
0236   reco::Particle::LorentzVector p4Muon;
0237   reco::Particle::LorentzVector p4JPsi;
0238   //   std::ostringstream stream3;
0239   for (auto& selectedMuonRef : selectedMuonRefs) {
0240     const reco::RecoChargedCandidate& muon = *selectedMuonRef;
0241     int qMuon = muon.charge();
0242     p4Muon = muon.p4();
0243     for (auto& selectedTrackRef : selectedTrackRefs) {
0244       const reco::RecoChargedCandidate& track = *selectedTrackRef;
0245       //       stream3 << "Combination " << im << " / " << it << " with dz / q / mass = "
0246       //          << muon.track()->dz(beamspot)-track.track()->dz(beamspot) << " "
0247       //          << track.charge()+qMuon << " "
0248       //          << (p4Muon+track.p4()).mass() << "\n";
0249
0250       //       if ( fabs(muon.track()->dz(beamspot)-track.track()->dz(beamspot))>
0251       //       maxDzMuonTrack_ )  continue;
0252       //       ++nDz;
0253       if (checkCharge_ && track.charge() != -qMuon)
0254         continue;
0255       ++nQ;
0256
0257       ///
0258
0259       // DCA between the two muons
0260
0261       reco::TrackRef tk1 = muon.track();
0262       reco::TrackRef tk2 = track.track();
0263
0264       reco::TransientTrack mu1TT(*tk1, &(*bFieldHandle));
0265       reco::TransientTrack mu2TT(*tk2, &(*bFieldHandle));
0266       TrajectoryStateClosestToPoint mu1TS = mu1TT.impactPointTSCP();
0267       TrajectoryStateClosestToPoint mu2TS = mu2TT.impactPointTSCP();
0268       if (mu1TS.isValid() && mu2TS.isValid()) {
0269         ClosestApproachInRPhi cApp;
0270         cApp.calculate(mu1TS.theState(), mu2TS.theState());
0271         if (!cApp.status() || cApp.distance() > max_DCAMuonTrack_)
0272           continue;
0273       }
0274
0275       ///
0276       // if cutting on cowboys reject muons that bend towards each other
0277       if (cutCowboys_ && (qMuon * deltaPhi(p4Muon.phi(), track.phi()) > 0.))
0278         continue;
0279       ++nCowboy;
0280
0281       if (checkPrevTracks) {
0282         if (!pairMatched(prevMuonRefs, prevTrackRefs, selectedMuonRef, selectedTrackRef))
0283           continue;
0284       }
0285       double mass = (p4Muon + track.p4()).mass();
0286       for (unsigned int j = 0; j < minMasses_.size(); ++j) {
0287         if (mass > minMasses_[j] && mass < maxMasses_[j]) {
0288           ++nComb;
0289           filterproduct.addObject(trigger::TriggerMuon, selectedMuonRef);
0290           filterproduct.addObject(trigger::TriggerTrack, selectedTrackRef);
0291           //      stream3 << "... accepted\n";
0292           break;
0293         }
0294       }
0295     }
0296   }
0297   //   LogDebug("HLTMuonTrackMassFilter") << stream3.str();
0298
0299   if (edm::isDebugEnabled()) {
0300     std::ostringstream stream;
0301     stream << "Total number of combinations = "
0302            //      << selectedMuonRefs.size()*selectedTrackRefs.size() << " , after dz " << nDz
0303            << " , after charge " << nQ << " , after cutCowboy " << nCowboy << " , after mass " << nComb << std::endl;
0304     stream << "Found " << nComb << " jpsi candidates with # / mass / q / pt / eta" << std::endl;
0305     std::vector<reco::RecoChargedCandidateRef> muRefs;
0306     std::vector<reco::RecoChargedCandidateRef> tkRefs;
0307     filterproduct.getObjects(trigger::TriggerMuon, muRefs);
0308     filterproduct.getObjects(trigger::TriggerTrack, tkRefs);
0309     reco::Particle::LorentzVector p4Mu;
0310     reco::Particle::LorentzVector p4Tk;
0311     reco::Particle::LorentzVector p4JPsi;
0312     if (muRefs.size() == tkRefs.size()) {
0313       for (unsigned int i = 0; i < muRefs.size(); ++i) {
0314         p4Mu = muRefs[i]->p4();
0315         p4Tk = tkRefs[i]->p4();
0316         p4JPsi = p4Mu + p4Tk;
0317         stream << "  " << i << " " << p4JPsi.M() << " " << muRefs[i]->charge() + tkRefs[i]->charge() << " "
0318                << p4JPsi.P() << " " << p4JPsi.Eta() << "\n";
0319       }
0320       LogDebug("HLTMuonTrackMassFilter") << stream.str();
0321     } else {
0322       LogDebug("HLTMuonTrackMassFilter") << "different sizes for muon and track containers!!!";
0323     }
0324   }
0325
0326   return nComb > 0;
0327 }
0328
0329 bool HLTMuonTrackMassFilter::pairMatched(std::vector<reco::RecoChargedCandidateRef>& prevMuonRefs,
0330                                          std::vector<reco::RecoChargedCandidateRef>& prevTrackRefs,
0331                                          const reco::RecoChargedCandidateRef& muonRef,
0332                                          const reco::RecoChargedCandidateRef& trackRef) const {
0333   //
0334   // check only if references to tracks are available
0335   //
0336   if (prevTrackRefs.empty())
0337     return true;
0338   //
0339   // validity
0340   //
0341   if (prevMuonRefs.size() != prevTrackRefs.size())
0342     return false;
0343   edm::RefToBase<TrajectorySeed> seedRef = trackRef->track()->seedRef();
0344   if (seedRef.isNull())
0345     return false;
0346   //
0347   // comparison by hits of TrajectorySeed of the new track
0348   // with the previous candidate
0349   //
0350   const TrajectorySeed::RecHitRange seedHits = seedRef->recHits();
0351   trackingRecHit_iterator prevTrackHitEnd;
0352   trackingRecHit_iterator iprev;
0353   for (size_t i = 0; i < prevMuonRefs.size(); ++i) {
0354     // identity of muon
0355     if (prevMuonRefs[i] != muonRef)
0356       continue;
0357     // validity of Ref to previous track candidate
0358     reco::TrackRef prevTrackRef = prevTrackRefs[i]->track();
0359     if (prevTrackRef.isNull())
0360       continue;
0361     // if the references are the same then found and return true otherwise compare by hits
0362     if (prevTrackRef == trackRef->track())
0363       return true;
0364     // same #hits
0365     if (seedRef->nHits() != prevTrackRef->recHitsSize())
0366       continue;
0367     // hit-by-hit comparison based on the sharesInput method
0368     auto iseed = seedHits.begin();
0369     iprev = prevTrackRef->recHitsBegin();
0370     prevTrackHitEnd = prevTrackRef->recHitsEnd();
0371     bool identical(true);
0372     for (; iseed != seedHits.end() && iprev != prevTrackHitEnd; ++iseed, ++iprev) {
0373       if ((*iseed).isValid() != (**iprev).isValid() || !(*iseed).sharesInput(&**iprev, TrackingRecHit::all)) {
0374         // terminate loop over hits on first mismatch
0375         identical = false;
0376         break;
0377       }
0378     }
0379     // if seed and previous track candidates are identical : return success
0380     if (identical)
0381       return true;
0382   }
0383   // no match found
0384   return false;
0385 }
0386
0387 // declare this class as a framework plugin
0388 DEFINE_FWK_MODULE(HLTMuonTrackMassFilter);