12 #include "TMatrixDSym.h"
13 #include "TMatrixDSymEigen.h"
23 : cutBased_(ps.getParameter<
bool>(
"cutBased")),
24 etaBinnedWeights_(
false),
27 label_(ps.getParameter<
std::
string>(
"label")),
31 std::vector<edm::FileInPath> tmvaEtaWeights;
36 etaBinnedWeights_ = ps.getParameter<
bool>(
"etaBinnedWeights");
37 if (etaBinnedWeights_) {
38 const std::vector<edm::ParameterSet>&
trainings = ps.getParameter<std::vector<edm::ParameterSet> >(
"trainings");
39 nEtaBins_ = ps.getParameter<
int>(
"nEtaBins");
40 for (
int v = 0;
v < nEtaBins_;
v++) {
42 jEtaMin_.push_back(
trainings.at(
v).getParameter<
double>(
"jEtaMin"));
43 jEtaMax_.push_back(
trainings.at(
v).getParameter<
double>(
"jEtaMax"));
45 for (
int v = 0;
v < nEtaBins_;
v++) {
46 tmvaEtaVariables_.push_back(
trainings.at(
v).getParameter<std::vector<std::string> >(
"tmvaVariables"));
49 tmvaVariables_ = ps.getParameter<std::vector<std::string> >(
"tmvaVariables");
51 tmvaMethod_ = ps.getParameter<
std::string>(
"tmvaMethod");
52 tmvaSpectators = ps.getParameter<std::vector<std::string> >(
"tmvaSpectators");
53 version = ps.getParameter<
int>(
"version");
59 for (
int i0 = 0; i0 < 3; i0++) {
68 for (
int i1 = 0;
i1 < nCut;
i1++) {
70 if (cutBased_ &&
i1 == 0)
71 lFullCutType =
"BetaStar" + lCutType;
72 if (cutBased_ &&
i1 == 1)
73 lFullCutType =
"RMS" + lCutType;
74 std::vector<double> pt010 = jetConfig.
getParameter<std::vector<double> >((
"Pt010_" + lFullCutType).c_str());
75 std::vector<double> pt1020 = jetConfig.
getParameter<std::vector<double> >((
"Pt1020_" + lFullCutType).c_str());
76 std::vector<double> pt2030 = jetConfig.
getParameter<std::vector<double> >((
"Pt2030_" + lFullCutType).c_str());
77 std::vector<double> pt3040 = jetConfig.
getParameter<std::vector<double> >((
"Pt3040_" + lFullCutType).c_str());
78 std::vector<double> pt4050 = jetConfig.
getParameter<std::vector<double> >((
"Pt4050_" + lFullCutType).c_str());
80 for (
int i2 = 0;
i2 < 4;
i2++)
81 mvacut_[i0][0][
i2] = pt010[
i2];
82 for (
int i2 = 0;
i2 < 4;
i2++)
83 mvacut_[i0][1][
i2] = pt1020[
i2];
84 for (
int i2 = 0;
i2 < 4;
i2++)
85 mvacut_[i0][2][
i2] = pt2030[
i2];
86 for (
int i2 = 0;
i2 < 4;
i2++)
87 mvacut_[i0][3][
i2] = pt3040[
i2];
88 for (
int i2 = 0;
i2 < 4;
i2++)
89 mvacut_[i0][4][
i2] = pt4050[
i2];
91 if (cutBased_ &&
i1 == 0) {
92 for (
int i2 = 0;
i2 < 4;
i2++)
93 betaStarCut_[i0][0][
i2] = pt010[
i2];
94 for (
int i2 = 0;
i2 < 4;
i2++)
95 betaStarCut_[i0][1][
i2] = pt1020[
i2];
96 for (
int i2 = 0;
i2 < 4;
i2++)
97 betaStarCut_[i0][2][
i2] = pt2030[
i2];
98 for (
int i2 = 0;
i2 < 4;
i2++)
99 betaStarCut_[i0][3][
i2] = pt3040[
i2];
100 for (
int i2 = 0;
i2 < 4;
i2++)
101 betaStarCut_[i0][4][
i2] = pt4050[
i2];
103 if (cutBased_ &&
i1 == 1) {
104 for (
int i2 = 0;
i2 < 4;
i2++)
105 rmsCut_[i0][0][
i2] = pt010[
i2];
106 for (
int i2 = 0;
i2 < 4;
i2++)
107 rmsCut_[i0][1][
i2] = pt1020[
i2];
108 for (
int i2 = 0;
i2 < 4;
i2++)
109 rmsCut_[i0][2][
i2] = pt2030[
i2];
110 for (
int i2 = 0;
i2 < 4;
i2++)
111 rmsCut_[i0][3][
i2] = pt3040[
i2];
112 for (
int i2 = 0;
i2 < 4;
i2++)
113 rmsCut_[i0][4][
i2] = pt4050[
i2];
119 assert(tmvaMethod_.empty() || ((!tmvaVariables_.empty() || (!tmvaEtaVariables_.empty())) &&
version ==
USER));
122 if ((!cutBased_) && (runMvas_)) {
123 if (etaBinnedWeights_) {
124 for (
int v = 0;
v < nEtaBins_;
v++) {
139 void assign(
const std::vector<float>& vec,
float&
a,
float&
b,
float&
c,
float&
d) {
140 size_t sz = vec.size();
141 a = (sz > 0 ? vec[0] : 0.);
142 b = (sz > 1 ? vec[1] : 0.);
143 c = (sz > 2 ? vec[2] : 0.);
144 d = (sz > 3 ? vec[3] : 0.);
159 const std::unique_ptr<const GBRForest>&
reader) {
160 std::vector<float>
vars;
161 for (std::vector<std::string>::const_iterator it = varList.begin(); it != varList.end(); ++it) {
216 return std::pair<int, int>(ptId, etaId);
221 float betaStarModified = betaStarClassic /
log(nvtx - 0.64);
271 const pat::Jet* patjet = dynamic_cast<const pat::Jet*>(
jet);
273 assert(patjet !=
nullptr || pfjet !=
nullptr);
274 if (patjet !=
nullptr &&
jec == 0.) {
281 const reco::Candidate *lLead =
nullptr, *lSecond =
nullptr, *lLeadNeut =
nullptr, *lLeadEm =
nullptr,
282 *lLeadCh =
nullptr, *lTrail =
nullptr;
283 std::vector<float>
frac, fracCh, fracEm, fracNeut;
284 float cones[] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7};
285 size_t ncones =
sizeof(cones) /
sizeof(
float);
314 TMatrixDSym covMatrix(2);
317 float sumPt = 0., sumPt2 = 0., sumTkPt = 0., sumPtCh = 0, sumPtNe = 0;
318 float multNeut = 0.0;
327 for (
unsigned i = 0;
i <
jet->numberOfSourceCandidatePtrs(); ++
i) {
333 bool isPacked =
true;
334 if (lPack ==
nullptr) {
337 float candPuppiWeight = 1.0;
340 float candPt = (icand->
pt()) * candPuppiWeight;
341 float candPtFrac = candPt /
jetPt;
343 float candDeta = icand->
eta() -
jet->eta();
345 float candPtDr = candPt * candDr;
346 size_t icone =
std::lower_bound(&cones[0], &cones[ncones], candDr) - &cones[0];
352 if (lLead ==
nullptr || candPt > lLead->
pt()) {
355 }
else if ((lSecond ==
nullptr || candPt > lSecond->pt()) && (candPt < lLead->
pt())) {
362 covMatrix(0, 0) += candPt * candPt * candDeta * candDeta;
363 covMatrix(0, 1) += candPt * candPt * candDeta * candDphi;
364 covMatrix(1, 1) += candPt * candPt * candDphi * candDphi;
367 sumPt2 += candPt * candPt;
370 frac.push_back(candPtFrac);
372 if (icone < ncones) {
373 *coneFracs[icone] += candPt;
378 if (lLeadNeut ==
nullptr || candPt > lLeadNeut->pt()) {
382 fracNeut.push_back(candPtFrac);
383 if (icone < ncones) {
384 *coneNeutFracs[icone] += candPt;
388 multNeut += candPuppiWeight;
391 if (icand->
pdgId() == 22) {
392 if (lLeadEm ==
nullptr || candPt > lLeadEm->pt()) {
396 fracEm.push_back(candPtFrac);
397 if (icone < ncones) {
398 *coneEmFracs[icone] += candPt;
402 multNeut += candPuppiWeight;
405 multNeut += candPuppiWeight;
408 if (icand->
charge() != 0) {
409 if (lLeadCh ==
nullptr || candPt > lLeadCh->pt()) {
413 if (lPF &&
std::abs(icand->
pdgId()) == 13 && pfTrk ==
nullptr) {
419 <<
"Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack ";
422 if (pfTrk ==
nullptr) {
427 }
else if (lPF !=
nullptr) {
439 fracCh.push_back(candPtFrac);
440 if (icone < ncones) {
441 *coneChFracs[icone] += candPt;
446 if (icand->
charge() != 0) {
452 bool inVtx0 =
vtx->trackWeight(lPF->
trackRef()) > 0;
454 bool inAnyOther =
false;
458 for (reco::VertexCollection::const_iterator vi = allvtx.begin(); vi != allvtx.end(); ++vi) {
460 if (
iv.isFake() ||
iv.ndof() < 4) {
464 bool isVtx0 = (
iv.position() -
vtx->position()).
r() < 0.02;
467 if (!isVtx0 && !inAnyOther) {
468 inAnyOther =
vtx->trackWeight(lPF->
trackRef()) <= 0;
474 if (inVtx0 && !inAnyOther) {
476 }
else if (!inVtx0 && inAnyOther) {
482 }
else if (dZ < 0.2) {
490 bool inVtxOther =
false;
492 double dZ_tmp = 9999.;
493 for (
unsigned vtx_i = 0; vtx_i < allvtx.size(); vtx_i++) {
494 auto iv = allvtx[vtx_i];
500 bool isVtx0 = (
iv.position() -
vtx->position()).
r() < 0.02;
505 if (lPack->
fromPV(vtx_i) == 0)
507 dZ0 = lPack->
dz(
iv.position());
510 if (fabs(lPack->
dz(
iv.position())) < fabs(dZ_tmp)) {
511 dZ_tmp = lPack->
dz(
iv.position());
516 }
else if (inVtxOther) {
519 if (fabs(dZ0) < 0.2) {
521 }
else if (fabs(dZ_tmp) < 0.2) {
528 if (lTrail ==
nullptr || candPt < lTrail->
pt()) {
534 assert(!(lLead ==
nullptr));
536 if (lSecond ==
nullptr) {
539 if (lLeadNeut ==
nullptr) {
542 if (lLeadEm ==
nullptr) {
545 if (lLeadCh ==
nullptr) {
549 if (patjet !=
nullptr) {
570 float sum_deta(0.0), sum_dphi(0.0);
571 float ave_deta(0.0), ave_dphi(0.0);
572 for (
size_t j = 0;
j <
jet->numberOfDaughters();
j++) {
573 const auto&
part =
jet->daughterPtr(
j);
574 if (!(
part.isAvailable() &&
part.isNonnull())) {
578 float partPuppiWeight = 1.0;
581 if (partpack !=
nullptr)
585 float weight = (
part->pt()) * partPuppiWeight;
588 float deta =
part->eta() -
jet->eta();
590 sum_deta += deta * weight2;
591 sum_dphi += dphi * weight2;
593 ave_deta = sum_deta / sumW2;
594 ave_dphi = sum_dphi / sumW2;
598 float ddetaR_sum(0.0), ddphiR_sum(0.0), pull_tmp(0.0);
599 for (
size_t i = 0;
i <
jet->numberOfDaughters();
i++) {
600 const auto&
part =
jet->daughterPtr(
i);
601 if (!(
part.isAvailable() &&
part.isNonnull())) {
605 float partPuppiWeight = 1.0;
608 if (partpack !=
nullptr)
612 float weight = partPuppiWeight * (
part->pt()) * partPuppiWeight * (
part->pt());
613 float deta =
part->eta() -
jet->eta();
615 float ddeta, ddphi, ddR;
616 ddeta = deta - ave_deta;
617 ddphi = dphi - ave_dphi;
618 ddR =
sqrt(ddeta * ddeta + ddphi * ddphi);
619 ddetaR_sum += ddR * ddeta *
weight;
620 ddphiR_sum += ddR * ddphi *
weight;
623 float ddetaR_ave = ddetaR_sum / sumW2;
624 float ddphiR_ave = ddphiR_sum / sumW2;
625 pull_tmp =
sqrt(ddetaR_ave * ddetaR_ave + ddphiR_ave * ddphiR_ave);
637 std::sort(fracCh.begin(), fracCh.end(), std::greater<float>());
638 std::sort(fracEm.begin(), fracEm.end(), std::greater<float>());
639 std::sort(fracNeut.begin(), fracNeut.end(), std::greater<float>());
651 covMatrix(0, 0) /= sumPt2;
652 covMatrix(0, 1) /= sumPt2;
653 covMatrix(1, 1) /= sumPt2;
654 covMatrix(1, 0) = covMatrix(0, 1);
659 eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
667 if (lSecond ==
nullptr) {
676 for (
size_t ic = 0; ic < ncones; ++ic) {
677 *coneFracs[ic] /=
jetPt;
678 *coneEmFracs[ic] /=
jetPt;
679 *coneNeutFracs[ic] /=
jetPt;
680 *coneChFracs[ic] /=
jetPt;
685 for (
unsigned int i0 = 0; i0 <
frac.size(); i0++) {
686 ptRMS += (
frac[i0] - ptMean) * (
frac[i0] - ptMean);
724 std::stringstream
out;
726 out << std::setw(15) << it->first << std::setw(3) <<
"=" << std::setw(5) << *it->second.first <<
" ("
727 << std::setw(5) << it->second.second <<
")" << std::endl;
736 *it->second.first = it->second.second;
741 #define INIT_VARIABLE(NAME, TMVANAME, VAL) \
742 internalId_.NAME##_ = VAL; \
743 variables_[#NAME] = std::make_pair(&internalId_.NAME##_, VAL);