File indexing completed on 2023-03-17 11:09:38
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "DataFormats/Common/interface/Handle.h"
0010
0011 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0012 #include "DataFormats/HLTReco/interface/TriggerRefsCollections.h"
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0018 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0019 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0020 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeed.h"
0021 #include "DataFormats/MuonSeed/interface/L2MuonTrajectorySeedCollection.h"
0022 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
0023 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
0024
0025 #include "HLTMuonDimuonL2Filter.h"
0026 #include "HLTMuonL2ToL1Map.h"
0027 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0028
0029 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0030
0031 using namespace edm;
0032 using namespace std;
0033 using namespace reco;
0034 using namespace trigger;
0035 using namespace l1extra;
0036
0037
0038
0039 HLTMuonDimuonL2Filter::HLTMuonDimuonL2Filter(const edm::ParameterSet& iConfig)
0040 : HLTFilter(iConfig),
0041 beamspotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0042 beamspotToken_(consumes<reco::BeamSpot>(beamspotTag_)),
0043 candTag_(iConfig.getParameter<edm::InputTag>("CandTag")),
0044 candToken_(consumes<reco::RecoChargedCandidateCollection>(candTag_)),
0045 previousCandTag_(iConfig.getParameter<InputTag>("PreviousCandTag")),
0046 previousCandToken_(consumes<trigger::TriggerFilterObjectWithRefs>(previousCandTag_)),
0047 seedMapTag_(iConfig.getParameter<InputTag>("SeedMapTag")),
0048 seedMapToken_(consumes<SeedMap>(seedMapTag_)),
0049 fast_Accept_(iConfig.getParameter<bool>("FastAccept")),
0050 max_Eta_(iConfig.getParameter<double>("MaxEta")),
0051 min_Nhits_(iConfig.getParameter<int>("MinNhits")),
0052 min_Nstations_(iConfig.getParameter<int>("MinNstations")),
0053 min_Nchambers_(iConfig.getParameter<int>("MinNchambers")),
0054 max_Dr_(iConfig.getParameter<double>("MaxDr")),
0055 max_Dz_(iConfig.getParameter<double>("MaxDz")),
0056 chargeOpt_(iConfig.getParameter<int>("ChargeOpt")),
0057 min_PtPair_(iConfig.getParameter<double>("MinPtPair")),
0058 min_PtMax_(iConfig.getParameter<double>("MinPtMax")),
0059 min_PtMin_(iConfig.getParameter<double>("MinPtMin")),
0060 min_InvMass_(iConfig.getParameter<double>("MinInvMass")),
0061 max_InvMass_(iConfig.getParameter<double>("MaxInvMass")),
0062 min_Acop_(iConfig.getParameter<double>("MinAcop")),
0063 max_Acop_(iConfig.getParameter<double>("MaxAcop")),
0064 min_Angle_(iConfig.getParameter<double>("MinAngle")),
0065 max_Angle_(iConfig.getParameter<double>("MaxAngle")),
0066 min_PtBalance_(iConfig.getParameter<double>("MinPtBalance")),
0067 max_PtBalance_(iConfig.getParameter<double>("MaxPtBalance")),
0068 nsigma_Pt_(iConfig.getParameter<double>("NSigmaPt")) {
0069 LogDebug("HLTMuonDimuonL2Filter")
0070 << " CandTag/MinN/MaxEta/MinNhits/MinNstations/MinNchambers/MaxDr/MaxDz/MinPt1/MinPt2/MinInvMass/MaxInvMass/"
0071 "MinAcop/MaxAcop/MinAngle/MaxAngle/MinPtBalance/MaxPtBalance/NSigmaPt : "
0072 << candTag_.encode() << " " << fast_Accept_ << " " << max_Eta_ << " " << min_Nhits_ << " " << min_Nstations_
0073 << " " << min_Nchambers_ << " " << max_Dr_ << " " << max_Dz_ << " " << chargeOpt_ << " " << min_PtPair_ << " "
0074 << min_PtMax_ << " " << min_PtMin_ << " " << min_InvMass_ << " " << max_InvMass_ << " " << min_Acop_ << " "
0075 << max_Acop_ << " " << min_Angle_ << " " << max_Angle_ << " " << min_PtBalance_ << " " << max_PtBalance_ << " "
0076 << nsigma_Pt_;
0077 }
0078
0079 HLTMuonDimuonL2Filter::~HLTMuonDimuonL2Filter() = default;
0080
0081
0082
0083
0084
0085 void HLTMuonDimuonL2Filter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0086 edm::ParameterSetDescription desc;
0087 makeHLTFilterDescription(desc);
0088 desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOfflineBeamSpot"));
0089 desc.add<edm::InputTag>("CandTag", edm::InputTag("hltL2MuonCandidates"));
0090 desc.add<edm::InputTag>("PreviousCandTag", edm::InputTag(""));
0091 desc.add<edm::InputTag>("SeedMapTag", edm::InputTag("hltL2Muons"));
0092 desc.add<bool>("FastAccept", false);
0093 desc.add<double>("MaxEta", 2.5);
0094 desc.add<int>("MinNhits", 0);
0095 desc.add<int>("MinNstations", 0);
0096 desc.add<int>("MinNchambers", 2);
0097 desc.add<double>("MaxDr", 100.0);
0098 desc.add<double>("MaxDz", 9999.0);
0099 desc.add<int>("ChargeOpt", 0);
0100 desc.add<double>("MinPtPair", 0.0);
0101 desc.add<double>("MinPtMax", 3.0);
0102 desc.add<double>("MinPtMin", 3.0);
0103 desc.add<double>("MinInvMass", 1.6);
0104 desc.add<double>("MaxInvMass", 5.6);
0105 desc.add<double>("MinAcop", -1.0);
0106 desc.add<double>("MaxAcop", 3.15);
0107 desc.add<double>("MinAngle", -999.0);
0108 desc.add<double>("MaxAngle", 2.5);
0109 desc.add<double>("MinPtBalance", -1.0);
0110 desc.add<double>("MaxPtBalance", 999999.0);
0111 desc.add<double>("NSigmaPt", 0.0);
0112 descriptions.add("hltMuonDimuonL2Filter", desc);
0113 }
0114
0115
0116 bool HLTMuonDimuonL2Filter::hltFilter(edm::Event& iEvent,
0117 const edm::EventSetup& iSetup,
0118 trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0119 double const MuMass = 0.106;
0120 double const MuMass2 = MuMass * MuMass;
0121
0122
0123
0124
0125
0126 RecoChargedCandidateRef ref1;
0127 RecoChargedCandidateRef ref2;
0128
0129
0130 Handle<RecoChargedCandidateCollection> mucands;
0131 iEvent.getByToken(candToken_, mucands);
0132 if (saveTags())
0133 filterproduct.addCollectionTag(candTag_);
0134
0135 BeamSpot beamSpot;
0136 Handle<BeamSpot> recoBeamSpotHandle;
0137 iEvent.getByToken(beamspotToken_, recoBeamSpotHandle);
0138 beamSpot = *recoBeamSpotHandle;
0139
0140
0141 HLTMuonL2ToL1Map mapL2ToL1(previousCandToken_, seedMapToken_, iEvent);
0142
0143
0144 int n = 0;
0145 double e1, e2;
0146 Particle::LorentzVector p, p1, p2;
0147
0148 RecoChargedCandidateCollection::const_iterator cand1;
0149 RecoChargedCandidateCollection::const_iterator cand2;
0150 for (cand1 = mucands->begin(); cand1 != mucands->end(); cand1++) {
0151 TrackRef tk1 = cand1->get<TrackRef>();
0152
0153 LogDebug("HLTMuonDimuonL2Filter") << " 1st muon in loop: q*pt= " << tk1->charge() * tk1->pt()
0154 << ", eta= " << tk1->eta() << ", hits= " << tk1->numberOfValidHits();
0155
0156 if (!mapL2ToL1.isTriggeredByL1(tk1))
0157 continue;
0158
0159 if (fabs(tk1->eta()) > max_Eta_)
0160 continue;
0161
0162
0163 if (tk1->numberOfValidHits() < min_Nhits_)
0164 continue;
0165
0166
0167 if (tk1->hitPattern().muonStationsWithAnyHits() < min_Nstations_)
0168 continue;
0169
0170
0171 if (tk1->hitPattern().dtStationsWithAnyHits() + tk1->hitPattern().cscStationsWithAnyHits() < min_Nchambers_)
0172 continue;
0173
0174
0175
0176 if (fabs(tk1->dxy(beamSpot.position())) > max_Dr_)
0177 continue;
0178
0179
0180 if (fabs(tk1->dz()) > max_Dz_)
0181 continue;
0182
0183
0184 double pt1 = tk1->pt();
0185 double err1 = tk1->error(0);
0186 double abspar1 = fabs(tk1->parameter(0));
0187 double ptLx1 = pt1;
0188
0189 if (abspar1 > 0)
0190 ptLx1 += nsigma_Pt_ * err1 / abspar1 * pt1;
0191 LogDebug("HLTMuonDimuonL2Filter") << " ... 1st muon in loop, pt1= " << pt1 << ", ptLx1= " << ptLx1;
0192
0193 cand2 = cand1;
0194 cand2++;
0195 for (; cand2 != mucands->end(); cand2++) {
0196 TrackRef tk2 = cand2->get<TrackRef>();
0197
0198
0199 LogDebug("HLTMuonDimuonL2Filter") << " 2nd muon in loop: q*pt= " << tk2->charge() * tk2->pt()
0200 << ", eta= " << tk2->eta() << ", hits= " << tk2->numberOfValidHits()
0201 << ", d0= " << tk2->d0();
0202 if (!mapL2ToL1.isTriggeredByL1(tk2))
0203 continue;
0204
0205 if (fabs(tk2->eta()) > max_Eta_)
0206 continue;
0207
0208
0209 if (tk2->numberOfValidHits() < min_Nhits_)
0210 continue;
0211
0212
0213 if (tk2->hitPattern().muonStationsWithAnyHits() < min_Nstations_)
0214 continue;
0215
0216
0217 if (tk2->hitPattern().dtStationsWithAnyHits() + tk2->hitPattern().cscStationsWithAnyHits() < min_Nchambers_)
0218 continue;
0219
0220
0221
0222 if (fabs(tk2->dxy(beamSpot.position())) > max_Dr_)
0223 continue;
0224
0225
0226 if (fabs(tk2->dz()) > max_Dz_)
0227 continue;
0228
0229
0230 double pt2 = tk2->pt();
0231 double err2 = tk2->error(0);
0232 double abspar2 = fabs(tk2->parameter(0));
0233 double ptLx2 = pt2;
0234
0235 if (abspar2 > 0)
0236 ptLx2 += nsigma_Pt_ * err2 / abspar2 * pt2;
0237 LogDebug("HLTMuonDimuonL2Filter") << " ... 2nd muon in loop, pt2= " << pt2 << ", ptLx2= " << ptLx2;
0238
0239 if (ptLx1 > ptLx2) {
0240 if (ptLx1 < min_PtMax_)
0241 continue;
0242 if (ptLx2 < min_PtMin_)
0243 continue;
0244 } else {
0245 if (ptLx2 < min_PtMax_)
0246 continue;
0247 if (ptLx1 < min_PtMin_)
0248 continue;
0249 }
0250
0251 if (chargeOpt_ < 0) {
0252 if (tk1->charge() * tk2->charge() > 0)
0253 continue;
0254 } else if (chargeOpt_ > 0) {
0255 if (tk1->charge() * tk2->charge() < 0)
0256 continue;
0257 }
0258
0259
0260 double acop = fabs(tk1->phi() - tk2->phi());
0261 if (acop > M_PI)
0262 acop = 2 * M_PI - acop;
0263 acop = M_PI - acop;
0264 LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 acop= " << acop;
0265 if (acop < min_Acop_)
0266 continue;
0267 if (acop > max_Acop_)
0268 continue;
0269
0270
0271 double angle =
0272 acos((tk1->px() * tk2->px() + tk1->py() * tk2->py() + tk1->pz() * tk2->pz()) / (tk1->p() * tk2->p()));
0273 LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 angle= " << angle;
0274 if (angle < min_Angle_)
0275 continue;
0276 if (angle > max_Angle_)
0277 continue;
0278
0279
0280 double ptbalance = fabs(tk1->pt() - tk2->pt());
0281 if (ptbalance < min_PtBalance_)
0282 continue;
0283 if (ptbalance > max_PtBalance_)
0284 continue;
0285
0286
0287 e1 = sqrt(tk1->momentum().Mag2() + MuMass2);
0288 e2 = sqrt(tk2->momentum().Mag2() + MuMass2);
0289 p1 = Particle::LorentzVector(tk1->px(), tk1->py(), tk1->pz(), e1);
0290 p2 = Particle::LorentzVector(tk2->px(), tk2->py(), tk2->pz(), e2);
0291 p = p1 + p2;
0292
0293 double pt12 = p.pt();
0294 LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 pt12= " << pt12;
0295 if (pt12 < min_PtPair_)
0296 continue;
0297
0298 double invmass = abs(p.mass());
0299
0300 LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 invmass= " << invmass;
0301 if (invmass < min_InvMass_)
0302 continue;
0303 if (invmass > max_InvMass_)
0304 continue;
0305
0306
0307 n++;
0308 LogDebug("HLTMuonDimuonL2Filter") << " Track1 passing filter: pt= " << tk1->pt() << ", eta: " << tk1->eta();
0309 LogDebug("HLTMuonDimuonL2Filter") << " Track2 passing filter: pt= " << tk2->pt() << ", eta: " << tk2->eta();
0310 LogDebug("HLTMuonDimuonL2Filter") << " Invmass= " << invmass;
0311
0312 bool i1done = false;
0313 bool i2done = false;
0314 vector<RecoChargedCandidateRef> vref;
0315 filterproduct.getObjects(TriggerMuon, vref);
0316 for (auto& i : vref) {
0317 RecoChargedCandidateRef candref = RecoChargedCandidateRef(i);
0318 TrackRef tktmp = candref->get<TrackRef>();
0319 if (tktmp == tk1) {
0320 i1done = true;
0321 } else if (tktmp == tk2) {
0322 i2done = true;
0323 }
0324 if (i1done && i2done)
0325 break;
0326 }
0327
0328 if (!i1done) {
0329 ref1 = RecoChargedCandidateRef(Ref<RecoChargedCandidateCollection>(mucands, distance(mucands->begin(), cand1)));
0330 filterproduct.addObject(TriggerMuon, ref1);
0331 }
0332 if (!i2done) {
0333 ref2 = RecoChargedCandidateRef(Ref<RecoChargedCandidateCollection>(mucands, distance(mucands->begin(), cand2)));
0334 filterproduct.addObject(TriggerMuon, ref2);
0335 }
0336
0337 if (fast_Accept_)
0338 break;
0339 }
0340 }
0341
0342
0343 const bool accept(n >= 1);
0344
0345 LogDebug("HLTMuonDimuonL2Filter") << " >>>>> Result of HLTMuonDimuonL2Filter is " << accept
0346 << ", number of muon pairs passing thresholds= " << n;
0347
0348 return accept;
0349 }
0350
0351
0352 #include "FWCore/Framework/interface/MakerMacros.h"
0353 DEFINE_FWK_MODULE(HLTMuonDimuonL2Filter);