Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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