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

Preparations for future trainings of the MVA tauID #18430

Merged
merged 11 commits into from
Apr 29, 2017
+208 −118
Merged

Preparations for future trainings of the MVA tauID #18430

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
db62270
- Bug fix: new tau decay mode (3-prong+1Pi0, 11) is considered when r…
anehrkor Apr 6, 2017
bfcbe20
Small fix: decayDistMag contains already square root
anehrkor Apr 6, 2017
5748327
Merge pull request #43 from anehrkor/CMSSW_9_1_X_tau-pog_tauID-GJangle
roger-wolf Apr 20, 2017
0526fdc
Merge branch 'CMSSW_9_1_X_tau-pog' into CMSSW_9_1_X_tau-pog_tauID-GJa…
roger-wolf Apr 20, 2017
c4c2d0e
Implement Slava's requests: remove dependence on functions from TMath…
anehrkor Apr 25, 2017
d14f66a
Merge pull request #46 from anehrkor/CMSSW_9_1_X_tau-pog_tauID-GJangle
roger-wolf Apr 25, 2017
f947030
Add protection against arguments of asin and acos being outside of al…
anehrkor Apr 25, 2017
1a3f320
Merge pull request #47 from anehrkor/CMSSW_9_1_X_tau-pog_tauID-GJangle
roger-wolf Apr 25, 2017
2076be9
Replace std::fabs by std::abs; remove redundant std::abs statements
anehrkor Apr 26, 2017
6bfb858
Add a few line breaks to improve readability
anehrkor Apr 26, 2017
5be50cb
Merge pull request #48 from anehrkor/CMSSW_9_1_X_tau-pog_tauID-GJangle
roger-wolf Apr 28, 2017
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
46 changes: 42 additions & 4 deletions RecoTauTag/RecoTau/plugins/PATTauDiscriminationByMVAIsolationRun2.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,13 +97,16 @@ class PATTauDiscriminationByMVAIsolationRun2 : public PATTauDiscriminationProduc
else if ( mvaOpt_string == "DBnewDMwLT" ) mvaOpt_ = kDBnewDMwLT;
else if ( mvaOpt_string == "PWoldDMwLT" ) mvaOpt_ = kPWoldDMwLT;
else if ( mvaOpt_string == "PWnewDMwLT" ) mvaOpt_ = kPWnewDMwLT;
else if ( mvaOpt_string == "DBoldDMwLTwGJ" ) mvaOpt_ = kDBoldDMwLTwGJ;
else if ( mvaOpt_string == "DBnewDMwLTwGJ" ) mvaOpt_ = kDBnewDMwLTwGJ;
else throw cms::Exception("PATTauDiscriminationByMVAIsolationRun2")
<< " Invalid Configuration Parameter 'mvaOpt' = " << mvaOpt_string << " !!\n";

if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) mvaInput_ = new float[6];
else if ( mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT ) mvaInput_ = new float[12];
else if ( mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kDBnewDMwLT ||
mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kPWnewDMwLT) mvaInput_ = new float[23];
mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kPWnewDMwLT ||
mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ) mvaInput_ = new float[23];
else assert(0);

chargedIsoPtSums_ = cfg.getParameter<std::string>("srcChargedIsoPtSum");
Expand Down Expand Up @@ -142,7 +145,7 @@ class PATTauDiscriminationByMVAIsolationRun2 : public PATTauDiscriminationProduc
bool loadMVAfromDB_;
edm::FileInPath inputFileName_;
const GBRForest* mvaReader_;
enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT };
enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT, kDBoldDMwLTwGJ, kDBnewDMwLTwGJ };
int mvaOpt_;
float* mvaInput_;

Expand Down Expand Up @@ -185,8 +188,8 @@ double PATTauDiscriminationByMVAIsolationRun2::discriminate(const TauRef& tau) c

int tauDecayMode = tau->decayMode();

if ( ((mvaOpt_ == kOldDMwoLT || mvaOpt_ == kOldDMwLT || mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kPWoldDMwLT) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
((mvaOpt_ == kNewDMwoLT || mvaOpt_ == kNewDMwLT || mvaOpt_ == kDBnewDMwLT || mvaOpt_ == kPWnewDMwLT) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 || tauDecayMode == 10)) ) {
if ( ((mvaOpt_ == kOldDMwoLT || mvaOpt_ == kOldDMwLT || mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kDBoldDMwLTwGJ) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
((mvaOpt_ == kNewDMwoLT || mvaOpt_ == kNewDMwLT || mvaOpt_ == kDBnewDMwLT || mvaOpt_ == kPWnewDMwLT || mvaOpt_ == kDBnewDMwLTwGJ) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 || tauDecayMode == 10 || tauDecayMode == 11)) ) {
Copy link
Contributor

Choose a reason for hiding this comment

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

together with the asin update, it would be nice to add a few linebreaks here to improve readability.


float chargedIsoPtSum = tau->tauID(chargedIsoPtSums_);
float neutralIsoPtSum = tau->tauID(neutralIsoPtSums_);
Expand All @@ -209,6 +212,17 @@ double PATTauDiscriminationByMVAIsolationRun2::discriminate(const TauRef& tau) c
// ---
float leadingTrackChi2 = tau->leadingTrackNormChi2();
float eRatio = clusterVariables_.tau_Eratio(*tau);

// Difference between measured and maximally allowed Gottfried-Jackson angle
float gjAngleDiff = -999;
if ( tauDecayMode == 10 ) {
double mTau = 1.77682;
double mAOne = tau->p4().M();
double pAOneMag = tau->p();
double thetaGJmax = TMath::ASin( (TMath::Power(mTau,2) - TMath::Power(mAOne,2) ) / ( 2 * mTau * pAOneMag ) );
double thetaGJmeasured = TMath::ACos( ( tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ ) / ( pAOneMag * decayDistMag ) );
Copy link
Contributor

Choose a reason for hiding this comment

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

please use native c++ methods where available.
std::pow, std::acos, std::asin etc)

Copy link
Contributor

Choose a reason for hiding this comment

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

I'm missing what by construction restricts the argument of asin to have an absolute value less or equal 1.0?
It looks like there should be a protection.

Copy link
Contributor

Choose a reason for hiding this comment

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

Dear Slava,

the valua of mAOne can only be between 0.8 and 1.5 GeV. This cut is placed on 3-prong taus during the tau reconstruction by HPS. Thus, there is no need for protection.

Cheers,
Alex

Copy link
Contributor

Choose a reason for hiding this comment

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

asin( (1.77682**2 - 0.8**2)/2/1.77682/pAOneMag)) ~ asin(0.71/pAOneMag)
Does it mean that tau->p() can not possibly be below 0.71 ?

Copy link
Contributor

Choose a reason for hiding this comment

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

For that to happen, the transverse momentum of the tau would have to be below 0.71 GeV. In MiniAOD, only taus with pt greater than 18 GeV are saved.
However, in AOD there is no such cut and the minimum pt of the tracks used to reconstruct the tau is 0.5 GeV, I believe. So there it might be necessary.
I will add a statement requiring tau->p() to be greater than 1 although in typical analyses, this will never not be true since typically taus with pt greater than 20 GeV are required.

gjAngleDiff = thetaGJmeasured - thetaGJmax;
}

if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) {
mvaInput_[0] = std::log(std::max((float)1., (float)tau->pt()));
Expand Down Expand Up @@ -278,6 +292,30 @@ double PATTauDiscriminationByMVAIsolationRun2::discriminate(const TauRef& tau) c
mvaInput_[20] = ( tau->hasSecondaryVertex() ) ? 1. : 0.;
mvaInput_[21] = std::sqrt(decayDistMag);
mvaInput_[22] = std::min((float)10., tau->flightLengthSig());
} else if ( mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ ) {
mvaInput_[0] = std::log(std::max((float)1., (float)tau->pt()));
Copy link
Contributor

Choose a reason for hiding this comment

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

1.f literal is shorter than (float)1.
same applies to all values below

mvaInput_[1] = std::fabs(tau->eta());
mvaInput_[2] = std::log(std::max((float)1.e-2, chargedIsoPtSum));
mvaInput_[3] = std::log(std::max((float)1.e-2, neutralIsoPtSum));
mvaInput_[4] = std::log(std::max((float)1.e-2, puCorrPtSum));
mvaInput_[5] = std::log(std::max((float)1.e-2, photonPtSumOutsideSignalCone));
mvaInput_[6] = tauDecayMode;
mvaInput_[7] = std::min((float)30., nPhoton);
mvaInput_[8] = std::min((float)0.5, ptWeightedDetaStrip);
mvaInput_[9] = std::min((float)0.5, ptWeightedDphiStrip);
mvaInput_[10] = std::min((float)0.5, ptWeightedDrSignal);
mvaInput_[11] = std::min((float)0.5, ptWeightedDrIsolation);
mvaInput_[12] = std::min((float)1., eRatio);
mvaInput_[13] = TMath::Sign((float)+1., tau->dxy());
Copy link
Contributor

Choose a reason for hiding this comment

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

std::copysign(1.f, tau->dxy());

mvaInput_[14] = std::sqrt(std::min((float)1., std::fabs(tau->dxy())));
mvaInput_[15] = std::min((float)10., std::fabs(tau->dxy_Sig()));
mvaInput_[16] = TMath::Sign((float)+1., tau->ip3d());
mvaInput_[17] = std::sqrt(std::min((float)1., std::fabs(tau->ip3d())));
mvaInput_[18] = std::min((float)10., std::fabs(tau->ip3d_Sig()));
mvaInput_[19] = ( tau->hasSecondaryVertex() ) ? 1. : 0.;
mvaInput_[20] = std::sqrt(decayDistMag);
mvaInput_[21] = std::min((float)10., tau->flightLengthSig());
mvaInput_[22] = std::max((float)-1., gjAngleDiff);
}

double mvaValue = mvaReader_->GetClassifier(mvaInput_);
Expand Down
46 changes: 42 additions & 4 deletions RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByMVAIsolationRun2.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,13 +89,16 @@ class PFRecoTauDiscriminationByMVAIsolationRun2 : public PFTauDiscriminationProd
else if ( mvaOpt_string == "DBnewDMwLT" ) mvaOpt_ = kDBnewDMwLT;
else if ( mvaOpt_string == "PWoldDMwLT" ) mvaOpt_ = kPWoldDMwLT;
else if ( mvaOpt_string == "PWnewDMwLT" ) mvaOpt_ = kPWnewDMwLT;
else if ( mvaOpt_string == "DBoldDMwLTwGJ" ) mvaOpt_ = kDBoldDMwLTwGJ;
else if ( mvaOpt_string == "DBnewDMwLTwGJ" ) mvaOpt_ = kDBnewDMwLTwGJ;
else throw cms::Exception("PFRecoTauDiscriminationByMVAIsolationRun2")
<< " Invalid Configuration Parameter 'mvaOpt' = " << mvaOpt_string << " !!\n";

if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) mvaInput_ = new float[6];
else if ( mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT ) mvaInput_ = new float[12];
else if ( mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kDBnewDMwLT ||
mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kPWnewDMwLT) mvaInput_ = new float[23];
mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kPWnewDMwLT ||
mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ) mvaInput_ = new float[23];
else assert(0);

TauTransverseImpactParameters_token = consumes<PFTauTIPAssociationByRef>(cfg.getParameter<edm::InputTag>("srcTauTransverseImpactParameters"));
Expand Down Expand Up @@ -136,7 +139,7 @@ class PFRecoTauDiscriminationByMVAIsolationRun2 : public PFTauDiscriminationProd
bool loadMVAfromDB_;
edm::FileInPath inputFileName_;
const GBRForest* mvaReader_;
enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT };
enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT, kDBoldDMwLT, kDBnewDMwLT, kPWoldDMwLT, kPWnewDMwLT, kDBoldDMwLTwGJ, kDBnewDMwLTwGJ };
int mvaOpt_;
float* mvaInput_;

Expand Down Expand Up @@ -197,8 +200,8 @@ double PFRecoTauDiscriminationByMVAIsolationRun2::discriminate(const PFTauRef& t

int tauDecayMode = tau->decayMode();

if ( ((mvaOpt_ == kOldDMwoLT || mvaOpt_ == kOldDMwLT || mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kPWoldDMwLT) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
((mvaOpt_ == kNewDMwoLT || mvaOpt_ == kNewDMwLT || mvaOpt_ == kDBnewDMwLT || mvaOpt_ == kPWnewDMwLT) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 || tauDecayMode == 10)) ) {
if ( ((mvaOpt_ == kOldDMwoLT || mvaOpt_ == kOldDMwLT || mvaOpt_ == kDBoldDMwLT || mvaOpt_ == kPWoldDMwLT || mvaOpt_ == kDBoldDMwLTwGJ) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
((mvaOpt_ == kNewDMwoLT || mvaOpt_ == kNewDMwLT || mvaOpt_ == kDBnewDMwLT || mvaOpt_ == kPWnewDMwLT || mvaOpt_ == kDBnewDMwLTwGJ) && (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 || tauDecayMode == 10 || tauDecayMode == 11)) ) {
Copy link
Contributor

Choose a reason for hiding this comment

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

please add a couple of linebreaks here.

Copy link
Contributor

Choose a reason for hiding this comment

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

just to be sure: is the decayMode 11 important only in the "new" MVAs?

Copy link
Contributor

Choose a reason for hiding this comment

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

Correct. The decay mode 11 is only available in the "new decay modes" to which the "new" MVA trainings correspond. The "old" MVA trainings in turn correspond to the "old decay modes" (0, 1, 2 and 10).
I'll work on your other comments first thing tomorrow morning.


float chargedIsoPtSum = (*chargedIsoPtSums_)[tau];
float neutralIsoPtSum = (*neutralIsoPtSums_)[tau];
Expand All @@ -221,6 +224,17 @@ double PFRecoTauDiscriminationByMVAIsolationRun2::discriminate(const PFTauRef& t
float leadingTrackChi2 = clusterVariables_.tau_leadTrackChi2(*tau);
float eRatio = clusterVariables_.tau_Eratio(*tau);

// Difference between measured and maximally allowed Gottfried-Jackson angle
float gjAngleDiff = -999;
if ( tauDecayMode == 10 ) {
double mTau = 1.77682;
double mAOne = tau->p4().M();
double pAOneMag = tau->p();
double thetaGJmax = TMath::ASin( (TMath::Power(mTau,2) - TMath::Power(mAOne,2) ) / ( 2 * mTau * pAOneMag ) );
Copy link
Contributor

Choose a reason for hiding this comment

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

same comments as in PAT tau apply here.

double thetaGJmeasured = TMath::ACos( ( tau->p4().px() * decayDistX + tau->p4().py() * decayDistY + tau->p4().pz() * decayDistZ ) / ( pAOneMag * decayDistMag ) );
gjAngleDiff = thetaGJmeasured - thetaGJmax;
}

if ( mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT ) {
mvaInput_[0] = std::log(std::max((float)1., (float)tau->pt()));
mvaInput_[1] = std::fabs((float)tau->eta());
Expand Down Expand Up @@ -289,6 +303,30 @@ double PFRecoTauDiscriminationByMVAIsolationRun2::discriminate(const PFTauRef& t
mvaInput_[20] = ( tauLifetimeInfo.hasSecondaryVertex() ) ? 1. : 0.;
mvaInput_[21] = std::sqrt(decayDistMag);
mvaInput_[22] = std::min((float)10., (float)tauLifetimeInfo.flightLengthSig());
} else if ( mvaOpt_ == kDBoldDMwLTwGJ || mvaOpt_ == kDBnewDMwLTwGJ ) {
mvaInput_[0] = std::log(std::max((float)1., (float)tau->pt()));
mvaInput_[1] = std::fabs((float)tau->eta());
mvaInput_[2] = std::log(std::max((float)1.e-2, chargedIsoPtSum));
mvaInput_[3] = std::log(std::max((float)1.e-2, neutralIsoPtSum));
mvaInput_[4] = std::log(std::max((float)1.e-2, puCorrPtSum));
mvaInput_[5] = std::log(std::max((float)1.e-2, photonPtSumOutsideSignalCone));
mvaInput_[6] = tauDecayMode;
mvaInput_[7] = std::min((float)30., nPhoton);
mvaInput_[8] = std::min((float)0.5, ptWeightedDetaStrip);
mvaInput_[9] = std::min((float)0.5, ptWeightedDphiStrip);
mvaInput_[10] = std::min((float)0.5, ptWeightedDrSignal);
mvaInput_[11] = std::min((float)0.5, ptWeightedDrIsolation);
mvaInput_[12] = std::min((float)1., eRatio);
mvaInput_[13] = TMath::Sign((float)+1., (float)tauLifetimeInfo.dxy());
mvaInput_[14] = std::sqrt(std::fabs(std::min((float)1., std::fabs((float)tauLifetimeInfo.dxy()))));
mvaInput_[15] = std::min((float)10., std::fabs((float)tauLifetimeInfo.dxy_Sig()));
mvaInput_[16] = TMath::Sign((float)+1., (float)tauLifetimeInfo.ip3d());
mvaInput_[17] = std::sqrt(std::fabs(std::min((float)1., std::fabs((float)tauLifetimeInfo.ip3d()))));
mvaInput_[18] = std::min((float)10., std::fabs((float)tauLifetimeInfo.ip3d_Sig()));
mvaInput_[19] = ( tauLifetimeInfo.hasSecondaryVertex() ) ? 1. : 0.;
mvaInput_[20] = std::sqrt(decayDistMag);
mvaInput_[21] = std::min((float)10., (float)tauLifetimeInfo.flightLengthSig());
mvaInput_[22] = std::max((float)-1., gjAngleDiff);
}

double mvaValue = mvaReader_->GetClassifier(mvaInput_);
Expand Down