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

Tau cluster variable refactoring: Avoid code duplication #83

Merged
merged 2 commits into from
Sep 3, 2018
Merged

Tau cluster variable refactoring: Avoid code duplication #83

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
eef7b80
Address code review comments, including modernisation of code
steggema Aug 30, 2018
dad7231
Create template function to avoid code duplication
steggema Aug 30, 2018
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
140 changes: 121 additions & 19 deletions RecoTauTag/RecoTau/interface/PFRecoTauClusterVariables.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,38 +16,140 @@
#include "DataFormats/PatCandidates/interface/Tau.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"

namespace reco { namespace tau { namespace mva {
namespace {
template<class TauType, class PhotonVectorType>
PhotonVectorType getGammas(const TauType& tau, bool signal);

/// return pf photon candidates that are associated to signal
template<>
const std::vector<reco::PFCandidatePtr>& getGammas<reco::PFTau, const std::vector<reco::PFCandidatePtr>&>(const reco::PFTau& tau, bool signal) {
if (signal){
return tau.signalPFGammaCands();
}
return tau.isolationPFGammaCands();
}

template<>
reco::CandidatePtrVector getGammas<pat::Tau, reco::CandidatePtrVector>(const pat::Tau& tau, bool signal) {
if(signal){
return tau.signalGammaCands();
}
return tau.isolationGammaCands();
}

/// decide if photon candidate is inside the cone to be associated to the tau signal
bool isInside(float photon_pt, float deta, float dphi) {
constexpr double stripEtaAssociationDistance_0p95_p0 = 0.197077;
constexpr double stripEtaAssociationDistance_0p95_p1 = 0.658701;
constexpr double stripPhiAssociationDistance_0p95_p0 = 0.352476;
constexpr double stripPhiAssociationDistance_0p95_p1 = 0.707716;
if(photon_pt==0){
return false;
} if((dphi<0.3 && dphi<std::max(0.05, stripPhiAssociationDistance_0p95_p0*std::pow(photon_pt, -stripPhiAssociationDistance_0p95_p1))) && (deta<0.15 && deta<std::max(0.05, stripEtaAssociationDistance_0p95_p0*std::pow(photon_pt, -stripEtaAssociationDistance_0p95_p1)))){
return true;
}
return false;
}
}

namespace reco { namespace tau {
/// return chi2 of the leading track ==> deprecated? <==
float tau_leadTrackChi2(const reco::PFTau& tau);
float lead_track_chi2(const reco::PFTau& tau);
/// return ratio of energy in ECAL over sum of energy in ECAL and HCAL
float tau_Eratio(const reco::PFTau& tau);
float tau_Eratio(const pat::Tau& tau);
float eratio(const reco::PFTau& tau);
float eratio(const pat::Tau& tau);
/// return sum of pt weighted values of distance to tau candidate for all pf photon candidates,
/// which are associated to signal; depending on var the distance is in 0=:dr, 1=:deta, 2=:dphi
float pt_weighted_dx(const reco::PFTau& tau, int mode = 0, int var = 0, int decaymode = -1);
float pt_weighted_dx(const pat::Tau& tau, int mode = 0, int var = 0, int decaymode = -1);
template<class TauType, class PhotonVectorType>
float pt_weighted_dx(const TauType& tau, int mode = 0, int var = 0, int decaymode = -1) {
float sum_pt = 0.;
float sum_dx_pt = 0.;
float signalrad = std::max(0.05, std::min(0.1, 3./std::max(1., tau.pt())));
int is3prong = (decaymode==10);
const auto& cands = getGammas<TauType, PhotonVectorType>(tau, mode < 2);
for (const auto& cand : cands) {
// only look at electrons/photons with pT > 0.5
if (cand->pt() < 0.5){
continue;
}
float dr = reco::deltaR(cand->eta(), cand->phi(), tau.eta(), tau.phi());
float deta = std::abs(cand->eta() - tau.eta());
float dphi = std::abs(reco::deltaPhi(cand->phi(), tau.phi()));
float pt = cand->pt();
bool flag = isInside(pt, deta, dphi);
if(is3prong==0){
if (mode == 2 || (mode == 0 && dr < signalrad) || (mode == 1 && dr > signalrad)) {
sum_pt += pt;
if (var == 0)
sum_dx_pt += pt * dr;
else if (var == 1)
sum_dx_pt += pt * deta;
else if (var == 2)
sum_dx_pt += pt * dphi;
}
}
else if(is3prong==1){
if( (mode==2 && flag==false) || (mode==1 && flag==true) || mode==0){
sum_pt += pt;
if (var == 0)
sum_dx_pt += pt * dr;
else if (var == 1)
sum_dx_pt += pt * deta;
else if (var == 2)
sum_dx_pt += pt * dphi;
}
}
}
if (sum_pt > 0.){
return sum_dx_pt/sum_pt;
}
return 0.;
}
/// return sum of pt weighted values of dr relative to tau candidate for all pf photon candidates,
/// which are associated to signal
float tau_pt_weighted_dr_signal(const reco::PFTau& tau, int dm);
float tau_pt_weighted_dr_signal(const pat::Tau& tau, int dm);
inline float pt_weighted_dr_signal(const reco::PFTau& tau, int dm) {
return pt_weighted_dx<reco::PFTau, const std::vector<reco::PFCandidatePtr>&>(tau, 0, 0, dm);
}
inline float pt_weighted_dr_signal(const pat::Tau& tau, int dm) {
return pt_weighted_dx<pat::Tau, reco::CandidatePtrVector>(tau, 0, 0, dm);
}
/// return sum of pt weighted values of deta relative to tau candidate for all pf photon candidates,
/// which are associated to signal
float tau_pt_weighted_deta_strip(const reco::PFTau& tau, int dm);
float tau_pt_weighted_deta_strip(const pat::Tau& tau, int dm);
inline float pt_weighted_deta_strip(const reco::PFTau& tau, int dm) {
return pt_weighted_dx<reco::PFTau, const std::vector<reco::PFCandidatePtr>&>(tau, dm==10 ? 2 : 1, 1, dm);
}
inline float pt_weighted_deta_strip(const pat::Tau& tau, int dm) {
return pt_weighted_dx<pat::Tau, reco::CandidatePtrVector>(tau, dm==10 ? 2 : 1, 1, dm);
}
/// return sum of pt weighted values of dphi relative to tau candidate for all pf photon candidates,
/// which are associated to signal
float tau_pt_weighted_dphi_strip(const reco::PFTau& tau, int dm);
float tau_pt_weighted_dphi_strip(const pat::Tau& tau, int dm);
inline float pt_weighted_dphi_strip(const reco::PFTau& tau, int dm) {
return pt_weighted_dx<reco::PFTau, const std::vector<reco::PFCandidatePtr>&>(tau, dm==10 ? 2 : 1, 2, dm);
}
inline float pt_weighted_dphi_strip(const pat::Tau& tau, int dm) {
return pt_weighted_dx<pat::Tau, reco::CandidatePtrVector>(tau, dm==10 ? 2 : 1, 2, dm);
}
/// return sum of pt weighted values of dr relative to tau candidate for all pf photon candidates,
/// which are inside an isolation conde but not associated to signal
inline float pt_weighted_dr_iso(const reco::PFTau& tau, int dm) {
return pt_weighted_dx<reco::PFTau, const std::vector<reco::PFCandidatePtr>&>(tau, 2, 0, dm);
}
inline float pt_weighted_dr_iso(const pat::Tau& tau, int dm) {
return pt_weighted_dx<pat::Tau, reco::CandidatePtrVector>(tau, 2, 0, dm);
}
/// return sum of pt weighted values of dr relative to tau candidate for all pf photon candidates,
/// which are inside an isolation conde but not associated to signal
float tau_pt_weighted_dr_iso(const reco::PFTau& tau, int dm);
float tau_pt_weighted_dr_iso(const pat::Tau& tau, int dm);
float pt_weighted_dr_iso(const reco::PFTau& tau, int dm);
float pt_weighted_dr_iso(const pat::Tau& tau, int dm);
/// return total number of pf photon candidates with pT>500 MeV, which are associated to signal
unsigned int tau_n_photons_total(const reco::PFTau& tau);
unsigned int tau_n_photons_total(const pat::Tau& tau);
unsigned int n_photons_total(const reco::PFTau& tau);
unsigned int n_photons_total(const pat::Tau& tau);

enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT, kDBoldDMwLTwGJ, kDBnewDMwLTwGJ };
bool fillMVAInputs(float* mvaInput, const pat::Tau& tau, int mvaOpt, const std::string nameCharged, const std::string nameNeutral, const std::string namePu, const std::string nameOutside, const std::string nameFootprint);
}}} // namespaces
enum {kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT,
kDBoldDMwLTwGJ, kDBnewDMwLTwGJ};
bool fillIsoMVARun2Inputs(float* mvaInput, const pat::Tau& tau, int mvaOpt, const std::string& nameCharged,
const std::string& nameNeutral, const std::string& namePu,
const std::string& nameOutside, const std::string& nameFootprint);
}} // namespaces

#endif
7 changes: 5 additions & 2 deletions RecoTauTag/RecoTau/plugins/PATTauDiscriminationByMVAIsolationRun2.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,8 @@ namespace
}
}

namespace reco { namespace tau {

class PATTauDiscriminationByMVAIsolationRun2 : public PATTauDiscriminationProducerBase
{
public:
Expand Down Expand Up @@ -143,7 +145,6 @@ class PATTauDiscriminationByMVAIsolationRun2 : public PATTauDiscriminationProduc
bool loadMVAfromDB_;
edm::FileInPath inputFileName_;
const GBRForest* mvaReader_;
enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT, kDBoldDMwLTwGJ, kDBnewDMwLTwGJ };
int mvaOpt_;
float* mvaInput_;

Expand Down Expand Up @@ -183,7 +184,7 @@ double PATTauDiscriminationByMVAIsolationRun2::discriminate(const TauRef& tau) c
// CV: computation of MVA value requires presence of leading charged hadron
if ( tau->leadChargedHadrCand().isNull() ) return 0.;

if (reco::tau::mva::fillMVAInputs(mvaInput_, *tau, mvaOpt_, chargedIsoPtSums_, neutralIsoPtSums_, puCorrPtSums_, photonPtSumOutsideSignalCone_, footprintCorrection_)) {
if (reco::tau::fillIsoMVARun2Inputs(mvaInput_, *tau, mvaOpt_, chargedIsoPtSums_, neutralIsoPtSums_, puCorrPtSums_, photonPtSumOutsideSignalCone_, footprintCorrection_)) {
double mvaValue = mvaReader_->GetClassifier(mvaInput_);
if ( verbosity_ ) {
edm::LogPrint("PATTauDiscByMVAIsolRun2") << "<PATTauDiscriminationByMVAIsolationRun2::discriminate>:";
Expand All @@ -207,3 +208,5 @@ void PATTauDiscriminationByMVAIsolationRun2::endEvent(edm::Event& evt)
}

DEFINE_FWK_MODULE(PATTauDiscriminationByMVAIsolationRun2);

}} //namespace
19 changes: 11 additions & 8 deletions RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByMVAIsolationRun2.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@ namespace
}
}

namespace reco { namespace tau {

class PFRecoTauDiscriminationByMVAIsolationRun2 : public PFTauDiscriminationProducerBase
{
public:
Expand Down Expand Up @@ -138,7 +140,6 @@ class PFRecoTauDiscriminationByMVAIsolationRun2 : public PFTauDiscriminationProd
bool loadMVAfromDB_;
edm::FileInPath inputFileName_;
const GBRForest* mvaReader_;
enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT, kDBoldDMwLTwGJ, kDBnewDMwLTwGJ };
int mvaOpt_;
float* mvaInput_;

Expand Down Expand Up @@ -218,13 +219,13 @@ double PFRecoTauDiscriminationByMVAIsolationRun2::discriminate(const PFTauRef& t
float decayDistZ = tauLifetimeInfo.flightLength().z();
float decayDistMag = std::sqrt(decayDistX*decayDistX + decayDistY*decayDistY + decayDistZ*decayDistZ);

float nPhoton = (float)reco::tau::mva::tau_n_photons_total(*tau);
float ptWeightedDetaStrip = reco::tau::mva::tau_pt_weighted_deta_strip(*tau, tauDecayMode);
float ptWeightedDphiStrip = reco::tau::mva::tau_pt_weighted_dphi_strip(*tau, tauDecayMode);
float ptWeightedDrSignal = reco::tau::mva::tau_pt_weighted_dr_signal(*tau, tauDecayMode);
float ptWeightedDrIsolation = reco::tau::mva::tau_pt_weighted_dr_iso(*tau, tauDecayMode);
float leadingTrackChi2 = reco::tau::mva::tau_leadTrackChi2(*tau);
float eRatio = reco::tau::mva::tau_Eratio(*tau);
float nPhoton = (float)reco::tau::n_photons_total(*tau);
float ptWeightedDetaStrip = reco::tau::pt_weighted_deta_strip(*tau, tauDecayMode);
float ptWeightedDphiStrip = reco::tau::pt_weighted_dphi_strip(*tau, tauDecayMode);
float ptWeightedDrSignal = reco::tau::pt_weighted_dr_signal(*tau, tauDecayMode);
float ptWeightedDrIsolation = reco::tau::pt_weighted_dr_iso(*tau, tauDecayMode);
float leadingTrackChi2 = reco::tau::lead_track_chi2(*tau);
float eRatio = reco::tau::eratio(*tau);

// Difference between measured and maximally allowed Gottfried-Jackson angle
float gjAngleDiff = -999;
Expand Down Expand Up @@ -359,3 +360,5 @@ void PFRecoTauDiscriminationByMVAIsolationRun2::endEvent(edm::Event& evt)
}

DEFINE_FWK_MODULE(PFRecoTauDiscriminationByMVAIsolationRun2);

}} //namespace
Loading