Skip to content

Commit

Permalink
Ecal rechit iso with GT thresholds, calotow Hcal/H0 iso
Browse files Browse the repository at this point in the history
  • Loading branch information
@24LopezR
24LopezR committed Jul 30, 2024
1 parent b943403 commit d3c29f2
Show file tree
Hide file tree
Showing 12 changed files with 317 additions and 61 deletions.
4 changes: 2 additions & 2 deletions RecoMuon/MuonIdentification/python/displacedMuons_cfi.py
View file
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
View file
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
View file
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
View file
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
View file
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

0 comments on commit d3c29f2

Please sign in to comment.