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.

Sign up for GitHub

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

Add track-based isolation, recHit energy sum to IsolatedTrack class #31399

Closed
wants to merge 5 commits into from
Closed

Add track-based isolation, recHit energy sum to IsolatedTrack class #31399

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
de6661c
Add track-based isolation, recHit energy sum to IsolatedTrack class
Jun 16, 2020
a9e80ed
changes for code-format check
czkaiweb Sep 11, 2020
090a358
changes for code format change
czkaiweb Sep 11, 2020
e7acbcb
modified to deltaR2
czkaiweb Sep 14, 2020
28fd4ce
modified to std::abs
czkaiweb Sep 14, 2020
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
28 changes: 27 additions & 1 deletion DataFormats/PatCandidates/interface/IsolatedTrack.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,11 @@ namespace pat {
: LeafCandidate(0, LorentzVector(0, 0, 0, 0)),
pfIsolationDR03_(pat::PFIsolation()),
miniIsolation_(pat::PFIsolation()),
trackIsolation_(0.),
matchedCaloJetEmEnergy_(0.),
matchedCaloJetHadEnergy_(0.),
associatedEcalEnergy_(0.),
associatedHcalEnergy_(0.),
pfLepOverlap_(false),
pfNeutralSum_(0.),
dz_(0.),
Expand All @@ -47,8 +50,11 @@ namespace pat {

explicit IsolatedTrack(const PFIsolation& iso,
const PFIsolation& miniiso,
float trackiso,
float caloJetEm,
float caloJetHad,
float ecalEnergy,
float hcalEnergy,
bool pfLepOverlap,
float pfNeutralSum,
const LorentzVector& p4,
Expand All @@ -73,8 +79,11 @@ namespace pat {
: LeafCandidate(charge, p4, Point(0., 0., 0.), id),
pfIsolationDR03_(iso),
miniIsolation_(miniiso),
trackIsolation_(trackiso),
matchedCaloJetEmEnergy_(caloJetEm),
matchedCaloJetHadEnergy_(caloJetHad),
associatedEcalEnergy_(ecalEnergy),
associatedHcalEnergy_(hcalEnergy),
pfLepOverlap_(pfLepOverlap),
pfNeutralSum_(pfNeutralSum),
dz_(dz),
Expand All @@ -97,12 +106,15 @@ namespace pat {
~IsolatedTrack() override {}

const PFIsolation& pfIsolationDR03() const { return pfIsolationDR03_; }

const PFIsolation& miniPFIsolation() const { return miniIsolation_; }
float trackIsolation() const { return trackIsolation_; }

float matchedCaloJetEmEnergy() const { return matchedCaloJetEmEnergy_; }
float matchedCaloJetHadEnergy() const { return matchedCaloJetHadEnergy_; }

float associatedEcalEnergy() const { return associatedEcalEnergy_; }
float associatedHcalEnergy() const { return associatedHcalEnergy_; }

bool pfLepOverlap() const { return pfLepOverlap_; }
float pfNeutralSum() const { return pfNeutralSum_; }

Expand All @@ -121,6 +133,17 @@ namespace pat {

const reco::HitPattern& hitPattern() const { return hitPattern_; }

// helper functions for PhysicsTools/PatAlgos/python/slimming/isolatedTracks_cfi.py
int numMissingInnerHits() const {
return hitPattern_.trackerLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS);
}
int numMissingMiddleHits() const {
return hitPattern_.trackerLayersWithoutMeasurement(reco::HitPattern::TRACK_HITS);
}
int numMissingOuterHits() const {
return hitPattern_.trackerLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS);
}

float dEdxStrip() const { return dEdxStrip_; }
float dEdxPixel() const { return dEdxPixel_; }

Expand All @@ -141,8 +164,11 @@ namespace pat {
protected:
PFIsolation pfIsolationDR03_;
PFIsolation miniIsolation_;
float trackIsolation_;
float matchedCaloJetEmEnergy_; //energy of nearest calojet within a given dR;
float matchedCaloJetHadEnergy_;
float associatedEcalEnergy_; // sum of EBEE recHits within a given dR
float associatedHcalEnergy_; // sum of HBHE recHits within a given dR
bool pfLepOverlap_;
float pfNeutralSum_;
float dz_, dxy_, dzError_, dxyError_;
Expand Down
3 changes: 2 additions & 1 deletion DataFormats/PatCandidates/src/classes_def_objects.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,8 @@
<version ClassVersion="10" checksum="2244564938"/>
</class>

<class name="pat::IsolatedTrack" ClassVersion="5">
<class name="pat::IsolatedTrack" ClassVersion="6">
<version ClassVersion="6" checksum="2605383357"/>
<version ClassVersion="5" checksum="139818330"/>
<version ClassVersion="4" checksum="71166093"/>
<version ClassVersion="3" checksum="550880582"/>
Expand Down
148 changes: 139 additions & 9 deletions PhysicsTools/PatAlgos/plugins/PATIsolatedTrackProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@
#include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
#include "DataFormats/TrackReco/interface/DeDxData.h"
#include "DataFormats/TrackReco/interface/DeDxHitInfo.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "RecoTracker/DeDx/interface/DeDxTools.h"
Expand All @@ -39,6 +41,11 @@

#include "MagneticField/Engine/interface/MagneticField.h"

#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"

namespace pat {

class PATIsolatedTrackProducer : public edm::stream::EDProducer<> {
Expand All @@ -58,6 +65,11 @@ namespace pat {
pat::PFIsolation& iso,
pat::PFIsolation& miniiso) const;

float getTrackIsolation(const PolarLorentzVector& track_p4,
const math::XYZPoint& track_vtx,
const float& track_dzError,
const reco::TrackCollection* tracks) const;

bool getPFLeptonOverlap(const PolarLorentzVector& p4, const pat::PackedCandidateCollection* pc) const;

float getPFNeutralSum(const PolarLorentzVector& p4, const pat::PackedCandidateCollection* pc, int pc_idx) const;
Expand All @@ -73,6 +85,15 @@ namespace pat {

void getCaloJetEnergy(const PolarLorentzVector&, const reco::CaloJetCollection*, float&, float&) const;

void getAssociatedCaloEnergy(const PolarLorentzVector&,
const EBRecHitCollection&,
const EERecHitCollection&,
const HBHERecHitCollection&,
float&,
float&) const;

bool insideCone(const PolarLorentzVector&, const DetId&, const double) const;

private:
const edm::EDGetTokenT<pat::PackedCandidateCollection> pc_;
const edm::EDGetTokenT<pat::PackedCandidateCollection> lt_;
Expand All @@ -87,17 +108,24 @@ namespace pat {
const edm::EDGetTokenT<reco::DeDxHitInfoAss> gt2dedxHitInfo_;
const bool addPrescaledDeDxTracks_;
const edm::EDGetTokenT<edm::ValueMap<int>> gt2dedxHitInfoPrescale_;
const edm::EDGetTokenT<EBRecHitCollection> EBRecHits_;
const edm::EDGetTokenT<EERecHitCollection> EERecHits_;
const edm::EDGetTokenT<HBHERecHitCollection> HBHERecHits_;

const bool usePrecomputedDeDxStrip_;
const bool usePrecomputedDeDxPixel_;
const float pT_cut_; // only save cands with pT>pT_cut_
const float pT_cut_noIso_; // above this pT, don't apply any iso cut
const float pfIsolation_DR_; // isolation radius
const float pfIsolation_DZ_; // used in determining if pfcand is from PV or PU
const float absIso_cut_; // save if ANY of absIso, relIso, or miniRelIso pass the cuts
const float pT_cut_; // only save cands with pT>pT_cut_
const float pT_cut_noIso_; // above this pT, don't apply any iso cut
const float pfIsolation_DR_; // isolation radius
const float pfIsolation_DZ_; // used in determining if pfcand is from PV or PU
const float trackIsolation_DR_; // isolation radius
const float trackIsolation_maxDZSig_; // used in determining if two tracks are from different vertices
const float absIso_cut_; // save if ANY of absIso, relIso, or miniRelIso pass the cuts
const float relIso_cut_;
const float miniRelIso_cut_;
const float caloJet_DR_; // save energy of nearest calojet within caloJet_DR_
const float pflepoverlap_DR_; // pf lepton overlap radius
const float caloJet_DR_; // save energy of nearest calojet within caloJet_DR_
const float associatedCaloEnergy_DR_; // sum recHit energy within DR of track
const float pflepoverlap_DR_; // pf lepton overlap radius
const float pflepoverlap_pTmin_; // pf lepton overlap min pT (only look at PF candidates with pT>pflepoverlap_pTmin_)
const float pcRefNearest_DR_; // radius for nearest charged packed candidate
const float pcRefNearest_pTmin_; // min pT for nearest charged packed candidate
Expand All @@ -109,6 +137,8 @@ namespace pat {

TrackDetectorAssociator trackAssociator_;
TrackAssociatorParameters trackAssocParameters_;

edm::ESHandle<CaloGeometry> caloGeometry_;
};
} // namespace pat

Expand All @@ -131,16 +161,22 @@ pat::PATIsolatedTrackProducer::PATIsolatedTrackProducer(const edm::ParameterSet&
gt2dedxHitInfoPrescale_(addPrescaledDeDxTracks_ ? consumes<edm::ValueMap<int>>(
iConfig.getParameter<edm::InputTag>("dEdxHitInfoPrescale"))
: edm::EDGetTokenT<edm::ValueMap<int>>()),
EBRecHits_(consumes<EBRecHitCollection>(iConfig.getParameter<edm::InputTag>("EBRecHits"))),
EERecHits_(consumes<EERecHitCollection>(iConfig.getParameter<edm::InputTag>("EERecHits"))),
HBHERecHits_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("HBHERecHits"))),
usePrecomputedDeDxStrip_(iConfig.getParameter<bool>("usePrecomputedDeDxStrip")),
usePrecomputedDeDxPixel_(iConfig.getParameter<bool>("usePrecomputedDeDxPixel")),
pT_cut_(iConfig.getParameter<double>("pT_cut")),
pT_cut_noIso_(iConfig.getParameter<double>("pT_cut_noIso")),
pfIsolation_DR_(iConfig.getParameter<double>("pfIsolation_DR")),
pfIsolation_DZ_(iConfig.getParameter<double>("pfIsolation_DZ")),
trackIsolation_DR_(iConfig.getParameter<double>("trackIsolation_DR")),
trackIsolation_maxDZSig_(iConfig.getParameter<double>("trackIsolation_maxDZSig")),
absIso_cut_(iConfig.getParameter<double>("absIso_cut")),
relIso_cut_(iConfig.getParameter<double>("relIso_cut")),
miniRelIso_cut_(iConfig.getParameter<double>("miniRelIso_cut")),
caloJet_DR_(iConfig.getParameter<double>("caloJet_DR")),
associatedCaloEnergy_DR_(iConfig.getParameter<double>("associatedCaloEnergy_DR")),
pflepoverlap_DR_(iConfig.getParameter<double>("pflepoverlap_DR")),
pflepoverlap_pTmin_(iConfig.getParameter<double>("pflepoverlap_pTmin")),
pcRefNearest_DR_(iConfig.getParameter<double>("pcRefNearest_DR")),
Expand Down Expand Up @@ -229,6 +265,21 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event
iSetup.get<EcalChannelStatusRcd>().get(ecalS_h);
const EcalChannelStatus* ecalS = ecalS_h.product();

// get recHits collections for associated calo energy sums
edm::Handle<EBRecHitCollection> EBRecHits;
iEvent.getByToken(EBRecHits_, EBRecHits);

edm::Handle<EERecHitCollection> EERecHits;
iEvent.getByToken(EERecHits_, EERecHits);

edm::Handle<HBHERecHitCollection> HBHERecHits;
iEvent.getByToken(HBHERecHits_, HBHERecHits);

// get calorimeter geometry for recHit positions in associated calo energy sum calculation
iSetup.get<CaloGeometryRecord>().get(caloGeometry_);
if (!caloGeometry_.isValid())
throw cms::Exception("FatalError") << "Unable to find CaloGeometryRecord in event.\n";

auto outDeDxC = std::make_unique<reco::DeDxHitInfoCollection>();
std::vector<int> dEdXass;

Expand All @@ -254,8 +305,9 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event
pat::PackedCandidateRef refToCand;
int pdgId, charge, fromPV;
float dz, dxy, dzError, dxyError;
int pfCandInd; //to avoid counting packedPFCands in their own isolation
int ltCandInd; //to avoid pointing lost track to itself when looking for closest
int pfCandInd; //to avoid counting packedPFCands in their own isolation
int ltCandInd; //to avoid pointing lost track to itself when looking for closest
math::XYZPoint vtx; // for use in getTrackIsolation

// get the four-momentum and charge
if (isInPackedCands) {
Expand All @@ -264,12 +316,14 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event
charge = pfCand->charge();
pfCandInd = pcref.key();
ltCandInd = -1;
vtx = pfCand->vertex();
} else if (isInLostTracks) {
p4 = lostTrack->p4();
polarP4 = lostTrack->p4();
charge = lostTrack->charge();
pfCandInd = -1;
ltCandInd = ltref.key();
vtx = lostTrack->vertex();
} else {
double m = 0.13957018; //assume pion mass
double E = sqrt(m * m + gentk.p() * gentk.p());
Expand All @@ -278,6 +332,7 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event
charge = gentk.charge();
pfCandInd = -1;
ltCandInd = -1;
vtx = gentk.vertex();
}

int prescaled = 0;
Expand Down Expand Up @@ -332,9 +387,15 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event
refToCand = pat::PackedCandidateRef(); //NULL reference
}

// get the track isolation of the track
float trackIso = getTrackIsolation(polarP4, vtx, dzError, generalTracks);

float caloJetEm, caloJetHad;
getCaloJetEnergy(polarP4, caloJets.product(), caloJetEm, caloJetHad);

float ecalEnergy, hcalEnergy;
getAssociatedCaloEnergy(polarP4, *EBRecHits, *EERecHits, *HBHERecHits, ecalEnergy, hcalEnergy);

bool pfLepOverlap = getPFLeptonOverlap(polarP4, pc);
float pfNeutralSum = getPFNeutralSum(polarP4, pc, pfCandInd);

Expand Down Expand Up @@ -393,8 +454,11 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event

outPtrP->push_back(pat::IsolatedTrack(isolationDR03,
miniIso,
trackIso,
caloJetEm,
caloJetHad,
ecalEnergy,
hcalEnergy,
pfLepOverlap,
pfNeutralSum,
p4,
Expand Down Expand Up @@ -478,9 +542,15 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event
fromPV = pfCand.fromPV();
refToCand = pcref;

// get the track isolation of the track
float trackIso = getTrackIsolation(polarP4, pfCand.vertex(), dzError, generalTracks);

float caloJetEm, caloJetHad;
getCaloJetEnergy(polarP4, caloJets.product(), caloJetEm, caloJetHad);

float ecalEnergy, hcalEnergy;
getAssociatedCaloEnergy(polarP4, *EBRecHits, *EERecHits, *HBHERecHits, ecalEnergy, hcalEnergy);

bool pfLepOverlap = getPFLeptonOverlap(polarP4, pc);
float pfNeutralSum = getPFNeutralSum(polarP4, pc, ipc);

Expand All @@ -507,8 +577,11 @@ void pat::PATIsolatedTrackProducer::produce(edm::Event& iEvent, const edm::Event

outPtrP->push_back(pat::IsolatedTrack(isolationDR03,
miniIso,
trackIso,
caloJetEm,
caloJetHad,
ecalEnergy,
hcalEnergy,
pfLepOverlap,
pfNeutralSum,
pfCand.p4(),
Expand Down Expand Up @@ -592,6 +665,25 @@ void pat::PATIsolatedTrackProducer::getIsolation(const PolarLorentzVector& p4,
miniiso = pat::PFIsolation(chmiso, nhmiso, phmiso, pumiso);
}

float pat::PATIsolatedTrackProducer::getTrackIsolation(const PolarLorentzVector& track_p4,
const math::XYZPoint& track_vtx,
const float& track_dzError,
const reco::TrackCollection* tracks) const {
float sumPt = 0.;
for (const auto& t : *tracks) {
// exclude tracks from different vertices by dz signficance
// similar to primary vertex assignment, but includes dzError from both tracks
if (fabs(t.dz(track_vtx)) >= trackIsolation_maxDZSig_ * hypot(t.dzError(), track_dzError))
Copy link
Contributor

Choose a reason for hiding this comment

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

std::abs should be preferred to fabs
Could you please take the opportunity to replace also the other (two) fabs in this class?

Copy link
Author

Choose a reason for hiding this comment

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

changes made

continue;

float dR = deltaR(t, track_p4);
Copy link
Contributor

Choose a reason for hiding this comment

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

Please replace with deltaR2 (from the same header, it avoids uselessly computing a sqrt every time)

if (dR < trackIsolation_DR_ && dR > 1.e-12) // exclude candidate itself with dR>1.e-12
sumPt += t.pt();
}

return sumPt;
}

//get overlap of isolated track with a PF lepton
bool pat::PATIsolatedTrackProducer::getPFLeptonOverlap(const PolarLorentzVector& p4,
const pat::PackedCandidateCollection* pc) const {
Expand Down Expand Up @@ -777,6 +869,44 @@ void pat::PATIsolatedTrackProducer::getCaloJetEnergy(const PolarLorentzVector& p
}
}

void pat::PATIsolatedTrackProducer::getAssociatedCaloEnergy(const PolarLorentzVector& p4,
const EBRecHitCollection& EBRecHits,
const EERecHitCollection& EERecHits,
const HBHERecHitCollection& HBHERecHits,
float& eEM,
float& eHad) const {
eEM = 0.;
for (const auto& hit : EBRecHits) {
if (insideCone(p4, hit.detid(), associatedCaloEnergy_DR_))
eEM += hit.energy();
}
for (const auto& hit : EERecHits) {
if (insideCone(p4, hit.detid(), associatedCaloEnergy_DR_))
eEM += hit.energy();
}

eHad = 0.;
for (const auto& hit : HBHERecHits) {
if (insideCone(p4, hit.detid(), associatedCaloEnergy_DR_))
eHad += hit.energy();
}
}

bool pat::PATIsolatedTrackProducer::insideCone(const PolarLorentzVector& p4, const DetId& id, const double dR) const {
if (!caloGeometry_.isValid() || !caloGeometry_->getSubdetectorGeometry(id) ||
!caloGeometry_->getSubdetectorGeometry(id)->getGeometry(id)) {
throw cms::Exception("FatalError") << "Failed to access geometry for DetId: " << id.rawId();
return false;
}

const GlobalPoint& idPosition = caloGeometry_->getSubdetectorGeometry(id)->getGeometry(id)->getPosition();
if (idPosition.mag() < 0.01)
return false;

math::XYZVector idPositionRoot(idPosition.x(), idPosition.y(), idPosition.z());
return deltaR(p4, idPositionRoot) < dR;
Copy link
Contributor

Choose a reason for hiding this comment

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

This can also get simplified by using deltaR2 instead: just pass as argument the square of the dR
The *_DR parameters could be replaced by some *_DR2 ones, and as such the same simplification could be applied also (I think) in all other cases in which deltaR is used in this producer.

Copy link
Author

Choose a reason for hiding this comment

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

deltaR2 are used in our changes of the producer. *DR parameters are kept for convenience, explicit multiplication is made for DR parameters while comparing to dR2. Other cases from original code in which deltaR is used is left
for now.

}

using pat::PATIsolatedTrackProducer;
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(PATIsolatedTrackProducer);
Loading