45 std::vector<TFile*>& inputFilesToDelete) {
47 throw cms::Exception(
"PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
48 <<
" Failed to find File = " << inputFileName <<
" !!\n";
54 throw cms::Exception(
"PFRecoTauDiscriminationByIsolationMVA2::loadMVA")
55 <<
" Failed to load MVA = " << mvaName.data() <<
" from file = " << inputFileName.
fullPath().data()
58 inputFilesToDelete.push_back(inputFile);
75 loadMVAfromDB_ = cfg.
getParameter<
bool>(
"loadMVAfromDB");
76 if (!loadMVAfromDB_) {
82 if (mvaOpt_string ==
"oldDMwoLT")
84 else if (mvaOpt_string ==
"oldDMwLT")
86 else if (mvaOpt_string ==
"newDMwoLT")
88 else if (mvaOpt_string ==
"newDMwLT")
90 else if (mvaOpt_string ==
"DBoldDMwLT")
92 else if (mvaOpt_string ==
"DBnewDMwLT")
94 else if (mvaOpt_string ==
"PWoldDMwLT")
96 else if (mvaOpt_string ==
"PWnewDMwLT")
98 else if (mvaOpt_string ==
"DBoldDMwLTwGJ")
100 else if (mvaOpt_string ==
"DBnewDMwLTwGJ")
104 <<
" Invalid Configuration Parameter 'mvaOpt' = " << mvaOpt_string <<
" !!\n";
107 mvaInput_ =
new float[6];
109 mvaInput_ =
new float[12];
112 mvaInput_ =
new float[23];
116 TauTransverseImpactParameters_token =
119 basicTauDiscriminators_token =
134 for (std::vector<TFile*>::iterator it = inputFilesToDelete_.begin(); it != inputFilesToDelete_.end(); ++it) {
159 int chargedIsoPtSum_index_ = 0;
160 int neutralIsoPtSum_index_ = 0;
161 int pucorrPtSum_index_ = 0;
162 int photonPtSumOutsideSignalCone_index_ = 0;
163 int footprintCorrection_index_ = 0;
176 if (loadMVAfromDB_) {
177 mvaReader_ = &es.
getData(mvaToken_);
179 mvaReader_ = loadMVAfromFile(inputFileName_, mvaName_, inputFilesToDelete_);
183 evt.
getByToken(TauTransverseImpactParameters_token, tauLifetimeInfos);
185 evt.
getByToken(basicTauDiscriminators_token, basicTauDiscriminators_);
194 const std::vector<edm::ParameterSet> psetsFromProvenance =
196 .getParameter<std::vector<edm::ParameterSet>>(
"IDdefinitions");
197 for (
uint i = 0;
i < psetsFromProvenance.size();
i++) {
198 if (psetsFromProvenance[
i].getParameter<std::string>(
"IDname") ==
"ChargedIsoPtSum" + input_id_name_suffix_)
199 chargedIsoPtSum_index_ =
i;
200 else if (psetsFromProvenance[
i].getParameter<std::string>(
"IDname") ==
201 "NeutralIsoPtSum" + input_id_name_suffix_)
202 neutralIsoPtSum_index_ =
i;
203 else if (psetsFromProvenance[
i].getParameter<std::string>(
"IDname") ==
"PUcorrPtSum" + input_id_name_suffix_)
204 pucorrPtSum_index_ =
i;
205 else if (psetsFromProvenance[
i].getParameter<std::string>(
"IDname") ==
206 "PhotonPtSumOutsideSignalCone" + input_id_name_suffix_)
207 photonPtSumOutsideSignalCone_index_ =
i;
208 else if (psetsFromProvenance[
i].getParameter<std::string>(
"IDname") ==
209 "TauFootprintCorrection" + input_id_name_suffix_)
210 footprintCorrection_index_ =
i;
221 if (tau->leadChargedHadrCand().
isNull())
224 int tauDecayMode = tau->decayMode();
228 (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
231 (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 ||
232 tauDecayMode == 10 || tauDecayMode == 11))) {
233 auto const rawValues = (*basicTauDiscriminators_)[
tau].rawValues;
234 float chargedIsoPtSum =
rawValues.at(chargedIsoPtSum_index_);
235 float neutralIsoPtSum =
rawValues.at(neutralIsoPtSum_index_);
236 float puCorrPtSum =
rawValues.at(pucorrPtSum_index_);
237 float photonPtSumOutsideSignalCone =
rawValues.at(photonPtSumOutsideSignalCone_index_);
238 float footprintCorrection =
rawValues.at(footprintCorrection_index_);
245 float decayDistMag =
std::sqrt(decayDistX * decayDistX + decayDistY * decayDistY + decayDistZ * decayDistZ);
256 float gjAngleDiff = -999;
257 if (tauDecayMode == 10) {
258 double mTau = 1.77682;
259 double mAOne = tau->p4().M();
260 double pAOneMag = tau->p();
261 double argumentThetaGJmax = (
std::pow(mTau, 2) -
std::pow(mAOne, 2)) / (2 * mTau * pAOneMag);
262 double argumentThetaGJmeasured =
263 (tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ) /
264 (pAOneMag * decayDistMag);
265 if (
std::abs(argumentThetaGJmax) <= 1. &&
std::abs(argumentThetaGJmeasured) <= 1.) {
266 double thetaGJmax = std::asin(argumentThetaGJmax);
267 double thetaGJmeasured = std::acos(argumentThetaGJmeasured);
268 gjAngleDiff = thetaGJmeasured - thetaGJmax;
274 mvaInput_[1] =
std::abs((
float)tau->eta());
278 mvaInput_[5] = tauDecayMode;
281 mvaInput_[1] =
std::abs((
float)tau->eta());
285 mvaInput_[5] = tauDecayMode;
286 mvaInput_[6] = std::copysign(+1.
f, (
float)tauLifetimeInfo.
dxy());
294 mvaInput_[1] =
std::abs((
float)tau->eta());
299 mvaInput_[6] = tauDecayMode;
301 mvaInput_[8] =
std::min(0.5
f, ptWeightedDetaStrip);
302 mvaInput_[9] =
std::min(0.5
f, ptWeightedDphiStrip);
303 mvaInput_[10] =
std::min(0.5
f, ptWeightedDrSignal);
304 mvaInput_[11] =
std::min(0.5
f, ptWeightedDrIsolation);
305 mvaInput_[12] =
std::min(100.
f, leadingTrackChi2);
307 mvaInput_[14] = std::copysign(+1.
f, (
float)tauLifetimeInfo.
dxy());
310 mvaInput_[17] = std::copysign(+1.
f, (
float)tauLifetimeInfo.
ip3d());
318 mvaInput_[1] =
std::abs((
float)tau->eta());
323 mvaInput_[6] = tauDecayMode;
325 mvaInput_[8] =
std::min(0.5
f, ptWeightedDetaStrip);
326 mvaInput_[9] =
std::min(0.5
f, ptWeightedDphiStrip);
327 mvaInput_[10] =
std::min(0.5
f, ptWeightedDrSignal);
328 mvaInput_[11] =
std::min(0.5
f, ptWeightedDrIsolation);
329 mvaInput_[12] =
std::min(100.
f, leadingTrackChi2);
331 mvaInput_[14] = std::copysign(+1.
f, (
float)tauLifetimeInfo.
dxy());
334 mvaInput_[17] = std::copysign(+1.
f, (
float)tauLifetimeInfo.
ip3d());
342 mvaInput_[1] =
std::abs((
float)tau->eta());
347 mvaInput_[6] = tauDecayMode;
349 mvaInput_[8] =
std::min(0.5
f, ptWeightedDetaStrip);
350 mvaInput_[9] =
std::min(0.5
f, ptWeightedDphiStrip);
351 mvaInput_[10] =
std::min(0.5
f, ptWeightedDrSignal);
352 mvaInput_[11] =
std::min(0.5
f, ptWeightedDrIsolation);
354 mvaInput_[13] = std::copysign(+1.
f, (
float)tauLifetimeInfo.
dxy());
357 mvaInput_[16] = std::copysign(+1.
f, (
float)tauLifetimeInfo.
ip3d());
363 mvaInput_[22] =
std::max(-1.
f, gjAngleDiff);
366 double mvaValue = mvaReader_->GetClassifier(mvaInput_);
368 edm::LogPrint(
"PFTauDiscByMVAIsol2") <<
"<PFRecoTauDiscriminationByMVAIsolationRun2::discriminate>:";
369 edm::LogPrint(
"PFTauDiscByMVAIsol2") <<
" tau: Pt = " << tau->pt() <<
", eta = " << tau->eta();
370 edm::LogPrint(
"PFTauDiscByMVAIsol2") <<
" isolation: charged = " << chargedIsoPtSum
371 <<
", neutral = " << neutralIsoPtSum <<
", PUcorr = " << puCorrPtSum;
372 edm::LogPrint(
"PFTauDiscByMVAIsol2") <<
" decay mode = " << tauDecayMode;
373 edm::LogPrint(
"PFTauDiscByMVAIsol2") <<
" impact parameter: distance = " << tauLifetimeInfo.
dxy()
374 <<
", significance = " << tauLifetimeInfo.
dxy_Sig();
377 <<
" distance = " << decayDistMag <<
", significance = " << tauLifetimeInfo.
flightLengthSig();
378 edm::LogPrint(
"PFTauDiscByMVAIsol2") <<
"--> mvaValue = " << mvaValue;
390 desc.
add<
bool>(
"loadMVAfromDB",
true);
398 desc.
add<
int>(
"verbosity", 0);
400 fillProducerDescriptions(desc);
402 descriptions.
add(
"pfRecoTauDiscriminationByMVAIsolationRun2", desc);
edm::EDGetTokenT< reco::TauDiscriminatorContainer > basicTauDiscriminators_token
unsigned int n_photons_total(const reco::PFTau &tau)
return total number of pf photon candidates with pT>500 MeV, which are associated to signal ...
~PFRecoTauDiscriminationByMVAIsolationRun2() override
bool hasSecondaryVertex() const
static std::vector< std::string > checklist log
ParameterDescriptionBase * addOptional(U const &iLabel, T const &value)
const GBRForest * mvaReader_
edm::Handle< PFTauTIPAssociationByRef > tauLifetimeInfos
float pt_weighted_dr_signal(const reco::PFTau &tau, int dm)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
#define DEFINE_FWK_MODULE(type)
StableProvenance const & stable() const
PFRecoTauDiscriminationByMVAIsolationRun2(const edm::ParameterSet &cfg)
bool getData(T &iHolder) const
double flightLengthSig() const
edm::Handle< TauCollection > taus_
const Vector & flightLength() const
float pt_weighted_deta_strip(const reco::PFTau &tau, int dm)
ParameterSet const & parameterSet(StableProvenance const &provenance, ProcessHistory const &history)
Abs< T >::type abs(const T &t)
float pt_weighted_dr_iso(const reco::PFTau &tau, int dm)
LocationCode location() const
Where was the file found?
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Log< level::Warning, true > LogPrint
bool isNull() const
Checks for null.
edm::AssociationVector< reco::PFTauRefProd, std::vector< reco::PFTauTransverseImpactParameterRef > > PFTauTIPAssociationByRef
edm::EDGetTokenT< PFTauTIPAssociationByRef > TauTransverseImpactParameters_token
ProcessHistory const & processHistory() const override
edm::Handle< reco::TauDiscriminatorContainer > basicTauDiscriminators_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
float pt_weighted_dphi_strip(const reco::PFTau &tau, int dm)
float eratio(const reco::PFTau &tau)
return ratio of energy in ECAL over sum of energy in ECAL and HCAL
T getParameter(std::string const &) const
void add(std::string const &label, ParameterSetDescription const &psetDescription)
ProcessHistoryID const & processHistoryID() const
edm::ProcessHistoryID phID_
std::string fullPath() const
float lead_track_chi2(const reco::PFTau &tau)
return chi2 of the leading track ==> deprecated? <==
std::string input_id_name_suffix_
edm::FileInPath inputFileName_
edm::ESGetToken< GBRForest, GBRWrapperRcd > mvaToken_
Power< A, B >::type pow(const A &a, const B &b)
std::vector< TFile * > inputFilesToDelete_
std::vector< float > rawValues