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

Muon Ecal rechit isolation thresh from GT, revert Hcal/H0 iso to calotower #45613

Merged
merged 4 commits into from
Aug 1, 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
7 changes: 7 additions & 0 deletions RecoMuon/MuonIdentification/plugins/MuonIdProducer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1532,6 +1532,13 @@ void MuonIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio
edm::ParameterSetDescription descCalo;
descCalo.setAllowAnything();
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this line setAllowAnything still needed, given you added the relevant parameters with values below?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I added below only the new parameters introduced in this PR, but there are many more. I do not know if it is safe to remove the setAllowAnything.

Copy link
Contributor

Choose a reason for hiding this comment

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

Ah OK - you'd have to add all of them ...

descCalo.add<edm::ParameterSetDescription>("TrackAssociatorParameters", descTrkAsoPar);
descCalo.add<bool>("UseEcalRecHitsFlag", false);
descCalo.add<bool>("UseHcalRecHitsFlag", false);
descCalo.add<bool>("UseHORecHitsFlag", false);
descCalo.add<bool>("EcalRecHitThresh", false);
descCalo.add<bool>("HcalCutsFromDB", false);
descCalo.add<int>("MaxSeverityHB", 9);
descCalo.add<int>("MaxSeverityHE", 9);
desc.add<edm::ParameterSetDescription>("CaloExtractorPSet", descCalo);

descriptions.addDefault(desc);
Expand Down
4 changes: 2 additions & 2 deletions RecoMuon/MuonIdentification/python/displacedMuons_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@

FillDetectorBasedIsolation = cms.bool(True),
EcalIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorHitsDisplaced","ecal"),
HcalIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorHitsDisplaced","hcal"),
HoIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorHitsDisplaced","ho"),
HcalIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorTowersDisplaced","hcal"),
HoIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorTowersDisplaced","ho"),
TrackIsoDeposits = cms.InputTag("muIsoDepositTkDisplaced"),
JetIsoDeposits = cms.InputTag("muIsoDepositJetsDisplaced"),

Expand Down
2 changes: 1 addition & 1 deletion RecoMuon/MuonIdentification/python/isolation_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
from RecoMuon.MuonIsolationProducers.jetExtractorBlock_cff import *
MIdIsoExtractorPSetBlock = cms.PSet(
CaloExtractorPSet = cms.PSet(
MIsoCaloExtractorByAssociatorHitsBlock
MIsoCaloExtractorByAssociatorMixedBlock
),
TrackExtractorPSet = cms.PSet(
MIsoTrackExtractorBlock
Expand Down
4 changes: 2 additions & 2 deletions RecoMuon/MuonIdentification/python/muons_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,8 @@

FillDetectorBasedIsolation = cms.bool(True),
EcalIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorHits","ecal"),
HcalIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorHits","hcal"),
HoIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorHits","ho"),
HcalIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorTowers","hcal"),
HoIsoDeposits = cms.InputTag("muIsoDepositCalByAssociatorTowers","ho"),
TrackIsoDeposits = cms.InputTag("muIsoDepositTk"),
JetIsoDeposits = cms.InputTag("muIsoDepositJets"),

Expand Down
6 changes: 6 additions & 0 deletions RecoMuon/MuonIsolation/plugins/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
<use name="FWCore/ParameterSet"/>
<use name="FWCore/PluginManager"/>
<use name="Geometry/CaloGeometry"/>
<use name="Geometry/CaloTopology"/>
<use name="Geometry/Records"/>
<use name="MagneticField/Engine"/>
<use name="MagneticField/Records"/>
Expand All @@ -24,5 +25,10 @@
<use name="TrackingTools/TrackAssociator"/>
<use name="TrackingTools/TransientTrack"/>
<use name="DataFormats/PatCandidates"/>
<use name="CondFormats/EcalObjects"/>
<use name="CondFormats/HcalObjects"/>
<use name="CondFormats/DataRecord"/>
<use name="CondTools/Hcal"/>
<use name="RecoLocalCalo/HcalRecAlgos"/>
<flags EDM_PLUGIN="1"/>
</library>
151 changes: 103 additions & 48 deletions RecoMuon/MuonIsolation/plugins/CaloExtractorByAssociator.cc
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,15 @@

#include "DataFormats/Math/interface/deltaR.h"

#include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
#include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
#include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
#include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
#include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
#include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"

using namespace edm;
using namespace std;
using namespace reco;
Expand All @@ -29,13 +38,17 @@ namespace {
}

CaloExtractorByAssociator::CaloExtractorByAssociator(const ParameterSet& par, edm::ConsumesCollector&& iC)
: theUseRecHitsFlag(par.getParameter<bool>("UseRecHitsFlag")),
: theUseEcalRecHitsFlag(par.getParameter<bool>("UseEcalRecHitsFlag")),
theUseHcalRecHitsFlag(par.getParameter<bool>("UseHcalRecHitsFlag")),
theUseHORecHitsFlag(par.getParameter<bool>("UseHORecHitsFlag")),
theDepositLabel(par.getUntrackedParameter<string>("DepositLabel")),
theDepositInstanceLabels(par.getParameter<std::vector<std::string> >("DepositInstanceLabels")),
thePropagatorName(par.getParameter<std::string>("PropagatorName")),
theThreshold_E(par.getParameter<double>("Threshold_E")),
theThreshold_H(par.getParameter<double>("Threshold_H")),
theThreshold_HO(par.getParameter<double>("Threshold_HO")),
theMaxSeverityHB(par.getParameter<int>("MaxSeverityHB")),
theMaxSeverityHE(par.getParameter<int>("MaxSeverityHE")),
theDR_Veto_E(par.getParameter<double>("DR_Veto_E")),
theDR_Veto_H(par.getParameter<double>("DR_Veto_H")),
theDR_Veto_HO(par.getParameter<double>("DR_Veto_HO")),
Expand All @@ -60,9 +73,15 @@ CaloExtractorByAssociator::CaloExtractorByAssociator(const ParameterSet& par, ed
theAssociatorParameters->loadParameters(par.getParameter<edm::ParameterSet>("TrackAssociatorParameters"), iC);
theAssociator = new TrackDetectorAssociator();

if (theUseRecHitsFlag) {
caloGeomToken_ = iC.esConsumes();
}
ecalRecHitThresh_ = par.getParameter<bool>("EcalRecHitThresh");
hcalCutsFromDB_ = par.getParameter<bool>("HcalCutsFromDB");

caloGeomToken_ = iC.esConsumes();
ecalPFRechitThresholdsToken_ = iC.esConsumes();
hcalCutsToken_ = iC.esConsumes();
hcalTopologyToken_ = iC.esConsumes();
hcalChannelQualityToken_ = iC.esConsumes(edm::ESInputTag("", "withTopo"));
hcalSevLvlComputerToken_ = iC.esConsumes();
}

CaloExtractorByAssociator::~CaloExtractorByAssociator() {
Expand Down Expand Up @@ -102,6 +121,12 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
theService->update(eventSetup);
theAssociator->setPropagator(&*(theService->propagator(thePropagatorName)));

const EcalPFRecHitThresholds* ecalThresholds = &eventSetup.getData(ecalPFRechitThresholdsToken_);
const HcalPFCuts* hcalCuts = &eventSetup.getData(hcalCutsToken_);
const HcalTopology* hcalTopology_ = &eventSetup.getData(hcalTopologyToken_);
const HcalChannelQuality* hcalChStatus_ = &eventSetup.getData(hcalChannelQualityToken_);
const HcalSeverityLevelComputer* hcalSevLvlComputer_ = &eventSetup.getData(hcalSevLvlComputerToken_);

//! check configuration consistency
//! could've been made at construction stage (fix later?)
if (theDepositInstanceLabels.size() != 3) {
Expand Down Expand Up @@ -155,29 +180,33 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
depHOcal.setVeto(Veto(dirTmp, dRtmp));
}

if (theUseRecHitsFlag) {
//! do things based on rec-hits here
//! too much copy-pasting now (refactor later?)
auto const& caloGeom = eventSetup.getData(caloGeomToken_);

if (theUseEcalRecHitsFlag) {
//Ecal
auto const& caloGeom = eventSetup.getData(caloGeomToken_);
std::vector<const EcalRecHit*>::const_iterator eHitCI = mInfo.ecalRecHits.begin();
for (; eHitCI != mInfo.ecalRecHits.end(); ++eHitCI) {
const EcalRecHit* eHitCPtr = *eHitCI;
GlobalPoint eHitPos = caloGeom.getPosition(eHitCPtr->detid());
double deltar0 = reco::deltaR(muon, eHitPos);
double deltaR2 = reco::deltaR2(muon, eHitPos);
double cosTheta = 1. / cosh(eHitPos.eta());
double energy = eHitCPtr->energy();
double et = energy * cosTheta;
if (deltar0 > std::max(dRMax_CandDep, theDR_Max) ||
!(et > theThreshold_E && energy > 3 * noiseRecHit(eHitCPtr->detid())))
if (deltaR2 > std::max(dRMax_CandDep * dRMax_CandDep, theDR_Max * theDR_Max))
continue;

if (ecalThresholds != nullptr) { // use thresholds from rechit
float rhThres = (ecalThresholds != nullptr) ? (*ecalThresholds)[eHitCPtr->detid()] : 0.f;
if (energy <= rhThres)
continue;
} else { // use thresholds from config
if (et <= theThreshold_E || energy <= 3 * noiseRecHit(eHitCPtr->detid()))
continue;
}

bool vetoHit = false;
double deltar = reco::deltaR(mInfo.trkGlobPosAtEcal, eHitPos);
//! first check if the hit is inside the veto cone by dR-alone
if (deltar < theDR_Veto_E) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto ECAL hit: Calo deltaR= " << deltar;
if (deltaR2 < std::pow(theDR_Veto_E, 2)) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto ECAL hit: Calo deltaR2= " << deltaR2;
LogDebug("RecoMuon|CaloExtractorByAssociator")
<< " >>> Calo eta phi ethcal: " << eHitPos.eta() << " " << eHitPos.phi() << " " << et;
vetoHit = true;
Expand All @@ -191,7 +220,7 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
}

//check theDR_Max only here to keep vetoHits being added to the veto energy
if (deltar0 > theDR_Max && !vetoHit)
if (deltaR2 > std::pow(theDR_Max, 2) && !vetoHit)
continue;

if (vetoHit) {
Expand All @@ -200,25 +229,48 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
depEcal.addDeposit(reco::isodeposit::Direction(eHitPos.eta(), eHitPos.phi()), et);
}
}
}

if (theUseHcalRecHitsFlag) {
//Hcal
auto const& caloGeom = eventSetup.getData(caloGeomToken_);
std::vector<const HBHERecHit*>::const_iterator hHitCI = mInfo.hcalRecHits.begin();
for (; hHitCI != mInfo.hcalRecHits.end(); ++hHitCI) {
const HBHERecHit* hHitCPtr = *hHitCI;
GlobalPoint hHitPos = caloGeom.getPosition(hHitCPtr->detid());
double deltar0 = reco::deltaR(muon, hHitPos);
double deltaR2 = reco::deltaR2(muon, hHitPos);
double cosTheta = 1. / cosh(hHitPos.eta());
double energy = hHitCPtr->energy();
double et = energy * cosTheta;
if (deltar0 > std::max(dRMax_CandDep, theDR_Max) ||
!(et > theThreshold_H && energy > 3 * noiseRecHit(hHitCPtr->detid())))
if (deltaR2 > std::max(dRMax_CandDep * dRMax_CandDep, theDR_Max * theDR_Max))
continue;

// check Hcal Cuts from DB
if (hcalCuts != nullptr) {
const HcalPFCut* item = hcalCuts->getValues(hHitCPtr->id().rawId());
if (energy <= item->noiseThreshold())
continue;
} else {
if (et <= theThreshold_H || energy <= 3 * noiseRecHit(hHitCPtr->detid()))
continue;
}

const HcalDetId hid(hHitCPtr->detid());
DetId did = hcalTopology_->idFront(hid);
const uint32_t flag = hHitCPtr->flags();
const uint32_t dbflag = hcalChStatus_->getValues(did)->getValue();
bool recovered = hcalSevLvlComputer_->recoveredRecHit(did, flag);
int severity = hcalSevLvlComputer_->getSeverityLevel(did, flag, dbflag);

const bool goodHB = hid.subdet() == HcalBarrel and (severity <= theMaxSeverityHB or recovered);
const bool goodHE = hid.subdet() == HcalEndcap and (severity <= theMaxSeverityHE or recovered);
if (!goodHB and !goodHE)
continue;

bool vetoHit = false;
double deltar = reco::deltaR(mInfo.trkGlobPosAtHcal, hHitPos);
//! first check if the hit is inside the veto cone by dR-alone
if (deltar < theDR_Veto_H) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HBHE hit: Calo deltaR= " << deltar;
if (deltaR2 < std::pow(theDR_Veto_H, 2)) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HBHE hit: Calo deltaR2= " << deltaR2;
LogDebug("RecoMuon|CaloExtractorByAssociator")
<< " >>> Calo eta phi ethcal: " << hHitPos.eta() << " " << hHitPos.phi() << " " << et;
vetoHit = true;
Expand All @@ -232,7 +284,7 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
}

//check theDR_Max only here to keep vetoHits being added to the veto energy
if (deltar0 > theDR_Max && !vetoHit)
if (deltaR2 > std::pow(theDR_Max, 2) && !vetoHit)
continue;

if (vetoHit) {
Expand All @@ -241,25 +293,27 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
depHcal.addDeposit(reco::isodeposit::Direction(hHitPos.eta(), hHitPos.phi()), et);
}
}
}

if (theUseHORecHitsFlag) {
//HOcal
auto const& caloGeom = eventSetup.getData(caloGeomToken_);
std::vector<const HORecHit*>::const_iterator hoHitCI = mInfo.hoRecHits.begin();
for (; hoHitCI != mInfo.hoRecHits.end(); ++hoHitCI) {
const HORecHit* hoHitCPtr = *hoHitCI;
GlobalPoint hoHitPos = caloGeom.getPosition(hoHitCPtr->detid());
double deltar0 = reco::deltaR(muon, hoHitPos);
double deltaR2 = reco::deltaR2(muon, hoHitPos);
double cosTheta = 1. / cosh(hoHitPos.eta());
double energy = hoHitCPtr->energy();
double et = energy * cosTheta;
if (deltar0 > std::max(dRMax_CandDep, theDR_Max) ||
if (deltaR2 > std::max(dRMax_CandDep * dRMax_CandDep, theDR_Max * theDR_Max) ||
!(et > theThreshold_HO && energy > 3 * noiseRecHit(hoHitCPtr->detid())))
continue;

bool vetoHit = false;
double deltar = reco::deltaR(mInfo.trkGlobPosAtHO, hoHitPos);
//! first check if the hit is inside the veto cone by dR-alone
if (deltar < theDR_Veto_HO) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HO hit: Calo deltaR= " << deltar;
if (deltaR2 < std::pow(theDR_Veto_HO, 2)) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HO hit: Calo deltaR2= " << deltaR2;
LogDebug("RecoMuon|CaloExtractorByAssociator")
<< " >>> Calo eta phi ethcal: " << hoHitPos.eta() << " " << hoHitPos.phi() << " " << et;
vetoHit = true;
Expand All @@ -273,7 +327,7 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
}

//check theDR_Max only here to keep vetoHits being added to the veto energy
if (deltar0 > theDR_Max && !vetoHit)
if (deltaR2 > std::pow(theDR_Max, 2) && !vetoHit)
continue;

if (vetoHit) {
Expand All @@ -282,14 +336,15 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
depHOcal.addDeposit(reco::isodeposit::Direction(hoHitPos.eta(), hoHitPos.phi()), et);
}
}
}

} else {
if (!theUseEcalRecHitsFlag or !theUseHcalRecHitsFlag or !theUseHORecHitsFlag) {
//! use calo towers
std::vector<const CaloTower*>::const_iterator calCI = mInfo.towers.begin();
for (; calCI != mInfo.towers.end(); ++calCI) {
const CaloTower* calCPtr = *calCI;
double deltar0 = reco::deltaR(muon, *calCPtr);
if (deltar0 > std::max(dRMax_CandDep, theDR_Max))
double deltaR2 = reco::deltaR2(muon, *calCPtr);
if (deltaR2 > std::max(dRMax_CandDep * dRMax_CandDep, theDR_Max * theDR_Max))
continue;

//even more copy-pasting .. need to refactor
Expand All @@ -306,28 +361,28 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
continue;

bool vetoTowerEcal = false;
double deltarEcal = reco::deltaR(mInfo.trkGlobPosAtEcal, *calCPtr);
double deltar2Ecal = reco::deltaR2(mInfo.trkGlobPosAtEcal, *calCPtr);
//! first check if the tower is inside the veto cone by dR-alone
if (deltarEcal < theDR_Veto_E) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto ecal tower: Calo deltaR= " << deltarEcal;
if (deltar2Ecal < std::pow(theDR_Veto_E, 2)) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto ecal tower: Calo deltaR= " << deltar2Ecal;
LogDebug("RecoMuon|CaloExtractorByAssociator")
<< " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
vetoTowerEcal = true;
}
bool vetoTowerHcal = false;
double deltarHcal = reco::deltaR(mInfo.trkGlobPosAtHcal, *calCPtr);
double deltar2Hcal = reco::deltaR2(mInfo.trkGlobPosAtHcal, *calCPtr);
//! first check if the tower is inside the veto cone by dR-alone
if (deltarHcal < theDR_Veto_H) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto hcal tower: Calo deltaR= " << deltarHcal;
if (deltar2Hcal < std::pow(theDR_Veto_H, 2)) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto hcal tower: Calo deltaR= " << deltar2Hcal;
LogDebug("RecoMuon|CaloExtractorByAssociator")
<< " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
vetoTowerHcal = true;
}
bool vetoTowerHOCal = false;
double deltarHOcal = reco::deltaR(mInfo.trkGlobPosAtHO, *calCPtr);
double deltar2HOcal = reco::deltaR2(mInfo.trkGlobPosAtHO, *calCPtr);
//! first check if the tower is inside the veto cone by dR-alone
if (deltarHOcal < theDR_Veto_HO) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HO tower: Calo deltaR= " << deltarHOcal;
if (deltar2HOcal < std::pow(theDR_Veto_HO, 2)) {
LogDebug("RecoMuon|CaloExtractorByAssociator") << " >>> Veto HO tower: Calo deltaR= " << deltar2HOcal;
LogDebug("RecoMuon|CaloExtractorByAssociator")
<< " >>> Calo eta phi ethcal: " << calCPtr->eta() << " " << calCPtr->phi() << " " << ethcal;
vetoTowerHOCal = true;
Expand All @@ -345,27 +400,27 @@ std::vector<IsoDeposit> CaloExtractorByAssociator::deposits(const Event& event,
}
}

if (deltar0 > theDR_Max && !(vetoTowerEcal || vetoTowerHcal || vetoTowerHOCal))
if (deltaR2 > std::pow(theDR_Max, 2) && !(vetoTowerEcal || vetoTowerHcal || vetoTowerHOCal))
continue;

reco::isodeposit::Direction towerDir(calCPtr->eta(), calCPtr->phi());
//! add the Et of the tower to deposits if it's not a vetoed; put into muonEnergy otherwise
if (doEcal) {
if (doEcal and !theUseEcalRecHitsFlag) {
if (vetoTowerEcal)
depEcal.addCandEnergy(etecal);
else if (deltar0 <= theDR_Max)
else if (deltaR2 <= std::pow(theDR_Max, 2))
depEcal.addDeposit(towerDir, etecal);
}
if (doHcal) {
if (doHcal and !theUseHcalRecHitsFlag) {
if (vetoTowerHcal)
depHcal.addCandEnergy(ethcal);
else if (deltar0 <= theDR_Max)
else if (deltaR2 <= std::pow(theDR_Max, 2))
depHcal.addDeposit(towerDir, ethcal);
}
if (doHOcal) {
if (doHOcal and !theUseHORecHitsFlag) {
if (vetoTowerHOCal)
depHOcal.addCandEnergy(ethocal);
else if (deltar0 <= theDR_Max)
else if (deltaR2 <= std::pow(theDR_Max, 2))
depHOcal.addDeposit(towerDir, ethocal);
}
}
Expand Down
Loading