Muon/plugins/HLTMuonL3andL2PreFilter.cc

0001 /** \class HLTMuonL3andL2PreFilter
0002  *
0003  *
0004  *  This class is an HLTFilter (-> EDFilter) implementing
0005  *  a simple filter to select L3 muons and L2 if they dont match with any L3
0006  *
0007  *  Original author:  A. Soto <alejandro.soto.rodriguez@cern.ch>
0008  *
0009  */
0010 #include "DataFormats/TrackReco/interface/Track.h"
0011 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0012 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 #include "HLTrigger/HLTcore/interface/HLTFilter.h"
0015 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0020
0021 #include <vector>
0022
0023 class HLTMuonL3andL2PreFilter : public HLTFilter {
0024 public:
0025   explicit HLTMuonL3andL2PreFilter(const edm::ParameterSet&);
0026   ~HLTMuonL3andL2PreFilter() override = default;
0027
0028   static void fillDescriptions(edm::ConfigurationDescriptions&);
0029   bool hltFilter(edm::Event&, const edm::EventSetup&, trigger::TriggerFilterObjectWithRefs&) const override;
0030
0031 private:
0032   class MuonSelectionCuts {
0033   public:
0034     explicit MuonSelectionCuts(const edm::ParameterSet&);
0035
0036     static void fillPSetDescription(edm::ParameterSetDescription&);
0037
0038     const int min_N;                   // minimum number of muons to fire the trigger
0039     const double max_DXYBeamSpot;      // cutoff in dxy from the beamspot
0040     const double min_DXYBeamSpot;      // minimum cut on dxy from the beamspot
0041     const double max_Dz;               // dz cut
0042     const double min_DxySig;           // dxy significance cut
0043     const double nsigma_Pt;            // pt uncertainty margin (in number of sigmas)
0044     const double max_PtDifference;     // cutoff in maximum different between global track and tracker track
0045     const double min_TrackPt;          // cutoff in tracker track pt
0046     const double max_NormalizedChi2;   // cutoff in normalized chi2
0047     const double min_Pt;               // pt threshold in GeV
0048     const double max_Eta;              // Eta cut
0049     const int min_NmuonHits;           // cutoff in minumum number of chi2 hits
0050     const std::vector<int> min_Nhits;  // threshold on number of hits on muon
0051     const std::vector<double>
0052         absetaBins;  // |eta| bins for minNstations cut (#bins must match #minNstations cuts and #minNhits cuts)
0053     const std::vector<int> minNstations;  // minimum number of muon stations used
0054     const std::vector<int> minNchambers;  // minimum number of valid chambers
0055   };
0056
0057   bool triggeredByPreviousLevel(const reco::RecoChargedCandidateRef&,
0058                                 const std::vector<reco::RecoChargedCandidateRef>&) const;
0059
0060   bool checkOverlap(const reco::RecoChargedCandidateRef&,
0061                     const std::vector<reco::RecoChargedCandidateRef>&,
0062                     const double) const;
0063
0064   bool applySelection(const reco::RecoChargedCandidateRef&, const reco::BeamSpot&, const MuonSelectionCuts&) const;
0065
0066   const edm::InputTag L3CandTag_;                                             // input tag identifying L3 muon container
0067   const edm::EDGetTokenT<reco::RecoChargedCandidateCollection> L3CandToken_;  // token identifying L3 muon container
0068
0069   const edm::InputTag L2CandTag_;                                             // input tag identifying L2 muon container
0070   const edm::EDGetTokenT<reco::RecoChargedCandidateCollection> L2CandToken_;  // token identifying L2 muon container
0071
0072   const edm::InputTag previousCandTag_;  // input tag identifying product contains muons passing the previous level
0073   const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs>
0074       previousCandToken_;  // token identifying product contains muons passing the previous level
0075
0076   const edm::InputTag beamspotTag_;                       // input tag identifying beamSpot container
0077   const edm::EDGetTokenT<reco::BeamSpot> beamspotToken_;  // token identifying beamSpot container
0078
0079   const bool matchPreviousCand_;  // flag to match L2 or L3 candidates to the previous level
0080   const double maxDeltaR_L2L3_;   // maximum deltaR between L2 and L3 muon to consider them the same
0081
0082   const MuonSelectionCuts cutParsL3_;  // container of all parameters to cut for L3
0083   const MuonSelectionCuts cutParsL2_;  // container of all parameters to cut for L2
0084 };
0085
0086 HLTMuonL3andL2PreFilter::MuonSelectionCuts::MuonSelectionCuts(const edm::ParameterSet& iConfig)
0087     : min_N(iConfig.getParameter<int>("MinN")),
0088       max_DXYBeamSpot(iConfig.getParameter<double>("MaxDXYBeamSpot")),
0089       min_DXYBeamSpot(iConfig.getParameter<double>("MinDXYBeamSpot")),
0090       max_Dz(iConfig.getParameter<double>("MaxDz")),
0091       min_DxySig(iConfig.getParameter<double>("MinDxySig")),
0092       nsigma_Pt(iConfig.getParameter<double>("NSigmaPt")),
0093       max_PtDifference(iConfig.getParameter<double>("MaxPtDifference")),
0094       min_TrackPt(iConfig.getParameter<double>("MinTrackPt")),
0095       max_NormalizedChi2(iConfig.getParameter<double>("MaxNormalizedChi2")),
0096       min_Pt(iConfig.getParameter<double>("MinPt")),
0097       max_Eta(iConfig.getParameter<double>("MaxEta")),
0098       min_NmuonHits(iConfig.getParameter<int>("MinNmuonHits")),
0099       min_Nhits(iConfig.getParameter<std::vector<int>>("MinNhits")),
0100       absetaBins(iConfig.getParameter<std::vector<double>>("AbsEtaBins")),
0101       minNstations(iConfig.getParameter<std::vector<int>>("MinNstations")),
0102       minNchambers(iConfig.getParameter<std::vector<int>>("MinNchambers")) {}
0103
0104 HLTMuonL3andL2PreFilter::HLTMuonL3andL2PreFilter(const edm::ParameterSet& iConfig)
0105     : HLTFilter(iConfig),
0106       L3CandTag_(iConfig.getParameter<edm::InputTag>("L3CandTag")),
0107       L3CandToken_(consumes(L3CandTag_)),
0108       L2CandTag_(iConfig.getParameter<edm::InputTag>("L2CandTag")),
0109       L2CandToken_(consumes(L2CandTag_)),
0110       previousCandTag_(iConfig.getParameter<edm::InputTag>("PreviousCandTag")),
0111       previousCandToken_(consumes(previousCandTag_)),
0112       beamspotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0113       beamspotToken_(consumes(beamspotTag_)),
0114       matchPreviousCand_(iConfig.getParameter<bool>("MatchToPreviousCand")),
0115       maxDeltaR_L2L3_(iConfig.getParameter<double>("MaxDeltaRL2L3")),
0116       cutParsL3_(iConfig.getParameter<edm::ParameterSet>("L3CandSelection")),
0117       cutParsL2_(iConfig.getParameter<edm::ParameterSet>("L2CandSelection")) {
0118   // check consistency of parameters
0119   if (cutParsL2_.absetaBins.size() != cutParsL2_.minNstations.size() ||
0120       cutParsL2_.absetaBins.size() != cutParsL2_.minNchambers.size() ||
0121       cutParsL2_.absetaBins.size() != cutParsL2_.min_Nhits.size()) {
0122     throw cms::Exception("Configuration")
0123         << "error in ParameterSet \"L2CandSelection\": size of \"AbsEtaBins\" (" << cutParsL2_.absetaBins.size()
0124         << "), \"MinNstations\" (" << cutParsL2_.minNstations.size() << "), \"MinNchambers\" ("
0125         << cutParsL2_.minNchambers.size() << ") and \"MinNhits\" (" << cutParsL2_.min_Nhits.size() << ") differ";
0126   }
0127   if (cutParsL3_.absetaBins.size() != cutParsL3_.minNstations.size() ||
0128       cutParsL3_.absetaBins.size() != cutParsL3_.minNchambers.size() ||
0129       cutParsL3_.absetaBins.size() != cutParsL3_.min_Nhits.size()) {
0130     throw cms::Exception("Configuration")
0131         << "error in ParameterSet \"L3CandSelection\": size of \"AbsEtaBins\" (" << cutParsL3_.absetaBins.size()
0132         << "), \"MinNstations\" (" << cutParsL3_.minNstations.size() << "), \"MinNchambers\" ("
0133         << cutParsL3_.minNchambers.size() << ") and \"MinNhits\" (" << cutParsL3_.min_Nhits.size() << ") differ";
0134   }
0135
0136   if (maxDeltaR_L2L3_ <= 0.) {
0137     throw cms::Exception("Configuration")
0138         << "invalid value for parameter \"MaxDeltaRL2L3\" (must be > 0): " << maxDeltaR_L2L3_;
0139   }
0140 }
0141
0142 void HLTMuonL3andL2PreFilter::MuonSelectionCuts::fillPSetDescription(edm::ParameterSetDescription& desc) {
0143   desc.add<int>("MinN", 1);
0144   desc.add<double>("MaxDXYBeamSpot", 9999.0);
0145   desc.add<double>("MinDXYBeamSpot", -1.0);
0146   desc.add<double>("MaxDz", 9999.0);
0147   desc.add<double>("MinDxySig", -1.0);
0148   desc.add<double>("NSigmaPt", 0.0);
0149   desc.add<double>("MaxPtDifference", 9999.0);
0150   desc.add<double>("MinTrackPt", 0.0);
0151   desc.add<double>("MaxNormalizedChi2", 9999.0);
0152   desc.add<double>("MinPt", 23.0);
0153   desc.add<double>("MaxEta", 2.5);
0154   desc.add<int>("MinNmuonHits", 0);
0155   desc.add<std::vector<int>>("MinNhits", {1, 0});
0156   desc.add<std::vector<double>>("AbsEtaBins", {1, 9999.});
0157   desc.add<std::vector<int>>("MinNstations", {1, 1});
0158   desc.add<std::vector<int>>("MinNchambers", {1, 0});
0159 };
0160
0161 void HLTMuonL3andL2PreFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0162   edm::ParameterSetDescription desc;
0163   makeHLTFilterDescription(desc);
0164   desc.add<edm::InputTag>("L3CandTag", edm::InputTag("hltL3MuonCandidates"));
0165   desc.add<edm::InputTag>("L2CandTag", edm::InputTag("hltL2MuonCandidates"));
0166   desc.add<edm::InputTag>("PreviousCandTag", edm::InputTag(""));
0167   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOnlineBeamSpot"));
0168   desc.add<bool>("MatchToPreviousCand", true);
0169   desc.add<double>("MaxDeltaRL2L3", 0.05);
0170
0171   edm::ParameterSetDescription descMuonSelL2;
0172   MuonSelectionCuts::fillPSetDescription(descMuonSelL2);
0173   desc.add<edm::ParameterSetDescription>("L2CandSelection", descMuonSelL2);
0174
0175   edm::ParameterSetDescription descMuonSelL3;
0176   MuonSelectionCuts::fillPSetDescription(descMuonSelL3);
0177   desc.add<edm::ParameterSetDescription>("L3CandSelection", descMuonSelL3);
0178
0179   descriptions.addWithDefaultLabel(desc);
0180 }
0181
0182 bool HLTMuonL3andL2PreFilter::hltFilter(edm::Event& iEvent,
0183                                         const edm::EventSetup& iSetup,
0184                                         trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0185   // All HLT filters must create and fill an HLT filter object,
0186   // recording any reconstructed physics objects satisfying (or not)
0187   // this HLT filter, and place it in the Event.
0188   if (saveTags()) {
0189     filterproduct.addCollectionTag(L3CandTag_);
0190     filterproduct.addCollectionTag(L2CandTag_);
0191   }
0192   // get hold of trks
0193   auto const L3mucands = iEvent.getHandle(L3CandToken_);
0194   auto const L2mucands = iEvent.getHandle(L2CandToken_);
0195   std::vector<reco::RecoChargedCandidateRef> vcands;
0196   if (matchPreviousCand_) {
0197     auto const& previousLevelCands = iEvent.get(previousCandToken_);
0198     previousLevelCands.getObjects(trigger::TriggerMuon, vcands);
0199   }
0200
0201   auto const& beamSpot = iEvent.get(beamspotToken_);
0202
0203   // Number of objects passing the L3 Trigger:
0204   int nL3 = 0;
0205   for (unsigned int iMuL3 = 0; iMuL3 < L3mucands->size(); iMuL3++) {
0206     reco::RecoChargedCandidateRef L3cand(L3mucands, iMuL3);
0207     LogDebug("HLTMuonL3andL2PreFilter") << "L3cand isNonnull " << L3cand.isNonnull();
0208     // did this candidate triggered at previous stage.
0209     if (matchPreviousCand_ && !triggeredByPreviousLevel(L3cand, vcands))
0210       continue;
0211     // apply all the cuts
0212     if (!applySelection(L3cand, beamSpot, cutParsL3_))
0213       continue;
0214     nL3++;
0215     filterproduct.addObject(trigger::TriggerMuon, L3cand);
0216   }
0217   // get the L3 muons that passed the selection
0218   std::vector<reco::RecoChargedCandidateRef> passingL3mucands;
0219   filterproduct.getObjects(trigger::TriggerMuon, passingL3mucands);
0220
0221   // Number of objects passing the L2 Trigger:
0222   int nL2 = 0;
0223   auto const maxDeltaR2_L2L3 = maxDeltaR_L2L3_ * maxDeltaR_L2L3_;
0224   for (unsigned int iMuL2 = 0; iMuL2 < L2mucands->size(); iMuL2++) {
0225     reco::RecoChargedCandidateRef L2cand(L2mucands, iMuL2);
0226     LogDebug("HLTMuonL3andL2PreFilter") << "L2cand isNonnull " << L2cand.isNonnull();
0227     // did this candidate triggered at previous stage.
0228     if (matchPreviousCand_ && !triggeredByPreviousLevel(L2cand, vcands))
0229       continue;
0230     // check if L2 and L3 candidates are the same
0231     if (checkOverlap(L2cand, passingL3mucands, maxDeltaR2_L2L3)) {
0232       continue;
0233     }
0234     // apply all the cuts
0235     if (!applySelection(L2cand, beamSpot, cutParsL2_))
0236       continue;
0237     nL2++;
0238     filterproduct.addObject(trigger::TriggerMuon, L2cand);
0239   }
0240
0241   // filter decision
0242   const bool acceptL3(nL3 >= cutParsL3_.min_N);
0243   const bool acceptL2(nL2 >= cutParsL2_.min_N);
0244
0245   LogDebug("HLTMuonL3andL2PreFilter") << " >>>>> Result of HLTMuonL3andL2PreFilter is " << (acceptL3 && acceptL2)
0246                                       << ", number of muons passing thresholds= " << nL2 + nL3;
0247
0248   return acceptL2 && acceptL3;
0249 }
0250
0251 bool HLTMuonL3andL2PreFilter::triggeredByPreviousLevel(const reco::RecoChargedCandidateRef& candref,
0252                                                        const std::vector<reco::RecoChargedCandidateRef>& vcands) const {
0253   return (std::find(vcands.begin(), vcands.end(), candref) != vcands.end());
0254 }
0255
0256 bool HLTMuonL3andL2PreFilter::checkOverlap(
0257     const reco::RecoChargedCandidateRef& r1,
0258     const std::vector<reco::RecoChargedCandidateRef>& R2,  // vector with the recoCharged candidates from L3
0259     const double maxDeltaR2) const {
0260   for (auto const& r2 : R2) {
0261     if (reco::deltaR2(*r1, *r2) < maxDeltaR2)
0262       return true;
0263   }
0264   return false;
0265 }
0266
0267 bool HLTMuonL3andL2PreFilter::applySelection(const reco::RecoChargedCandidateRef& candidate,
0268                                              const reco::BeamSpot& beamSpot,
0269                                              const MuonSelectionCuts& cutPars) const {
0270   // eta cut
0271   auto const candAbsEta = std::abs(candidate->eta());
0272   if (candAbsEta > cutPars.max_Eta)
0273     return false;
0274
0275   reco::TrackRef const& tk = candidate->track();
0276   LogDebug("HLTMuonL3andL2PreFilter") << " Muon in loop, q*pt= " << tk->charge() * tk->pt() << " ("
0277                                       << candidate->charge() * candidate->pt() << ") "
0278                                       << ", eta= " << tk->eta() << " (" << candidate->eta() << ") "
0279                                       << ", hits= " << tk->numberOfValidHits() << ", d0= " << tk->d0()
0280                                       << ", dz= " << tk->dz();
0281
0282   // normalizedChi2 cut
0283   if (tk->normalizedChi2() > cutPars.max_NormalizedChi2)
0284     return false;
0285
0286   // number of eta bins for cut on number of stations
0287   const auto nAbsetaBins = cutPars.absetaBins.size();
0288
0289   // cut on number of stations
0290   bool failNstations(false), failNhits(false), failNchambers(false);
0291   for (unsigned int i = 0; i < nAbsetaBins; ++i) {
0292     if (candAbsEta < cutPars.absetaBins[i]) {
0293       if (candidate->track()->hitPattern().muonStationsWithAnyHits() < cutPars.minNstations[i]) {
0294         failNstations = true;
0295       }
0296       if (candidate->track()->numberOfValidHits() < cutPars.min_Nhits[i]) {
0297         failNhits = true;
0298       }
0299       if ((candidate->track()->hitPattern().dtStationsWithAnyHits() +
0300            candidate->track()->hitPattern().cscStationsWithAnyHits()) < cutPars.minNchambers[i]) {
0301         failNchambers = true;
0302       }
0303       break;
0304     }
0305   }
0306   if (failNstations || failNhits || failNchambers)
0307     return false;
0308
0309   // dz cut (impact parameter of the muon in the z direction computed with respect to the beamspot position)
0310   auto const vz = candidate->vz() - beamSpot.z0();
0311   auto const vx = candidate->vx() - beamSpot.x0();
0312   auto const vy = candidate->vy() - beamSpot.y0();
0313   if (candidate->pt() <= 0 or
0314       (std::abs(vz - (vx * candidate->px() + vy * candidate->py()) / (candidate->pt() * candidate->pt()) *
0315                          candidate->pz()) > cutPars.max_Dz))
0316     return false;
0317
0318   // dxy significance cut (safeguard against bizarre values)
0319   auto const absDxy = std::abs(tk->dxy(beamSpot.position()));
0320   if (cutPars.min_DxySig > 0 && (tk->dxyError() <= 0 || (absDxy < cutPars.min_DxySig * tk->dxyError())))
0321     return false;
0322
0323   // dxy beamspot cut
0324   if (absDxy > cutPars.max_DXYBeamSpot || absDxy < cutPars.min_DXYBeamSpot)
0325     return false;
0326
0327   // min muon hits cut
0328   const reco::HitPattern& trackHits = tk->hitPattern();
0329   if (trackHits.numberOfValidMuonHits() < cutPars.min_NmuonHits)
0330     return false;
0331
0332   // pt difference cut
0333   auto const candPt = candidate->pt();
0334   auto const trackPt = tk->pt();
0335
0336   if (std::abs(candPt - trackPt) > cutPars.max_PtDifference)
0337     return false;
0338
0339   // track pt cut
0340   if (trackPt < cutPars.min_TrackPt)
0341     return false;
0342
0343   // pt threshold cut
0344   auto const err0 = tk->error(0);
0345   auto const abspar0 = std::abs(tk->parameter(0));
0346   // if abspar0 > 0, convert 50% efficiency threshold to 90% efficiency threshold
0347   auto const ptLx = abspar0 > 0 ? candPt + cutPars.nsigma_Pt * err0 / abspar0 * candPt : candPt;
0348   if (ptLx < cutPars.min_Pt)
0349     return false;
0350
0351   return true;
0352 }
0353
0354 DEFINE_FWK_MODULE(HLTMuonL3andL2PreFilter);