Muon/plugins/HLTMuonDimuonL3Filter.cc

0001 /** \class HLTMuonDimuonL2Filter
0002  *
0003  * See header file for documentation
0004  *
0005  *  \author J. Alcaraz, P. Garcia
0006  *
0007  */
0008
0009 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0010 #include "DataFormats/HLTReco/interface/TriggerRefsCollections.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "DataFormats/TrackReco/interface/Track.h"
0013 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0014 #include "HLTMuonDimuonL3Filter.h"
0015 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
0016 #include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeedCollection.h"
0017 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0018 #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
0019 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "DataFormats/Math/interface/deltaR.h"
0022
0023 using namespace edm;
0024 using namespace std;
0025 using namespace reco;
0026 using namespace trigger;
0027
0028 //
0029 // constructors and destructor
0030 //
0031 namespace {
0032   struct Out {
0033     Out(std::vector<double> const& v) : v_(v) {}
0034
0035     std::vector<double> const& v_;
0036   };
0037
0038 #if defined(EDM_ML_DEBUG)
0039   std::ostream& operator<<(std::ostream& iS, Out const& iO) {
0040     iS << "[";
0041     for (double v : iO.v_) {
0042       iS << v << " ";
0043     }
0044     iS << "]";
0045     return iS;
0046   }
0047 #endif
0048 }  // namespace
0049
0050 HLTMuonDimuonL3Filter::HLTMuonDimuonL3Filter(const edm::ParameterSet& iConfig)
0051     : HLTFilter(iConfig),
0052       propSetup_(iConfig, consumesCollector()),
0053       idealMagneticFieldRecordToken_(esConsumes()),
0054       beamspotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0055       beamspotToken_(consumes<reco::BeamSpot>(beamspotTag_)),
0056       candTag_(iConfig.getParameter<edm::InputTag>("CandTag")),
0057       candToken_(consumes<reco::RecoChargedCandidateCollection>(candTag_)),
0058       previousCandTag_(iConfig.getParameter<InputTag>("PreviousCandTag")),
0059       previousCandToken_(consumes<trigger::TriggerFilterObjectWithRefs>(previousCandTag_)),
0060       l1CandTag_(iConfig.getParameter<InputTag>("L1CandTag")),
0061       l1CandToken_(consumes<trigger::TriggerFilterObjectWithRefs>(l1CandTag_)),
0062       recoMuTag_(iConfig.getParameter<InputTag>("inputMuonCollection")),
0063       recoMuToken_(consumes<reco::MuonCollection>(recoMuTag_)),
0064       previousCandIsL2_(iConfig.getParameter<bool>("PreviousCandIsL2")),
0065       fast_Accept_(iConfig.getParameter<bool>("FastAccept")),
0066       min_N_(iConfig.getParameter<int>("MinN")),
0067       max_Eta_(iConfig.getParameter<double>("MaxEta")),
0068       min_Nhits_(iConfig.getParameter<int>("MinNhits")),
0069       max_Dr_(iConfig.getParameter<double>("MaxDr")),
0070       max_Dz_(iConfig.getParameter<double>("MaxDz")),
0071       chargeOpt_(iConfig.getParameter<int>("ChargeOpt")),
0072       min_PtPair_(iConfig.getParameter<vector<double> >("MinPtPair")),
0073       max_PtPair_(iConfig.getParameter<vector<double> >("MaxPtPair")),
0074       min_PtMax_(iConfig.getParameter<vector<double> >("MinPtMax")),
0075       min_PtMin_(iConfig.getParameter<vector<double> >("MinPtMin")),
0076       max_PtMin_(iConfig.getParameter<vector<double> >("MaxPtMin")),
0077       min_InvMass_(iConfig.getParameter<vector<double> >("MinInvMass")),
0078       max_InvMass_(iConfig.getParameter<vector<double> >("MaxInvMass")),
0079       applyMinDiMuonDeltaR2Cut_(iConfig.getParameter<double>("MinDiMuonDeltaR") > 0.),
0080       min_DiMuonDeltaR2_(iConfig.getParameter<double>("MinDiMuonDeltaR") *
0081                          iConfig.getParameter<double>("MinDiMuonDeltaR")),
0082       min_Acop_(iConfig.getParameter<double>("MinAcop")),
0083       max_Acop_(iConfig.getParameter<double>("MaxAcop")),
0084       min_PtBalance_(iConfig.getParameter<double>("MinPtBalance")),
0085       max_PtBalance_(iConfig.getParameter<double>("MaxPtBalance")),
0086       nsigma_Pt_(iConfig.getParameter<double>("NSigmaPt")),
0087       max_DCAMuMu_(iConfig.getParameter<double>("MaxDCAMuMu")),
0088       max_YPair_(iConfig.getParameter<double>("MaxRapidityPair")),
0089       cutCowboys_(iConfig.getParameter<bool>("CutCowboys")),
0090       theL3LinksLabel(iConfig.getParameter<InputTag>("InputLinks")),
0091       linkToken_(consumes<reco::MuonTrackLinksCollection>(theL3LinksLabel)),
0092       L1MatchingdR_(iConfig.getParameter<double>("L1MatchingdR")),
0093       L1MatchingdR2_(L1MatchingdR_ * L1MatchingdR_),
0094       matchPreviousCand_(iConfig.getParameter<bool>("MatchToPreviousCand")),
0095       MuMass2_(0.106 * 0.106) {
0096   // check consistency of parameters for mass-window cuts
0097   if (min_InvMass_.size() != min_PtPair_.size()) {
0098     throw cms::Exception("Configuration") << "size of \"MinInvMass\" (" << min_InvMass_.size()
0099                                           << ") and \"MinPtPair\" (" << min_PtPair_.size() << ") differ";
0100   }
0101   if (min_InvMass_.size() != max_PtPair_.size()) {
0102     throw cms::Exception("Configuration") << "size of \"MinInvMass\" (" << min_InvMass_.size()
0103                                           << ") and \"MaxPtPair\" (" << max_PtPair_.size() << ") differ";
0104   }
0105   if (min_InvMass_.size() != min_PtMax_.size()) {
0106     throw cms::Exception("Configuration") << "size of \"MinInvMass\" (" << min_InvMass_.size() << ") and \"MinPtMax\" ("
0107                                           << min_PtMax_.size() << ") differ";
0108   }
0109   if (min_InvMass_.size() != min_PtMin_.size()) {
0110     throw cms::Exception("Configuration") << "size of \"MinInvMass\" (" << min_InvMass_.size() << ") and \"MinPtMin\" ("
0111                                           << min_PtMin_.size() << ") differ";
0112   }
0113   if (min_InvMass_.size() != max_PtMin_.size()) {
0114     throw cms::Exception("Configuration") << "size of \"MinInvMass\" (" << min_InvMass_.size() << ") and \"MaxPtMin\" ("
0115                                           << max_PtMin_.size() << ") differ";
0116   }
0117   if (min_InvMass_.size() != max_InvMass_.size()) {
0118     throw cms::Exception("Configuration") << "size of \"MinInvMass\" (" << min_InvMass_.size()
0119                                           << ") and \"MaxInvMass\" (" << max_InvMass_.size() << ") differ";
0120   }
0121
0122   if (L1MatchingdR_ <= 0.) {
0123     throw cms::Exception("HLTMuonDimuonL3FilterConfiguration")
0124         << "invalid value for parameter \"L1MatchingdR\" (must be > 0): " << L1MatchingdR_;
0125   }
0126   LogDebug("HLTMuonDimuonL3Filter") << " CandTag/FastAccept/MinN/MaxEta/MinNhits/MaxDr/MaxDz/MinPt1/MinPt2/MinInvMass/"
0127                                        "MaxInvMass/applyMinDiMuonDeltaRCut/MinDiMuonDeltaR"
0128                                        "MinAcop/MaxAcop/MinPtBalance/MaxPtBalance/NSigmaPt/MaxDzMuMu/MaxRapidityPair : "
0129                                     << candTag_.encode() << " " << fast_Accept_ << " " << min_N_ << " " << max_Eta_
0130                                     << " " << min_Nhits_ << " " << max_Dr_ << " " << max_Dz_ << " " << chargeOpt_ << " "
0131                                     << Out(min_PtPair_) << " " << Out(min_PtMax_) << " " << Out(min_PtMin_) << " "
0132                                     << Out(min_InvMass_) << " " << Out(max_InvMass_) << " " << applyMinDiMuonDeltaR2Cut_
0133                                     << " " << sqrt(min_DiMuonDeltaR2_) << " " << min_Acop_ << " " << max_Acop_ << " "
0134                                     << min_PtBalance_ << " " << max_PtBalance_ << " " << nsigma_Pt_ << " "
0135                                     << max_DCAMuMu_ << " " << max_YPair_;
0136 }
0137
0138 HLTMuonDimuonL3Filter::~HLTMuonDimuonL3Filter() = default;
0139
0140 void HLTMuonDimuonL3Filter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0141   edm::ParameterSetDescription desc;
0142   makeHLTFilterDescription(desc);
0143   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOfflineBeamSpot"));
0144   desc.add<edm::InputTag>("CandTag", edm::InputTag("hltL3MuonCandidates"));
0145   //  desc.add<edm::InputTag>("PreviousCandTag",edm::InputTag("hltDiMuonL2PreFiltered0"));
0146   desc.add<edm::InputTag>("PreviousCandTag", edm::InputTag(""));
0147   desc.add<edm::InputTag>("L1CandTag", edm::InputTag(""));
0148   desc.add<edm::InputTag>("inputMuonCollection", edm::InputTag(""));
0149   desc.add<bool>("PreviousCandIsL2", true);
0150   desc.add<bool>("FastAccept", false);
0151   desc.add<int>("MinN", 1);
0152   desc.add<double>("MaxEta", 2.5);
0153   desc.add<int>("MinNhits", 0);
0154   desc.add<double>("MaxDr", 2.0);
0155   desc.add<double>("MaxDz", 9999.0);
0156   desc.add<int>("ChargeOpt", 0);
0157   vector<double> v1;
0158   v1.push_back(0.0);
0159   vector<double> v2;
0160   v2.push_back(1e125);
0161   vector<double> v3;
0162   v3.push_back(3.0);
0163   vector<double> v4;
0164   v4.push_back(3.0);
0165   vector<double> v5;
0166   v5.push_back(1e125);
0167   vector<double> v6;
0168   v6.push_back(2.8);
0169   vector<double> v7;
0170   v7.push_back(3.4);
0171   desc.add<vector<double> >("MinPtPair", v1);
0172   desc.add<vector<double> >("MaxPtPair", v2);
0173   desc.add<vector<double> >("MinPtMax", v3);
0174   desc.add<vector<double> >("MinPtMin", v4);
0175   desc.add<vector<double> >("MaxPtMin", v5);
0176   desc.add<vector<double> >("MinInvMass", v6);
0177   desc.add<vector<double> >("MaxInvMass", v7);
0178   desc.add<double>("MinDiMuonDeltaR", -1.);
0179   desc.add<double>("MinAcop", -1.0);
0180   desc.add<double>("MaxAcop", 3.15);
0181   desc.add<double>("MinPtBalance", -1.0);
0182   desc.add<double>("MaxPtBalance", 999999.0);
0183   desc.add<double>("NSigmaPt", 0.0);
0184   desc.add<double>("MaxDCAMuMu", 99999.9);
0185   desc.add<double>("MaxRapidityPair", 999999.0);
0186   desc.add<bool>("CutCowboys", false);
0187   desc.add<edm::InputTag>("InputLinks", edm::InputTag(""));
0188   desc.add<double>("L1MatchingdR", 0.3);
0189   desc.add<bool>("MatchToPreviousCand", true);
0190   PropagateToMuonSetup::fillPSetDescription(desc);
0191   descriptions.add("hltMuonDimuonL3Filter", desc);
0192 }
0193
0194 //
0195 // member functions
0196 //
0197
0198 // ------------ method called to produce the data  ------------
0199 bool HLTMuonDimuonL3Filter::hltFilter(edm::Event& iEvent,
0200                                       const edm::EventSetup& iSetup,
0201                                       trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0202   // All HLT filters must create and fill an HLT filter object,
0203   // recording any reconstructed physics objects satisfying (or not)
0204   // this HLT filter, and place it in the Event.
0205
0206   auto const prop = propSetup_.init(iSetup);
0207
0208   // Read RecoChargedCandidates from L3MuonCandidateProducer:
0209   Handle<RecoChargedCandidateCollection> mucands;
0210   if (saveTags())
0211     filterproduct.addCollectionTag(candTag_);  //?
0212   iEvent.getByToken(candToken_, mucands);
0213
0214   // Read L2 triggered objects:
0215   Handle<TriggerFilterObjectWithRefs> previousLevelCands;
0216   iEvent.getByToken(previousCandToken_, previousLevelCands);
0217   vector<RecoChargedCandidateRef> vl2cands;
0218   previousLevelCands->getObjects(TriggerMuon, vl2cands);
0219
0220   // Read BeamSpot information:
0221   Handle<BeamSpot> recoBeamSpotHandle;
0222   iEvent.getByToken(beamspotToken_, recoBeamSpotHandle);
0223   const BeamSpot& beamSpot = *recoBeamSpotHandle;
0224
0225   // sort them by L2Track
0226   std::map<reco::TrackRef, std::vector<RecoChargedCandidateRef> > L2toL3s;
0227   // map the L3 cands matched to a L1 to their position in the recoMuon collection
0228   std::map<unsigned int, RecoChargedCandidateRef> MuonToL3s;
0229
0230   // Test to see if we can use L3MuonTrajectorySeeds:
0231   if (mucands->empty())
0232     return false;
0233   auto const& tk = (*mucands)[0].track();
0234   bool useL3MTS = false;
0235
0236   if (tk->seedRef().isNonnull()) {
0237     auto a = dynamic_cast<const L3MuonTrajectorySeed*>(tk->seedRef().get());
0238     useL3MTS = a != nullptr;
0239   }
0240
0241   // If we can use L3MuonTrajectory seeds run the older code:
0242   if (useL3MTS) {
0243     unsigned int maxI = mucands->size();
0244     for (unsigned int i = 0; i != maxI; i++) {
0245       const TrackRef& tk = (*mucands)[i].track();
0246       if (previousCandIsL2_) {
0247         edm::Ref<L3MuonTrajectorySeedCollection> l3seedRef =
0248             tk->seedRef().castTo<edm::Ref<L3MuonTrajectorySeedCollection> >();
0249         TrackRef staTrack = l3seedRef->l2Track();
0250         L2toL3s[staTrack].push_back(RecoChargedCandidateRef(mucands, i));
0251       } else {
0252         L2toL3s[tk].push_back(RecoChargedCandidateRef(mucands, i));
0253       }
0254     }
0255   }
0256   // Using normal TrajectorySeeds:
0257   else {
0258     // Read Links collection:
0259     edm::Handle<reco::MuonTrackLinksCollection> links;
0260     iEvent.getByToken(linkToken_, links);
0261
0262     edm::Handle<trigger::TriggerFilterObjectWithRefs> level1Cands;
0263     std::vector<l1t::MuonRef> vl1cands;
0264
0265     bool check_l1match = true;
0266
0267     // Loop over RecoChargedCandidates:
0268     for (unsigned int i(0); i < mucands->size(); ++i) {
0269       RecoChargedCandidateRef cand(mucands, i);
0270       TrackRef tk = cand->track();  // is inner track
0271
0272       if (!matchPreviousCand_) {
0273         MuonToL3s[i] = RecoChargedCandidateRef(cand);
0274       } else {
0275         check_l1match = true;
0276         for (auto const& link : *links) {
0277           // Using the same method that was used to create the links between L3 and L2
0278           // ToDo: there should be a better way than dR,dPt matching
0279           const reco::Track& trackerTrack = *link.trackerTrack();
0280           if (tk->pt() == 0 or trackerTrack.pt() == 0)
0281             continue;
0282
0283           float dR2 = deltaR2(tk->eta(), tk->phi(), trackerTrack.eta(), trackerTrack.phi());
0284           float dPt = std::abs(tk->pt() - trackerTrack.pt()) / tk->pt();
0285
0286           if (dR2 < 0.02 * 0.02 and dPt < 0.001) {
0287             const TrackRef staTrack = link.standAloneTrack();
0288             L2toL3s[staTrack].push_back(RecoChargedCandidateRef(cand));
0289             check_l1match = false;
0290           }
0291         }  //MTL loop
0292
0293         if (not l1CandTag_.label().empty() and check_l1match) {
0294           auto const propagated = prop.extrapolate(*tk);
0295           auto const etaForMatch = propagated.isValid() ? propagated.globalPosition().eta() : cand->eta();
0296           auto const phiForMatch = propagated.isValid() ? (double)propagated.globalPosition().phi() : cand->phi();
0297           iEvent.getByToken(l1CandToken_, level1Cands);
0298           level1Cands->getObjects(trigger::TriggerL1Mu, vl1cands);
0299           const unsigned int nL1Muons(vl1cands.size());
0300           for (unsigned int il1 = 0; il1 != nL1Muons; ++il1) {
0301             if (deltaR2(etaForMatch, phiForMatch, vl1cands[il1]->eta(), vl1cands[il1]->phi()) <
0302                 L1MatchingdR2_) {  //was muon, non cand
0303               MuonToL3s[i] = RecoChargedCandidateRef(cand);
0304             }
0305           }
0306         }
0307       }
0308     }  //RCC loop
0309   }    //end of using normal TrajectorySeeds
0310
0311   // Needed for DCA calculation
0312   auto const& bFieldHandle = iSetup.getHandle(idealMagneticFieldRecordToken_);
0313
0314   // look at all mucands,  check cuts and add to filter object
0315   int n = 0;
0316
0317   // look at all mucands,  check cuts and add to filter object
0318   auto L2toL3s_it1 = L2toL3s.begin();
0319   auto L2toL3s_end = L2toL3s.end();
0320   bool atLeastOnePair = false;
0321   for (; L2toL3s_it1 != L2toL3s_end; ++L2toL3s_it1) {
0322     if (!triggeredByLevel2(L2toL3s_it1->first, vl2cands))
0323       continue;
0324
0325     //loop over the L3Tk reconstructed for this L2.
0326     unsigned int iTk1 = 0;
0327     unsigned int maxItk1 = L2toL3s_it1->second.size();
0328     for (; iTk1 != maxItk1; iTk1++) {
0329       bool thisL3Index1isDone = false;
0330       RecoChargedCandidateRef& cand1 = L2toL3s_it1->second[iTk1];
0331       TrackRef tk1 = cand1->get<TrackRef>();
0332
0333       LogDebug("HLTMuonDimuonL3Filter") << " 1st muon in loop: q*pt= " << tk1->charge() * tk1->pt() << " ("
0334                                         << cand1->charge() * cand1->pt() << ") "
0335                                         << ", eta= " << tk1->eta() << " (" << cand1->eta() << ") "
0336                                         << ", hits= " << tk1->numberOfValidHits();
0337
0338       // Run muon selection on first muon:
0339       if (!applyMuonSelection(cand1, beamSpot))
0340         continue;
0341
0342       // Pt threshold cut
0343       // Don't convert to 90% efficiency threshold
0344       LogDebug("HLTMuonDimuonL3Filter") << " ... 1st muon in loop, pt1= " << cand1->pt();
0345
0346       // Loop on 2nd muon cand
0347       auto L2toL3s_it2 = L2toL3s_it1;
0348       L2toL3s_it2++;
0349       for (; L2toL3s_it2 != L2toL3s_end; ++L2toL3s_it2) {
0350         if (!triggeredByLevel2(L2toL3s_it2->first, vl2cands))
0351           continue;
0352
0353         //loop over the L3Tk reconstructed for this L2.
0354         unsigned int iTk2 = 0;
0355         unsigned int maxItk2 = L2toL3s_it2->second.size();
0356         for (; iTk2 != maxItk2; iTk2++) {
0357           RecoChargedCandidateRef& cand2 = L2toL3s_it2->second[iTk2];
0358           TrackRef tk2 = cand2->get<TrackRef>();
0359
0360           LogDebug("HLTMuonDimuonL3Filter") << " 2nd muon in loop: q*pt= " << tk2->charge() * tk2->pt() << " ("
0361                                             << cand2->charge() * cand2->pt() << ") "
0362                                             << ", eta= " << tk2->eta() << " (" << cand2->eta() << ") "
0363                                             << ", hits= " << tk2->numberOfValidHits() << ", d0= " << tk2->d0();
0364
0365           // Run muon selection on second muon:
0366           if (!applyMuonSelection(cand2, beamSpot))
0367             continue;
0368
0369           // Pt threshold cut
0370           // Don't convert to 90% efficiency threshold
0371           LogDebug("HLTMuonDimuonL3Filter") << " ... 2nd muon in loop, pt2= " << cand2->pt();
0372
0373           // Run dimuon selection:
0374           if (!applyDiMuonSelection(cand1, cand2, beamSpot, bFieldHandle))
0375             continue;
0376
0377           // Add this pair
0378           n++;
0379           LogDebug("HLTMuonDimuonL3Filter")
0380               << " Track1 passing filter: pt= " << cand1->pt() << ", eta: " << cand1->eta();
0381           LogDebug("HLTMuonDimuonL3Filter")
0382               << " Track2 passing filter: pt= " << cand2->pt() << ", eta: " << cand2->eta();
0383
0384           bool i1done = false;
0385           bool i2done = false;
0386           vector<RecoChargedCandidateRef> vref;
0387           filterproduct.getObjects(TriggerMuon, vref);
0388           for (auto& i : vref) {
0389             RecoChargedCandidateRef candref = RecoChargedCandidateRef(i);
0390             TrackRef tktmp = candref->get<TrackRef>();
0391             if (tktmp == tk1)
0392               i1done = true;
0393             else if (tktmp == tk2)
0394               i2done = true;  //why is this an elif?
0395             if (i1done && i2done)
0396               break;
0397           }
0398           if (!i1done)
0399             filterproduct.addObject(TriggerMuon, cand1);
0400           if (!i2done)
0401             filterproduct.addObject(TriggerMuon, cand2);
0402
0403           //break anyway since a L3 track pair has been found matching the criteria
0404           thisL3Index1isDone = true;
0405           atLeastOnePair = true;
0406           break;
0407         }  //loop on the track of the second L2
0408         //break the loop if fast accept.
0409         if (atLeastOnePair && fast_Accept_)
0410           break;
0411       }  //loop on the second L2
0412       //break the loop if fast accept.
0413       if (atLeastOnePair && fast_Accept_)
0414         break;
0415       if (thisL3Index1isDone)
0416         break;
0417
0418       //Loop over L3FromL1 collection see if we get a pair that way
0419       auto MuonToL3s_it1 = MuonToL3s.begin();
0420       auto MuonToL3s_end = MuonToL3s.end();
0421       for (; MuonToL3s_it1 != MuonToL3s_end; ++MuonToL3s_it1) {
0422         const RecoChargedCandidateRef& cand2 = MuonToL3s_it1->second;
0423         if (!applyMuonSelection(cand2, beamSpot))
0424           continue;
0425         TrackRef tk2 = cand2->get<TrackRef>();
0426
0427         // Run dimuon selection:
0428         if (!applyDiMuonSelection(cand1, cand2, beamSpot, bFieldHandle))
0429           continue;
0430         n++;
0431         LogDebug("HLTMuonDimuonL3Filter")
0432             << " L3FromL2 Track1 passing filter: pt= " << cand1->pt() << ", eta: " << cand1->eta();
0433         LogDebug("HLTMuonDimuonL3Filter")
0434             << " L3FromL1 Track2 passing filter: pt= " << cand2->pt() << ", eta: " << cand2->eta();
0435
0436         bool i1done = false;
0437         bool i2done = false;
0438         vector<RecoChargedCandidateRef> vref;
0439         filterproduct.getObjects(TriggerMuon, vref);
0440         for (auto& i : vref) {
0441           RecoChargedCandidateRef candref = RecoChargedCandidateRef(i);
0442           TrackRef tktmp = candref->get<TrackRef>();
0443           if (tktmp == tk1)
0444             i1done = true;
0445           else if (tktmp == tk2)
0446             i2done = true;  //why is this an elif?
0447           if (i1done && i2done)
0448             break;
0449         }
0450         if (!i1done)
0451           filterproduct.addObject(TriggerMuon, cand1);
0452         if (!i2done)
0453           filterproduct.addObject(TriggerMuon, cand2);
0454
0455         //break anyway since a L3 track pair has been found matching the criteria
0456         thisL3Index1isDone = true;
0457         atLeastOnePair = true;
0458         break;
0459       }  //L3FromL1 loop
0460       //break the loop if fast accept.
0461       if (atLeastOnePair && fast_Accept_)
0462         break;
0463       if (thisL3Index1isDone)
0464         break;
0465
0466     }  //loop on tracks for first L2
0467     //break the loop if fast accept.
0468     if (atLeastOnePair && fast_Accept_)
0469       break;
0470   }  //loop on the first L2
0471
0472   // now loop on 1st L3 from L1
0473   auto MuonToL3s_it1 = MuonToL3s.begin();
0474   auto MuonToL3s_end = MuonToL3s.end();
0475   for (; MuonToL3s_it1 != MuonToL3s_end; ++MuonToL3s_it1) {
0476     bool thisL3Index1isDone = false;
0477     const RecoChargedCandidateRef& cand1 = MuonToL3s_it1->second;
0478     if (!applyMuonSelection(cand1, beamSpot))
0479       continue;
0480     TrackRef tk1 = cand1->get<TrackRef>();
0481
0482     // Loop on 2nd L3 from L1
0483     auto MuonToL3s_it2 = MuonToL3s_it1;
0484     for (; MuonToL3s_it2 != MuonToL3s_end; ++MuonToL3s_it2) {
0485       const RecoChargedCandidateRef& cand2 = MuonToL3s_it2->second;
0486       if (!applyMuonSelection(cand2, beamSpot))
0487         continue;
0488       TrackRef tk2 = cand2->get<TrackRef>();
0489
0490       // Run dimuon selection:
0491       if (!applyDiMuonSelection(cand1, cand2, beamSpot, bFieldHandle))
0492         continue;
0493
0494       n++;
0495       LogDebug("HLTMuonDimuonL3Filter") << " L3FromL1 Track1 passing filter: pt= " << cand1->pt()
0496                                         << ", eta: " << cand1->eta();
0497       LogDebug("HLTMuonDimuonL3Filter") << " L3FromL1 Track2 passing filter: pt= " << cand2->pt()
0498                                         << ", eta: " << cand2->eta();
0499
0500       bool i1done = false;
0501       bool i2done = false;
0502       vector<RecoChargedCandidateRef> vref;
0503       filterproduct.getObjects(TriggerMuon, vref);
0504       for (auto& i : vref) {
0505         RecoChargedCandidateRef candref = RecoChargedCandidateRef(i);
0506         TrackRef tktmp = candref->get<TrackRef>();
0507         if (tktmp == tk1)
0508           i1done = true;
0509         else if (tktmp == tk2)
0510           i2done = true;  //why is this an elif?
0511         if (i1done && i2done)
0512           break;
0513       }
0514       if (!i1done)
0515         filterproduct.addObject(TriggerMuon, cand1);
0516       if (!i2done)
0517         filterproduct.addObject(TriggerMuon, cand2);
0518
0519       //break anyway since a L3 track pair has been found matching the criteria
0520       thisL3Index1isDone = true;
0521       atLeastOnePair = true;
0522       break;
0523     }  //loop on 2nd muon
0524
0525     //break the loop if fast accept
0526     if (atLeastOnePair && fast_Accept_)
0527       break;
0528     if (thisL3Index1isDone)
0529       break;
0530   }  //loop on 1st muon
0531
0532   // filter decision
0533   const bool accept(n >= min_N_);
0534
0535   LogDebug("HLTMuonDimuonL3Filter") << " >>>>> Result of HLTMuonDimuonL3Filter is " << accept
0536                                     << ", number of muon pairs passing thresholds= " << n;
0537
0538   return accept;
0539 }
0540
0541 bool HLTMuonDimuonL3Filter::triggeredByLevel2(TrackRef const& staTrack, vector<RecoChargedCandidateRef> const& vcands) {
0542   bool ok = false;
0543   for (auto const& vcand : vcands) {
0544     if (vcand->get<TrackRef>() == staTrack) {
0545       ok = true;
0546       LogDebug("HLTMuonL3PreFilter") << "The L2 track triggered";
0547       break;
0548     }
0549   }
0550   return ok;
0551 }
0552
0553 bool HLTMuonDimuonL3Filter::applyMuonSelection(const RecoChargedCandidateRef& cand, const BeamSpot& beamSpot) const {
0554   // eta cut
0555   if (std::abs(cand->eta()) > max_Eta_)
0556     return false;
0557
0558   // cut on number of hits
0559   TrackRef tk = cand->track();
0560   if (tk->numberOfValidHits() < min_Nhits_)
0561     return false;
0562
0563   //dr cut
0564   if (std::abs((-(cand->vx() - beamSpot.x0()) * cand->py() + (cand->vy() - beamSpot.y0()) * cand->px()) / cand->pt()) >
0565       max_Dr_)
0566     return false;
0567
0568   //dz cut
0569   if (std::abs((cand->vz() - beamSpot.z0()) -
0570                ((cand->vx() - beamSpot.x0()) * cand->px() + (cand->vy() - beamSpot.y0()) * cand->py()) / cand->pt() *
0571                    cand->pz() / cand->pt()) > max_Dz_)
0572     return false;
0573
0574   return true;
0575 }
0576
0577 bool HLTMuonDimuonL3Filter::applyDiMuonSelection(const RecoChargedCandidateRef& cand1,
0578                                                  const RecoChargedCandidateRef& cand2,
0579                                                  const BeamSpot& beamSpot,
0580                                                  const ESHandle<MagneticField>& bFieldHandle) const {
0581   // Opposite Charge
0582   if (chargeOpt_ < 0 and (cand1->charge() * cand2->charge() > 0))
0583     return false;
0584   else if (chargeOpt_ > 0 and (cand1->charge() * cand2->charge() < 0))
0585     return false;
0586
0587   // Acoplanarity
0588   double acop = std::abs(cand1->phi() - cand2->phi());
0589   if (acop > M_PI)
0590     acop = 2 * M_PI - acop;
0591   acop = M_PI - acop;
0592   LogDebug("HLTMuonDimuonL3Filter") << " ... 1-2 acop= " << acop;
0593   if (acop < min_Acop_)
0594     return false;
0595   if (acop > max_Acop_)
0596     return false;
0597
0598   // Pt balance
0599   double ptbalance = std::abs(cand1->pt() - cand2->pt());
0600   if (ptbalance < min_PtBalance_)
0601     return false;
0602   if (ptbalance > max_PtBalance_)
0603     return false;
0604
0605   // Combined dimuon syste
0606   double e1, e2;
0607   Particle::LorentzVector p, p1, p2;
0608   e1 = sqrt(cand1->momentum().Mag2() + MuMass2_);
0609   e2 = sqrt(cand2->momentum().Mag2() + MuMass2_);
0610   p1 = Particle::LorentzVector(cand1->px(), cand1->py(), cand1->pz(), e1);
0611   p2 = Particle::LorentzVector(cand2->px(), cand2->py(), cand2->pz(), e2);
0612   p = p1 + p2;
0613
0614   double pt12 = p.pt();
0615   LogDebug("HLTMuonDimuonL3Filter") << " ... 1-2 pt12= " << pt12;
0616
0617   // Angle between the muons
0618   if (applyMinDiMuonDeltaR2Cut_ and reco::deltaR2(p1, p2) < min_DiMuonDeltaR2_)
0619     return false;
0620
0621   double ptLx1 = cand1->pt();
0622   double ptLx2 = cand2->pt();
0623   double invmass = abs(p.mass());
0624   // if (invmass>0) invmass = sqrt(invmass); else invmass = 0;
0625   LogDebug("HLTMuonDimuonL3Filter") << " ... 1-2 invmass= " << invmass;
0626   bool proceed = false;
0627   for (unsigned int iv = 0; iv < min_InvMass_.size(); iv++) {
0628     if (invmass < min_InvMass_[iv])
0629       continue;
0630     if (invmass > max_InvMass_[iv])
0631       continue;
0632     if (ptLx1 > ptLx2) {
0633       if (ptLx1 < min_PtMax_[iv])
0634         continue;
0635       if (ptLx2 < min_PtMin_[iv])
0636         continue;
0637       if (ptLx2 > max_PtMin_[iv])
0638         continue;
0639     } else {
0640       if (ptLx2 < min_PtMax_[iv])
0641         continue;
0642       if (ptLx1 < min_PtMin_[iv])
0643         continue;
0644       if (ptLx1 > max_PtMin_[iv])
0645         continue;
0646     }
0647     if (pt12 < min_PtPair_[iv])
0648       continue;
0649     if (pt12 > max_PtPair_[iv])
0650       continue;
0651     proceed = true;
0652     break;
0653   }
0654   if (!proceed)
0655     return false;
0656
0657   // Delta Z between the two muons
0658   //double DeltaZMuMu = std::abs(tk2->dz(beamSpot.position())-tk1->dz(beamSpot.position()));
0659   //if ( DeltaZMuMu > max_DzMuMu_) return false;
0660
0661   // DCA between the two muons
0662   TrackRef tk1 = cand1->track();
0663   TrackRef tk2 = cand2->track();
0664   TransientTrack mu1TT(*tk1, &(*bFieldHandle));
0665   TransientTrack mu2TT(*tk2, &(*bFieldHandle));
0666   TrajectoryStateClosestToPoint mu1TS = mu1TT.impactPointTSCP();
0667   TrajectoryStateClosestToPoint mu2TS = mu2TT.impactPointTSCP();
0668   if (mu1TS.isValid() && mu2TS.isValid()) {
0669     ClosestApproachInRPhi cApp;
0670     cApp.calculate(mu1TS.theState(), mu2TS.theState());
0671     if (!cApp.status() || cApp.distance() > max_DCAMuMu_)
0672       return false;
0673   }
0674
0675   // Max dimuon |rapidity|
0676   double rapidity = std::abs(p.Rapidity());
0677   if (rapidity > max_YPair_)
0678     return false;
0679
0680   // if cutting on cowboys reject muons that bend towards each other
0681   if (cutCowboys_ && (cand1->charge() * deltaPhi(cand1->phi(), cand2->phi()) > 0.))
0682     return false;
0683   return true;
0684 }
0685
0686 // declare this class as a framework plugin
0687 #include "FWCore/Framework/interface/MakerMacros.h"
0688 DEFINE_FWK_MODULE(HLTMuonDimuonL3Filter);