Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:09:38 

0001 /** \class HLTMuonDimuonL2Filter
0002  *
0003  * See header file for documentation
0004  *
0005  *  \author J. Alcaraz, P. Garcia, M.Vander Donckt
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 // constructors and destructor
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 // member functions
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 // ------------ method called to produce the data  ------------
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   // All HLT filters must create and fill an HLT filter object,
0122   // recording any reconstructed physics objects satisfying (or not)
0123   // this HLT filter, and place it in the Event.
0124 
0125   // Ref to Candidate object to be recorded in filter object
0126   RecoChargedCandidateRef ref1;
0127   RecoChargedCandidateRef ref2;
0128 
0129   // get hold of trks
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   // get the L2 to L1 map object for this event
0141   HLTMuonL2ToL1Map mapL2ToL1(previousCandToken_, seedMapToken_, iEvent);
0142 
0143   // look at all mucands,  check cuts and add to filter object
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     // eta cut
0153     LogDebug("HLTMuonDimuonL2Filter") << " 1st muon in loop: q*pt= " << tk1->charge() * tk1->pt()
0154                                       << ", eta= " << tk1->eta() << ", hits= " << tk1->numberOfValidHits();
0155     // find the L1 Particle corresponding to the L2 Track
0156     if (!mapL2ToL1.isTriggeredByL1(tk1))
0157       continue;
0158 
0159     if (fabs(tk1->eta()) > max_Eta_)
0160       continue;
0161 
0162     // cut on number of hits
0163     if (tk1->numberOfValidHits() < min_Nhits_)
0164       continue;
0165 
0166     // number of stations
0167     if (tk1->hitPattern().muonStationsWithAnyHits() < min_Nstations_)
0168       continue;
0169 
0170     // number of chambers
0171     if (tk1->hitPattern().dtStationsWithAnyHits() + tk1->hitPattern().cscStationsWithAnyHits() < min_Nchambers_)
0172       continue;
0173 
0174     //dr cut
0175     //if (fabs(tk1->d0())>max_Dr_) continue;
0176     if (fabs(tk1->dxy(beamSpot.position())) > max_Dr_)
0177       continue;
0178 
0179     //dz cut
0180     if (fabs(tk1->dz()) > max_Dz_)
0181       continue;
0182 
0183     // Pt threshold cut
0184     double pt1 = tk1->pt();
0185     double err1 = tk1->error(0);
0186     double abspar1 = fabs(tk1->parameter(0));
0187     double ptLx1 = pt1;
0188     // convert 50% efficiency threshold to 90% efficiency threshold
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       // eta cut
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       // cut on number of hits
0209       if (tk2->numberOfValidHits() < min_Nhits_)
0210         continue;
0211 
0212       // number of stations
0213       if (tk2->hitPattern().muonStationsWithAnyHits() < min_Nstations_)
0214         continue;
0215 
0216       // number of chambers
0217       if (tk2->hitPattern().dtStationsWithAnyHits() + tk2->hitPattern().cscStationsWithAnyHits() < min_Nchambers_)
0218         continue;
0219 
0220       //dr cut
0221       //if (fabs(tk2->d0())>max_Dr_) continue;
0222       if (fabs(tk2->dxy(beamSpot.position())) > max_Dr_)
0223         continue;
0224 
0225       //dz cut
0226       if (fabs(tk2->dz()) > max_Dz_)
0227         continue;
0228 
0229       // Pt threshold cut
0230       double pt2 = tk2->pt();
0231       double err2 = tk2->error(0);
0232       double abspar2 = fabs(tk2->parameter(0));
0233       double ptLx2 = pt2;
0234       // convert 50% efficiency threshold to 90% efficiency threshold
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       // Acoplanarity
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       // 3D angle
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       // Pt balance
0280       double ptbalance = fabs(tk1->pt() - tk2->pt());
0281       if (ptbalance < min_PtBalance_)
0282         continue;
0283       if (ptbalance > max_PtBalance_)
0284         continue;
0285 
0286       // Combined dimuon system
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       // if (invmass>0) invmass = sqrt(invmass); else invmass = 0;
0300       LogDebug("HLTMuonDimuonL2Filter") << " ... 1-2 invmass= " << invmass;
0301       if (invmass < min_InvMass_)
0302         continue;
0303       if (invmass > max_InvMass_)
0304         continue;
0305 
0306       // Add this pair
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   // filter decision
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 // declare this class as a framework plugin
0352 #include "FWCore/Framework/interface/MakerMacros.h"
0353 DEFINE_FWK_MODULE(HLTMuonDimuonL2Filter);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!