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