Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added signed GSF track variable collections to the HLT producer: 14_0 #43774

Merged
merged 7 commits into from
Jan 25, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 13 additions & 1 deletion HLTrigger/Configuration/python/customizeHLTforCMSSW.py
Original file line number Diff line number Diff line change
Expand Up @@ -261,6 +261,18 @@ def customizeHLTfor43549(process):

return process

def customizeHLTfor43774(process):
filt_types = ["HLTEgammaGenericFilter","HLTEgammaGenericQuadraticEtaFilter","HLTEgammaGenericQuadraticFilter","HLTElectronGenericFilter"]
absAbleVar = ["DEta","deta","DetaSeed","Dphi","OneOESuperMinusOneOP","OneOESeedMinusOneOP"]
for filt_type in filt_types:
for filt in filters_by_type(process, filt_type):
if filt.varTag.productInstanceLabel in absAbleVar:
filt.useAbs = cms.bool(True)

return process
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

voicing here for the record a concern relayed from @missirol:
New trigger proposals might introduce new modules that do not comply with this customisation, and we might miss it (say they 'migrate' their config, and they do not realise that they have to put True for certain modules).

In the "usual" case, the default value from fillDescription is the correct one. We do not need the customisation, and people get the right value after they migrate their config (because ConfDB gives them the default value, which in this scenario is the correct one).

In "this" case, it looks like the correct value depends on the specific module and the value of its varTag, so I suspect some people might miss this if they migrate their Paths to 14_0_X (the default will be False, and they may not realise that they need to change this).

Bottom line: to ensure an extra line of defense, would it make sense to enforce at the C++ level that for the affected filters, if varTag is among the ones in this list the useAbs parameter is checked to be True (and we log an error or even throw otherwise) ?




# CMSSW version specific customizations
def customizeHLTforCMSSW(process, menuType="GRun"):

Expand All @@ -271,5 +283,5 @@ def customizeHLTforCMSSW(process, menuType="GRun"):

process = customizeHLTfor43025(process)
process = customizeHLTfor43549(process)

process = customizeHLTfor43774(process)
return process
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ HLTEgammaGenericQuadraticEtaFilter::HLTEgammaGenericQuadraticEtaFilter(const edm
energyLowEdges_ = iConfig.getParameter<std::vector<double> >("energyLowEdges");
lessThan_ = iConfig.getParameter<bool>("lessThan");
useEt_ = iConfig.getParameter<bool>("useEt");
useAbs_ = iConfig.getParameter<bool>("useAbs");

etaBoundaryEB12_ = iConfig.getParameter<double>("etaBoundaryEB12");
etaBoundaryEE12_ = iConfig.getParameter<double>("etaBoundaryEE12");
Expand Down Expand Up @@ -101,6 +102,7 @@ void HLTEgammaGenericQuadraticEtaFilter::fillDescriptions(edm::ConfigurationDesc
desc.add<std::vector<double> >("energyLowEdges", {0.0}); // No energy-dependent cuts by default
desc.add<bool>("lessThan", true);
desc.add<bool>("useEt", true);
desc.add<bool>("useAbs", false);
desc.add<double>("etaBoundaryEB12", 1.0);
desc.add<double>("etaBoundaryEE12", 2.0);
desc.add<std::vector<double> >("thrRegularEB1", {4.0});
Expand Down Expand Up @@ -175,7 +177,7 @@ bool HLTEgammaGenericQuadraticEtaFilter::hltFilter(edm::Event& iEvent,
ref = recoecalcands[i];
reco::RecoEcalCandidateIsolationMap::const_iterator mapi = (*depMap).find(ref);

float vali = mapi->val;
float vali = useAbs_ ? std::abs(mapi->val) : mapi->val;
float EtaSC = ref->eta();

// Pick the right EA and do rhoCorr
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ class HLTEgammaGenericQuadraticEtaFilter : public HLTFilter {
edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> varToken_;
bool lessThan_; // the cut is "<" or ">" ?
bool useEt_; // use E or Et in relative isolation cuts
bool useAbs_; // use the standard abs of the variable (before any rho corr)
/* Barrel quadratic threshold function:
vali (<= or >=) thrRegularEB_ + (E or Et)*thrOverEEB_ + (E or Et)*(E or Et)*thrOverE2EB_
Endcap quadratic threshold function:
Expand Down
4 changes: 3 additions & 1 deletion HLTrigger/Egamma/plugins/HLTEgammaGenericQuadraticFilter.cc
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ HLTEgammaGenericQuadraticFilter::HLTEgammaGenericQuadraticFilter(const edm::Para
energyLowEdges_ = iConfig.getParameter<std::vector<double> >("energyLowEdges");
lessThan_ = iConfig.getParameter<bool>("lessThan");
useEt_ = iConfig.getParameter<bool>("useEt");
useAbs_ = iConfig.getParameter<bool>("useAbs");

thrRegularEB_ = iConfig.getParameter<std::vector<double> >("thrRegularEB");
thrRegularEE_ = iConfig.getParameter<std::vector<double> >("thrRegularEE");
Expand Down Expand Up @@ -89,6 +90,7 @@ void HLTEgammaGenericQuadraticFilter::fillDescriptions(edm::ConfigurationDescrip
desc.add<std::vector<double> >("energyLowEdges", {0.0}); // No energy-dependent cuts by default
desc.add<bool>("lessThan", true);
desc.add<bool>("useEt", false);
desc.add<bool>("useAbs", false);
desc.add<std::vector<double> >("thrRegularEB", {0.0});
desc.add<std::vector<double> >("thrRegularEE", {0.0});
desc.add<std::vector<double> >("thrOverEEB", {-1.0});
Expand Down Expand Up @@ -155,7 +157,7 @@ bool HLTEgammaGenericQuadraticFilter::hltFilter(edm::Event& iEvent,
ref = recoecalcands[i];
reco::RecoEcalCandidateIsolationMap::const_iterator mapi = (*depMap).find(ref);

float vali = mapi->val;
float vali = useAbs_ ? std::abs(mapi->val) : mapi->val;
float EtaSC = ref->eta();

// Pick the right EA and do rhoCorr
Expand Down
1 change: 1 addition & 0 deletions HLTrigger/Egamma/plugins/HLTEgammaGenericQuadraticFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ class HLTEgammaGenericQuadraticFilter : public HLTFilter {

bool lessThan_; // the cut is "<" or ">" ?
bool useEt_; // use E or Et in relative isolation cuts
bool useAbs_; // use the standard abs of the variable (before any rho corr)
/* Barrel quadratic threshold function:
vali (<= or >=) thrRegularEB_ + (E or Et)*thrOverEEB_ + (E or Et)*(E or Et)*thrOverE2EB_
Endcap quadratic threshold function:
Expand Down
9 changes: 7 additions & 2 deletions HLTrigger/Egamma/plugins/HLTGenericFilter.cc
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ HLTGenericFilter<T1>::HLTGenericFilter(const edm::ParameterSet& iConfig) : HLTFi
energyLowEdges_ = iConfig.template getParameter<std::vector<double>>("energyLowEdges");
lessThan_ = iConfig.template getParameter<bool>("lessThan");
useEt_ = iConfig.template getParameter<bool>("useEt");
useAbs_ = iConfig.template getParameter<bool>("useAbs");
thrRegularEB_ = iConfig.template getParameter<std::vector<double>>("thrRegularEB");
thrRegularEE_ = iConfig.template getParameter<std::vector<double>>("thrRegularEE");
thrOverEEB_ = iConfig.template getParameter<std::vector<double>>("thrOverEEB");
Expand Down Expand Up @@ -88,6 +89,7 @@ void HLTGenericFilter<T1>::fillDescriptions(edm::ConfigurationDescriptions& desc
desc.add<std::vector<double>>("energyLowEdges", {0.0}); // No energy-dependent cuts by default
desc.add<bool>("lessThan", true);
desc.add<bool>("useEt", false);
desc.add<bool>("useAbs", false);
desc.add<std::vector<double>>("thrRegularEB", {0.0});
desc.add<std::vector<double>>("thrRegularEE", {0.0});
desc.add<std::vector<double>>("thrOverEEB", {-1.0});
Expand Down Expand Up @@ -177,7 +179,11 @@ bool HLTGenericFilter<T1>::hltFilter(edm::Event& iEvent,
T1Ref ref = recoCands[i];
typename T1IsolationMap::const_iterator mapi = (*depMap).find(ref);

float vali = mapi->val;
//should we do the abs before or after the rho corr?
//note: its very unlikely to rho correct a variable that wants to be abs
//decided yes as it seems slightly better this way but can not think of a use case either way
float vali = useAbs_ ? std::abs(mapi->val) : mapi->val;

float EtaSC = ref->eta();

// Pick the right EA and do rhoCorr
Expand All @@ -191,7 +197,6 @@ bool HLTGenericFilter<T1>::hltFilter(edm::Event& iEvent,
energy = getEt(ref);
else
energy = getEnergy(ref);
//if (energy < 0.) energy = 0.; // require energy to be positive (needed?)

// Pick the right cut threshold
double cutRegularEB_ = 9999., cutRegularEE_ = 9999.;
Expand Down
1 change: 1 addition & 0 deletions HLTrigger/Egamma/plugins/HLTGenericFilter.h
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ class HLTGenericFilter : public HLTFilter {
std::vector<double> energyLowEdges_; // lower bin edges for energy-dependent cuts
bool lessThan_; // the cut is "<" or ">" ?
bool useEt_; // use E or Et in relative isolation cuts
bool useAbs_; // use the standard abs of the variable (before any rho corr)
std::vector<double> thrRegularEB_; // threshold for regular cut (x < thr) - ECAL barrel
std::vector<double> thrRegularEE_; // threshold for regular cut (x < thr) - ECAL endcap
std::vector<double> thrOverEEB_; // threshold for x/E < thr cut (isolations) - ECAL barrel
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,29 @@
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"

class EgammaHLTGsfTrackVarProducer : public edm::global::EDProducer<> {
public:
struct GsfTrackExtrapolations {
GsfTrackExtrapolations() {}
void operator()(const reco::GsfTrack& trk,
const reco::SuperCluster& sc,
const MultiTrajectoryStateTransform& mtsTransform);
TrajectoryStateOnSurface innTSOS;
TrajectoryStateOnSurface outTSOS;
TrajectoryStateOnSurface sclTSOS;

GlobalVector innMom, outMom;
GlobalPoint sclPos;
};

public:
explicit EgammaHLTGsfTrackVarProducer(const edm::ParameterSet&);
void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
void fillAbsAbleVar(float& existVal, const float newVal) const {
if (std::abs(newVal) < std::abs(existVal)) {
existVal = produceAbsValues_ ? std::abs(newVal) : newVal;
}
}

private:
const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandToken_;
Expand All @@ -57,6 +76,7 @@ class EgammaHLTGsfTrackVarProducer : public edm::global::EDProducer<> {
const int lowerTrackNrToRemoveCut_;
const bool useDefaultValuesForBarrel_;
const bool useDefaultValuesForEndcap_;
const bool produceAbsValues_;

const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dEtaMapPutToken_;
const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> dEtaSeedMapPutToken_;
Expand All @@ -67,8 +87,15 @@ class EgammaHLTGsfTrackVarProducer : public edm::global::EDProducer<> {
const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> validHitsMapPutToken_;
const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> nLayerITMapPutToken_;
const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> chi2MapPutToken_;
const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> fbremMapPutToken_;
};

namespace {

float calRelDelta(float a, float b, float defaultVal = 0.f) { return a != 0.f ? (a - b) / a : defaultVal; }

} // namespace

EgammaHLTGsfTrackVarProducer::EgammaHLTGsfTrackVarProducer(const edm::ParameterSet& config)
: recoEcalCandToken_(
consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
Expand All @@ -81,6 +108,7 @@ EgammaHLTGsfTrackVarProducer::EgammaHLTGsfTrackVarProducer(const edm::ParameterS
lowerTrackNrToRemoveCut_{config.getParameter<int>("lowerTrackNrToRemoveCut")},
useDefaultValuesForBarrel_{config.getParameter<bool>("useDefaultValuesForBarrel")},
useDefaultValuesForEndcap_{config.getParameter<bool>("useDefaultValuesForEndcap")},
produceAbsValues_{config.getParameter<bool>("produceAbsValues")},
dEtaMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Deta").setBranchAlias("deta")},
dEtaSeedMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("DetaSeed").setBranchAlias("detaseed")},
dPhiMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Dphi").setBranchAlias("dphi")},
Expand All @@ -90,7 +118,8 @@ EgammaHLTGsfTrackVarProducer::EgammaHLTGsfTrackVarProducer(const edm::ParameterS
produces<reco::RecoEcalCandidateIsolationMap>("MissingHits").setBranchAlias("missinghits")},
validHitsMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("ValidHits").setBranchAlias("validhits")},
nLayerITMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("NLayerIT").setBranchAlias("nlayerit")},
chi2MapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Chi2").setBranchAlias("chi2")} {}
chi2MapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("Chi2").setBranchAlias("chi2")},
fbremMapPutToken_{produces<reco::RecoEcalCandidateIsolationMap>("fbrem")} {}

void EgammaHLTGsfTrackVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
Expand All @@ -101,6 +130,7 @@ void EgammaHLTGsfTrackVarProducer::fillDescriptions(edm::ConfigurationDescriptio
desc.add<int>(("lowerTrackNrToRemoveCut"), -1);
desc.add<bool>(("useDefaultValuesForBarrel"), false);
desc.add<bool>(("useDefaultValuesForEndcap"), false);
desc.add<bool>(("produceAbsValues"), false);

descriptions.add("hltEgammaHLTGsfTrackVarProducer", desc);
}
Expand All @@ -124,6 +154,7 @@ void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, co
reco::RecoEcalCandidateIsolationMap validHitsMap(recoEcalCandHandle);
reco::RecoEcalCandidateIsolationMap nLayerITMap(recoEcalCandHandle);
reco::RecoEcalCandidateIsolationMap chi2Map(recoEcalCandHandle);
reco::RecoEcalCandidateIsolationMap fbremMap(recoEcalCandHandle);

for (unsigned int iRecoEcalCand = 0; iRecoEcalCand < recoEcalCandHandle->size(); ++iRecoEcalCand) {
reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle, iRecoEcalCand);
Expand Down Expand Up @@ -156,6 +187,7 @@ void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, co
float dPhiInValue = 999999;
float oneOverESuperMinusOneOverPValue = 999999;
float oneOverESeedMinusOneOverPValue = 999999;
float fbrem = 999999;

const int nrTracks = gsfTracks.size();
const bool rmCutsDueToNrTracks = nrTracks <= lowerTrackNrToRemoveCut_ || nrTracks >= upperTrackNrToRemoveCut_;
Expand All @@ -164,6 +196,9 @@ void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, co
? useDefaultValuesForBarrel_ && nrTracks >= 1
: useDefaultValuesForEndcap_ && nrTracks >= 1;

MultiTrajectoryStateTransform mtsTransform(&trackerGeometry, &magneticField);
GsfTrackExtrapolations gsfTrackExtrapolations;

if (rmCutsDueToNrTracks || useDefaultValues) {
nLayerITValue = 100;
dEtaInValue = 0;
Expand All @@ -174,30 +209,23 @@ void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, co
chi2Value = 0;
oneOverESuperMinusOneOverPValue = 0;
oneOverESeedMinusOneOverPValue = 0;
fbrem = 0;
} else {
for (size_t trkNr = 0; trkNr < gsfTracks.size(); trkNr++) {
GlobalPoint scPos(scRef->x(), scRef->y(), scRef->z());

GlobalPoint trackExtrapToSC;
{
auto innTSOS =
MultiTrajectoryStateTransform::innerStateOnSurface(*gsfTracks[trkNr], trackerGeometry, &magneticField);
auto posTSOS = extrapolator.extrapolate(innTSOS, scPos);
rappoccio marked this conversation as resolved.
Show resolved Hide resolved
multiTrajectoryStateMode::positionFromModeCartesian(posTSOS, trackExtrapToSC);
}
gsfTrackExtrapolations(*gsfTracks[trkNr], *scRef, mtsTransform);

EleRelPointPair scAtVtx(scRef->position(), trackExtrapToSC, beamSpotPosition);
EleRelPointPair scAtVtx(scRef->position(), gsfTrackExtrapolations.sclPos, beamSpotPosition);

fbrem = calRelDelta(gsfTrackExtrapolations.innMom.mag(), gsfTrackExtrapolations.outMom.mag(), fbrem);

float trkP = gsfTracks[trkNr]->p();
if (scRef->energy() != 0 && trkP != 0) {
if (std::abs(1 / scRef->energy() - 1 / trkP) < oneOverESuperMinusOneOverPValue) {
oneOverESuperMinusOneOverPValue = std::abs(1 / scRef->energy() - 1 / trkP);
}
fillAbsAbleVar(oneOverESuperMinusOneOverPValue, 1 / scRef->energy() - 1 / trkP);
}
if (scRef->seed().isNonnull() && scRef->seed()->energy() != 0 && trkP != 0) {
if (std::abs(1 / scRef->seed()->energy() - 1 / trkP) < oneOverESeedMinusOneOverPValue) {
oneOverESeedMinusOneOverPValue = std::abs(1 / scRef->seed()->energy() - 1 / trkP);
}
fillAbsAbleVar(oneOverESeedMinusOneOverPValue, 1 / scRef->seed()->energy() - 1 / trkP);
}

if (gsfTracks[trkNr]->missingInnerHits() < missingHitsValue) {
Expand All @@ -218,19 +246,9 @@ void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, co
chi2Value = gsfTracks[trkNr]->normalizedChi2();
}

if (std::abs(scAtVtx.dEta()) < dEtaInValue) {
// we are allowing them to come from different tracks
dEtaInValue = std::abs(scAtVtx.dEta());
}

if (std::abs(scAtVtx.dEta()) < dEtaSeedInValue) {
dEtaSeedInValue = std::abs(scAtVtx.dEta() - scRef->position().eta() + scRef->seed()->position().eta());
}

if (std::abs(scAtVtx.dPhi()) < dPhiInValue) {
// we are allowing them to come from different tracks
dPhiInValue = std::abs(scAtVtx.dPhi());
}
fillAbsAbleVar(dEtaInValue, scAtVtx.dEta());
fillAbsAbleVar(dEtaSeedInValue, scAtVtx.dEta() - scRef->position().eta() + scRef->seed()->position().eta());
rappoccio marked this conversation as resolved.
Show resolved Hide resolved
fillAbsAbleVar(dPhiInValue, scAtVtx.dPhi());
}
}

Expand All @@ -243,6 +261,7 @@ void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, co
validHitsMap.insert(recoEcalCandRef, validHitsValue);
nLayerITMap.insert(recoEcalCandRef, nLayerITValue);
chi2Map.insert(recoEcalCandRef, chi2Value);
fbremMap.insert(recoEcalCandRef, fbrem);
}

iEvent.emplace(dEtaMapPutToken_, dEtaMap);
Expand All @@ -254,6 +273,18 @@ void EgammaHLTGsfTrackVarProducer::produce(edm::StreamID, edm::Event& iEvent, co
iEvent.emplace(validHitsMapPutToken_, validHitsMap);
iEvent.emplace(nLayerITMapPutToken_, nLayerITMap);
iEvent.emplace(chi2MapPutToken_, chi2Map);
iEvent.emplace(fbremMapPutToken_, fbremMap);
}

void EgammaHLTGsfTrackVarProducer::GsfTrackExtrapolations::operator()(
const reco::GsfTrack& trk, const reco::SuperCluster& sc, const MultiTrajectoryStateTransform& mtsTransform) {
innTSOS = mtsTransform.innerStateOnSurface(trk);
outTSOS = mtsTransform.outerStateOnSurface(trk);
sclTSOS = mtsTransform.extrapolatedState(innTSOS, GlobalPoint(sc.x(), sc.y(), sc.z()));

multiTrajectoryStateMode::momentumFromModeCartesian(innTSOS, innMom);
multiTrajectoryStateMode::positionFromModeCartesian(sclTSOS, sclPos);
multiTrajectoryStateMode::momentumFromModeCartesian(outTSOS, outMom);
}

#include "FWCore/Framework/interface/MakerMacros.h"
Expand Down