Skip to content

Commit

Permalink
Address PR comments
Browse files Browse the repository at this point in the history
- Introduce more named constants for later maintainability
- Make the maxChi2 in the estimator customizable
- Use vdt trig functions and precompute all values when possible
- Fix code style: remove not-needed member variable, use DTChamberId
  constants instead of magic numbers
  • Loading branch information
Parsifal-2045 committed Dec 17, 2024
1 parent 4762785 commit 62c4e64
Show file tree
Hide file tree
Showing 6 changed files with 97 additions and 76 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -22,24 +22,28 @@
)

phase2HltL2MuonSeedsFromL1TkMuon = cms.EDProducer('Phase2L2MuonSeedCreator',
InputObjects = cms.InputTag('l1tTkMuonsGmt'),
CSCRecSegmentLabel = cms.InputTag('hltCscSegments'),
DTRecSegmentLabel = cms.InputTag('hltDt4DSegments'),
MinPL1Tk = cms.double(3.5),
MaxPL1Tk = cms.double(200),
inputObjects = cms.InputTag('l1tTkMuonsGmt'),
cscRecSegmentLabel = cms.InputTag('hltCscSegments'),
dtRecSegmentLabel = cms.InputTag('hltDt4DSegments'),
minPL1Tk = cms.double(3.5),
maxPL1Tk = cms.double(200),
stubMatchDPhi = cms.double(0.05),
stubMatchDTheta = cms.double(0.1),
extrapolationWindowClose = cms.double(0.2),
extrapolationWindowFar = cms.double(0.1),
maximumEtaBarrel = cms.double(0.7),
maximumEtaOverlap = cms.double(1.3),
Propagator = cms.string('SteppingHelixPropagatorAny'),
ServiceParameters = cms.PSet(
propagator = cms.string('SteppingHelixPropagatorAny'),
serviceParameters = cms.PSet(
Propagators = cms.untracked.vstring('SteppingHelixPropagatorAny'),
RPCLayers = cms.bool(True),
UseMuonNavigation = cms.untracked.bool(True)
)
),
estimatorMaxChi2 = cms.double(1000.0)
)

from Configuration.ProcessModifiers.phase2L2AndL3Muons_cff import phase2L2AndL3Muons
phase2L2AndL3Muons.toReplaceWith(hltL2MuonSeedsFromL1TkMuon, phase2HltL2MuonSeedsFromL1TkMuon)
phase2L2AndL3Muons.toReplaceWith(
hltL2MuonSeedsFromL1TkMuon,
phase2HltL2MuonSeedsFromL1TkMuon
)
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
import FWCore.ParameterSet.Config as cms

hltPhase2L3MuonFilter = cms.EDProducer("phase2HLTMuonSelectorForL3",
hltPhase2L3MuonFilter = cms.EDProducer("Phase2HLTMuonSelectorForL3",
l1TkMuons = cms.InputTag("l1tTkMuonsGmt"),
l2MuonsUpdVtx = cms.InputTag("hltL2MuonsFromL1TkMuon:UpdatedAtVtx"),
l3Tracks = cms.InputTag("hltIter2Phase2L3FromL1TkMuonMerged"),
Expand Down
116 changes: 67 additions & 49 deletions RecoMuon/L2MuonSeedGenerator/src/Phase2L2MuonSeedCreator.cc
Original file line number Diff line number Diff line change
Expand Up @@ -30,56 +30,59 @@
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "CLHEP/Vector/ThreeVector.h"
#include <vdt/atan.h>
#include <vdt/exp.h>
#include <vdt/sincos.h>

#include <vector>

// Constructor
Phase2L2MuonSeedCreator::Phase2L2MuonSeedCreator(const edm::ParameterSet& pset)
: l1TkMuCollToken_{consumes(pset.getParameter<edm::InputTag>("InputObjects"))},
cscSegmentCollToken_{consumes(pset.getParameter<edm::InputTag>("CSCRecSegmentLabel"))},
dtSegmentCollToken_{consumes(pset.getParameter<edm::InputTag>("DTRecSegmentLabel"))},
: l1TkMuCollToken_{consumes(pset.getParameter<edm::InputTag>("inputObjects"))},
cscSegmentCollToken_{consumes(pset.getParameter<edm::InputTag>("cscRecSegmentLabel"))},
dtSegmentCollToken_{consumes(pset.getParameter<edm::InputTag>("dtRecSegmentLabel"))},
cscGeometryToken_{esConsumes<CSCGeometry, MuonGeometryRecord>()},
dtGeometryToken_{esConsumes<DTGeometry, MuonGeometryRecord>()},
magneticFieldToken_{esConsumes<MagneticField, IdealMagneticFieldRecord>()},
minMomentum_{pset.getParameter<double>("MinPL1Tk")},
maxMomentum_{pset.getParameter<double>("MaxPL1Tk")},
minMomentum_{pset.getParameter<double>("minPL1Tk")},
maxMomentum_{pset.getParameter<double>("maxPL1Tk")},
matchingPhiWindow_{pset.getParameter<double>("stubMatchDPhi")},
matchingThetaWindow_{pset.getParameter<double>("stubMatchDTheta")},
extrapolationDeltaPhiClose_{pset.getParameter<double>("extrapolationWindowClose")},
extrapolationDeltaPhiFar_{pset.getParameter<double>("extrapolationWindowFar")},
maxEtaBarrel_{pset.getParameter<double>("maximumEtaBarrel")},
maxEtaOverlap_{pset.getParameter<double>("maximumEtaOverlap")},
propagatorName_{pset.getParameter<string>("Propagator")} {
propagatorName_{pset.getParameter<string>("propagator")} {
// Service parameters
edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters");
edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("serviceParameters");
// Services
service_ = std::make_unique<MuonServiceProxy>(serviceParameters, consumesCollector());
estimator_ = std::make_unique<Chi2MeasurementEstimator>(10000.);
estimator_ = std::make_unique<Chi2MeasurementEstimator>(pset.getParameter<double>("estimatorMaxChi2"));
produces<L2MuonTrajectorySeedCollection>();
}

void Phase2L2MuonSeedCreator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("InputObjects", edm::InputTag("l1tTkMuonsGmt"));
desc.add<edm::InputTag>("CSCRecSegmentLabel", edm::InputTag("hltCscSegments"));
desc.add<edm::InputTag>("DTRecSegmentLabel", edm::InputTag("hltDt4DSegments"));
desc.add<double>("MinPL1Tk", 3.5);
desc.add<double>("MaxPL1Tk", 200);
desc.add<edm::InputTag>("inputObjects", edm::InputTag("l1tTkMuonsGmt"));
desc.add<edm::InputTag>("cscRecSegmentLabel", edm::InputTag("hltCscSegments"));
desc.add<edm::InputTag>("dtRecSegmentLabel", edm::InputTag("hltDt4DSegments"));
desc.add<double>("minPL1Tk", 3.5);
desc.add<double>("maxPL1Tk", 200);
desc.add<double>("stubMatchDPhi", 0.05);
desc.add<double>("stubMatchDTheta", 0.1);
desc.add<double>("extrapolationWindowClose", 0.1);
desc.add<double>("extrapolationWindowFar", 0.05);
desc.add<double>("maximumEtaBarrel", 0.7);
desc.add<double>("maximumEtaOverlap", 1.3);
desc.add<string>("Propagator", "SteppingHelixPropagatorAny");
desc.add<string>("propagator", "SteppingHelixPropagatorAny");

// Service parameters
edm::ParameterSetDescription psd0;
psd0.addUntracked<std::vector<std::string>>("Propagators", {"SteppingHelixPropagatorAny"});
psd0.add<bool>("RPCLayers", true);
psd0.addUntracked<bool>("UseMuonNavigation", true);
desc.add<edm::ParameterSetDescription>("ServiceParameters", psd0);
desc.add<edm::ParameterSetDescription>("serviceParameters", psd0);
desc.add<double>("estimatorMaxChi2", 1000.0);
descriptions.add("Phase2L2MuonSeedCreator", desc);
}

Expand All @@ -97,7 +100,8 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&

cscGeometry_ = iSetup.getHandle(cscGeometryToken_);
dtGeometry_ = iSetup.getHandle(dtGeometryToken_);
magneticField_ = iSetup.getHandle(magneticFieldToken_);

auto magneticFieldHandle = iSetup.getHandle(magneticFieldToken_);

LogDebug(metname) << "Number of L1 Tracker Muons in the Event: " << l1TkMuColl->size();

Expand All @@ -106,14 +110,25 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
l1t::TrackerMuonRef l1TkMuRef(l1TkMuColl, l1TkMuIndex);

// Physical info of L1TkMu
float pt = l1TkMuRef->phPt();
const float pt = l1TkMuRef->phPt();
if (pt < minMomentum_) {
continue;
}
float eta = l1TkMuRef->phEta();
float theta = 2 * std::atan(std::exp(-eta));
float phi = l1TkMuRef->phPhi();
int charge = l1TkMuRef->phCharge();
const float eta = l1TkMuRef->phEta();
const float phi = l1TkMuRef->phPhi();
const int charge = l1TkMuRef->phCharge();

// Calculate theta once to use it for stub-segment matching
// theta == 2 * atan(e^(-eta))
const float theta = 2 * vdt::fast_atanf(vdt::fast_expf(-eta));

// Precompute trig functions for theta
float sinTheta, cosTheta;
vdt::fast_sincosf(theta, sinTheta, cosTheta);

// Precompute trig functions for phi
float sinPhi, cosPhi;
vdt::fast_sincosf(phi, sinPhi, cosPhi);

LogDebug(metname) << "L1TKMu pT: " << pt << ", eta: " << eta << ", phi: " << phi;
Type muonType = overlap;
Expand Down Expand Up @@ -165,7 +180,7 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
stub->phiRegion() + 1); // sector, online to offline
LogDebug(metname) << "Stub DT detId: " << stubId << ". RawId: " << stubId.rawId();

auto& tmpMatch = matchingStubSegment(stubId, stub, dtSegments, l1TkMuRef);
auto& tmpMatch = matchingStubSegment(stubId, stub, dtSegments, theta);

// Found a match -> update matching info
if (tmpMatch.first != -1) {
Expand Down Expand Up @@ -196,7 +211,7 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
stub->phiRegion()); // chamber
LogDebug(metname) << "Stub CSC detId: " << stubId << ". RawId: " << stubId.rawId();

auto& tmpMatch = matchingStubSegment(stubId, stub, cscSegments, l1TkMuRef);
auto& tmpMatch = matchingStubSegment(stubId, stub, cscSegments, theta);

// Found a match -> update matching info
if (tmpMatch.first != -1) {
Expand Down Expand Up @@ -225,7 +240,7 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
stub->phiRegion() + 1); // sector, online to offline
LogDebug(metname) << "Stub DT detId: " << stubId << ". RawId: " << stubId.rawId();

auto& tmpMatch = matchingStubSegment(stubId, stub, dtSegments, l1TkMuRef);
auto& tmpMatch = matchingStubSegment(stubId, stub, dtSegments, theta);
totalBarrelQuality += tmpMatch.second;

// Found a match -> update matching info
Expand Down Expand Up @@ -254,7 +269,7 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
stub->phiRegion()); // chamber
LogDebug(metname) << "Stub CSC detId: " << stubId << ". RawId: " << stubId.rawId();

auto& tmpMatch = matchingStubSegment(stubId, stub, cscSegments, l1TkMuRef);
auto& tmpMatch = matchingStubSegment(stubId, stub, cscSegments, theta);
totalEndcapQuality += tmpMatch.second;

// Found a match -> update matching info
Expand Down Expand Up @@ -322,22 +337,18 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
const DetLayer* detLayer = nullptr;
float radius = 0.;

CLHEP::Hep3Vector vec(0., 1., 0.);
vec.setTheta(theta);
vec.setPhi(phi);

DetId propagateToId;

edm::OwnVector<TrackingRecHit> container;
if (bestInDt) {
// Found at least one matching segment in DT -> propagate to MB2
// Found at least one matching segment in DT -> propagate to Muon Barrel 2 (MB2) for track finding
LogDebug(metname) << "Found matching segment(s) in DTs, propagating L1TkMu info to MB2 to seed";
// MB2
// MB2 detId
propagateToId = DTChamberId(0, 2, 0);
detLayer = service_->detLayerGeometry()->idToLayer(propagateToId);
const BoundSurface* sur = &(detLayer->surface());
const BoundCylinder* bc = dynamic_cast<const BoundCylinder*>(sur);
radius = std::abs(bc->radius() / sin(theta));
radius = std::abs(bc->radius() / sinTheta);

// Propagate matched segments to the seed and try to extrapolate in unmatched chambers
std::vector<int> matchedStations;
Expand All @@ -348,7 +359,8 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
container.push_back(dtSegments[matchingPair.first]);
matchedStations.push_back(detId.station());
}
for (int station = 1; station != 5; ++station) {
// Loop over all barrel muon stations (1-4)
for (int station = DTChamberId::minStationId; station <= DTChamberId::maxStationId; ++station) {
if (std::find(matchedStations.begin(), matchedStations.end(), station) == matchedStations.end()) {
// Try to extrapolate from stations with a match to the ones without
LogDebug(metname) << "No matching DT segment found in station " << station;
Expand All @@ -361,25 +373,24 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
}
}
} else if (!bestInDt) {
// Found a matching segment in CSC -> propagate to ME2
// Found a matching segment in CSC -> propagate to Muon Endcap 2 (ME2) for track finding
LogDebug(metname) << "Found matching segment(s) in CSCs, propagating L1TkMu info to ME2 to seed";
// ME2
// ME2 detId
propagateToId = eta > 0 ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
detLayer = service_->detLayerGeometry()->idToLayer(propagateToId);
radius = std::abs(detLayer->position().z() / cos(theta));
radius = std::abs(detLayer->position().z() / cosTheta);

// Fill seed with matched segment(s)
for (auto& [detId, matchingPair] : matchesInEndcap) {
LogDebug(metname) << "Adding matched CSC segment in station " << detId.station() << " to the seed";
container.push_back(cscSegments[matchingPair.first]);
}
}
vec.setMag(radius);
// Get Global point and direction
GlobalPoint pos(vec.x(), vec.y(), vec.z());
GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
GlobalPoint pos(radius * cosPhi * sinTheta, radius * sinPhi * sinTheta, radius * cosTheta);
GlobalVector mom(pt * cosPhi, pt * sinPhi, pt * cosTheta / sinTheta);

GlobalTrajectoryParameters param(pos, mom, charge, &*magneticField_);
GlobalTrajectoryParameters param(pos, mom, charge, &*magneticFieldHandle);

AlgebraicSymMatrix55 mat;

Expand Down Expand Up @@ -407,8 +418,7 @@ void Phase2L2MuonSeedCreator::produce(edm::Event& iEvent, const edm::EventSetup&
// Propagation to MB2 failed, fallback to ME2 (might be an overlap edge case)
LogDebug(metname) << "Warning: detsWithStates collection is empty for a barrel collection. Falling back to ME2";
// Get ME2 DetLayer
DetId fallback_id;
theta < Geom::pi() / 2. ? fallback_id = CSCDetId(1, 2, 0, 0, 0) : fallback_id = CSCDetId(2, 2, 0, 0, 0);
DetId fallback_id = eta > 0 ? CSCDetId(1, 2, 0, 0, 0) : CSCDetId(2, 2, 0, 0, 0);
const DetLayer* ME2DetLayer = service_->detLayerGeometry()->idToLayer(fallback_id);
// Trajectory state on ME2 DetLayer
tsos = service_->propagator(propagatorName_)->propagate(state, ME2DetLayer->surface());
Expand Down Expand Up @@ -465,7 +475,7 @@ const std::vector<DTChamberId> Phase2L2MuonSeedCreator::matchingIds(const DTCham
const std::pair<int, int> Phase2L2MuonSeedCreator::matchingStubSegment(const DTChamberId& stubId,
const l1t::MuonStubRef stub,
const DTRecSegment4DCollection& segments,
const l1t::TrackerMuonRef l1TkMuRef) const {
const float l1TkMuTheta) const {
const std::string metname = "RecoMuon|Phase2L2MuonSeedCreator";

int bestSegIndex = -1;
Expand All @@ -492,7 +502,7 @@ const std::pair<int, int> Phase2L2MuonSeedCreator::matchingStubSegment(const DTC
double deltaPhi = std::abs(segPos.phi() - stub->offline_coord1());
LogDebug(metname) << "deltaPhi: " << deltaPhi;

double deltaTheta = std::abs(segPos.theta() - 2 * std::atan(std::exp(-l1TkMuRef->phEta())));
double deltaTheta = std::abs(segPos.theta() - l1TkMuTheta);
LogDebug(metname) << "deltaTheta: " << deltaTheta;

// Skip segments outside phi window or very far in the theta view
Expand Down Expand Up @@ -573,7 +583,7 @@ const std::vector<CSCDetId> Phase2L2MuonSeedCreator::matchingIds(const CSCDetId&
const std::pair<int, int> Phase2L2MuonSeedCreator::matchingStubSegment(const CSCDetId& stubId,
const l1t::MuonStubRef stub,
const CSCSegmentCollection& segments,
const l1t::TrackerMuonRef l1TkMuRef) const {
const float l1TkMuTheta) const {
const std::string metname = "RecoMuon|Phase2L2MuonSeedCreator";

int bestSegIndex = -1;
Expand All @@ -597,10 +607,14 @@ const std::pair<int, int> Phase2L2MuonSeedCreator::matchingStubSegment(const CSC
double deltaPhi = std::abs(segPos.phi() - stub->offline_coord1());
LogDebug(metname) << "deltaPhi: " << deltaPhi;

double deltaTheta = std::abs(segPos.theta() - 2 * std::atan(std::exp(-l1TkMuRef->phEta())));
double deltaTheta = std::abs(segPos.theta() - l1TkMuTheta);
LogDebug(metname) << "deltaTheta: " << deltaTheta;

if (deltaPhi > matchingPhiWindow_ or deltaTheta > 0.4) {
// Theta mainly used in cases where multiple matches are found
// to keep only the best one. Still skip segments way outside
// a reasonable window
const double roughThetaWindow = 0.4;
if (deltaPhi > matchingPhiWindow_ or deltaTheta > roughThetaWindow) {
continue; // skip segments outside phi window
}

Expand Down Expand Up @@ -791,7 +805,11 @@ const std::pair<int, int> Phase2L2MuonSeedCreator::extrapolateMatch(const int be
double matchingDeltaPhi =
std::abs(matchId.station() - endingStation) == 1 ? extrapolationDeltaPhiClose_ : extrapolationDeltaPhiFar_;

if (deltaPhi > matchingDeltaPhi or deltaTheta > 0.4) {
// Theta mainly used in cases where multiple matches are found
// to keep only the best one. Still skip segments way outside
// a reasonable window
const double roughThetaWindow = 0.4;
if (deltaPhi > matchingDeltaPhi or deltaTheta > roughThetaWindow) {
continue;
}

Expand Down
5 changes: 2 additions & 3 deletions RecoMuon/L2MuonSeedGenerator/src/Phase2L2MuonSeedCreator.h
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,6 @@ class Phase2L2MuonSeedCreator : public edm::stream::EDProducer<> {
const double maxEtaOverlap_; // overlap with |eta| < 1.3, endcap after that

// Handles
edm::ESHandle<MagneticField> magneticField_;
edm::ESHandle<CSCGeometry> cscGeometry_;
edm::ESHandle<DTGeometry> dtGeometry_;

Expand All @@ -132,13 +131,13 @@ class Phase2L2MuonSeedCreator : public edm::stream::EDProducer<> {
const std::pair<int, int> matchingStubSegment(const DTChamberId& stubId,
const l1t::MuonStubRef stub,
const DTRecSegment4DCollection& segments,
const l1t::TrackerMuonRef l1TkMuRef) const;
const float l1TkMuTheta) const;

// Logic to match L1 stubs to CSC segments
const std::pair<int, int> matchingStubSegment(const CSCDetId& stubId,
const l1t::MuonStubRef stub,
const CSCSegmentCollection& segments,
const l1t::TrackerMuonRef l1TkMuRef) const;
const float l1TkMuTheta) const;

// Logic to extrapolate from nearby stations in the barrel
const std::pair<int, int> extrapolateToNearbyStation(const int endingStation,
Expand Down
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#ifndef RecoMuon_L3TrackFinder_phase2HLTMuonSelectorForL3_H
#define RecoMuon_L3TrackFinder_phase2HLTMuonSelectorForL3_H

/** \class phase2HLTMuonSelectorForL3
/** \class Phase2HLTMuonSelectorForL3
*
* Phase-2 L3 selector for Muons
* This module allows to choose whether to perform
Expand Down Expand Up @@ -49,13 +49,13 @@ namespace edm {
class EventSetup;
} // End namespace edm

class phase2HLTMuonSelectorForL3 : public edm::stream::EDProducer<> {
class Phase2HLTMuonSelectorForL3 : public edm::stream::EDProducer<> {
public:
// Constructor
phase2HLTMuonSelectorForL3(const edm::ParameterSet&);
Phase2HLTMuonSelectorForL3(const edm::ParameterSet&);

// Destructor
~phase2HLTMuonSelectorForL3() override = default;
~Phase2HLTMuonSelectorForL3() override = default;

// Default values
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
Expand Down
Loading

0 comments on commit 62c4e64

Please sign in to comment.