10 #include "TMatrixDSym.h"
11 #include "TMatrixDSymEigen.h"
26 for (
int i0 = 0; i0 <
NWPs; i0++) {
32 for (
int i1 = 0; i1 < 1; i1++) {
33 std::vector<double> pt010 = jetConfig.
getParameter<std::vector<double> >((
"Pt010_" + lCutType).c_str());
34 std::vector<double> pt1020 = jetConfig.
getParameter<std::vector<double> >((
"Pt1020_" + lCutType).c_str());
35 std::vector<double> pt2030 = jetConfig.
getParameter<std::vector<double> >((
"Pt2030_" + lCutType).c_str());
36 std::vector<double> pt3050 = jetConfig.
getParameter<std::vector<double> >((
"Pt3050_" + lCutType).c_str());
37 for (
int i2 = 0; i2 <
NPts; i2++)
39 for (
int i2 = 0; i2 <
NPts; i2++)
40 mvacut_[i0][1][i2] = pt1020[i2];
41 for (
int i2 = 0; i2 <
NPts; i2++)
42 mvacut_[i0][2][i2] = pt2030[i2];
43 for (
int i2 = 0; i2 <
NPts; i2++)
44 mvacut_[i0][3][i2] = pt3050[i2];
53 Float_t impactParTkThreshod,
111 void Assign(
const std::vector<float> &vec,
float &
a,
float &
b,
float &
c,
float &
d) {
112 size_t sz = vec.size();
113 a = (sz > 0 ? vec[0] : 0.);
114 b = (sz > 1 ? vec[1] : 0.);
115 c = (sz > 2 ? vec[2] : 0.);
116 d = (sz > 3 ? vec[3] : 0.);
126 reader_ =
new TMVA::Reader(
"!Color:Silent");
158 if (jetPt > 10 && jetPt < 20)
160 if (jetPt >= 20 && jetPt < 30)
166 if (fabs(jetEta) > 2.5 && fabs(jetEta) < 2.75)
168 if (fabs(jetEta) >= 2.75 && fabs(jetEta) < 3.0)
170 if (fabs(jetEta) >= 3.0 && fabs(jetEta) < 5.0)
172 return std::pair<int, int>(ptId, etaId);
176 std::pair<int, int> jetIdKey =
getJetIdKey(jetPt, jetEta);
202 typedef std::vector<reco::PFCandidatePtr> constituents_type;
203 typedef std::vector<reco::PFCandidatePtr>::iterator constituents_iterator;
216 std::vector<float>
frac, fracCh, fracEm, fracNeut;
217 constexpr
int ncones = 4;
218 std::array<float, ncones> cones{{0.1, 0.2, 0.3, 0.4}};
219 std::array<float *, ncones> coneFracs{
221 TMatrixDSym covMatrix(2);
226 float sumPt = 0., sumPt2 = 0., sumTkPt = 0., sumPtCh = 0, sumPtNe = 0;
230 float sum_detadphi = 0;
236 for (constituents_iterator it = constituents.begin(); it != constituents.end(); ++it) {
238 float candPt = icand->pt();
239 float candPtFrac = candPt /
jetPt;
241 float candDeta = fabs((*it)->eta() - jet->
eta());
243 float candPtDr = candPt * candDr;
244 size_t icone =
std::lower_bound(&cones[0], &cones[ncones], candDr) - &cones[0];
245 float weight2 = candPt * candPt;
247 if (lLead.
isNull() || candPt > lLead->pt()) {
250 }
else if ((lSecond.
isNull() || candPt > lSecond->pt()) && (candPt < lLead->
pt())) {
259 sumPt2 += candPt * candPt;
260 sum_deta += candDeta * weight2;
261 sum_dphi += candDphi * weight2;
262 sum_deta2 += candDeta * candDeta * weight2;
263 sum_detadphi += candDeta * candDphi * weight2;
264 sum_dphi2 += candDphi * candDphi * weight2;
268 frac.push_back(candPtFrac);
269 if (icone < ncones) {
270 *coneFracs[icone] += candPt;
274 if (lLeadNeut.
isNull() || candPt > lLeadNeut->pt()) {
278 fracNeut.push_back(candPtFrac);
283 if (lLeadEm.
isNull() || candPt > lLeadEm->pt()) {
287 fracEm.push_back(candPtFrac);
292 if (lLeadCh.
isNull() || candPt > lLeadCh->pt()) {
296 fracCh.push_back(candPtFrac);
301 float tkpt = icand->trackRef()->pt();
305 bool inAnyOther =
false;
307 double dZ0 = fabs(icand->trackRef()->dz(vtx->
position()));
309 for (reco::VertexCollection::const_iterator vi = allvtx.begin(); vi != allvtx.end(); ++vi) {
317 if (!isVtx0 && !inAnyOther) {
324 if (inVtx0 && !inAnyOther) {
326 }
else if (!inVtx0 && inAnyOther) {
332 }
else if (dZ < 0.2) {
337 if (lTrail.
isNull() || candPt < lTrail->pt()) {
355 impactTrack = lLeadCh->trackRef();
374 std::sort(frac.begin(), frac.end(), std::greater<float>());
375 std::sort(fracCh.begin(), fracCh.end(), std::greater<float>());
376 std::sort(fracEm.begin(), fracEm.end(), std::greater<float>());
377 std::sort(fracNeut.begin(), fracNeut.end(), std::greater<float>());
401 for (
size_t ic = 0; ic < ncones; ++ic) {
402 *coneFracs[ic] /=
jetPt;
407 for (
unsigned int i0 = 0; i0 < frac.size(); i0++) {
408 ptRMS += (frac[i0] - ptMean) * (frac[i0] - ptMean);
428 float ave_deta = sum_deta / sumPt2;
429 float ave_dphi = sum_dphi / sumPt2;
430 float ave_deta2 = sum_deta2 / sumPt2;
431 float ave_dphi2 = sum_dphi2 / sumPt2;
432 float a = ave_deta2 - ave_deta * ave_deta;
433 float b = ave_dphi2 - ave_dphi * ave_dphi;
434 float c = -(sum_detadphi / sumPt2 - ave_deta * ave_dphi);
437 if ((((a - b) * (a - b) + 4 * c * c)) > 0) {
438 float delta =
sqrt(((a - b) * (a - b) + 4 * c * c));
439 if (a + b + delta > 0) {
440 axis1 =
sqrt(0.5 * (a + b + delta));
442 if (a + b - delta > 0) {
443 axis2 =
sqrt(0.5 * (a + b - delta));
453 float ddetaR_sum(0.0), ddphiR_sum(0.0);
456 float weight = part->pt() * part->pt();
457 float deta = part->eta() - jet->
eta();
459 float ddeta, ddphi, ddR;
460 ddeta = deta - ave_deta;
462 ddR =
sqrt(ddeta * ddeta + ddphi * ddphi);
463 ddetaR_sum += ddR * ddeta *
weight;
464 ddphiR_sum += ddR * ddphi *
weight;
467 float ddetaR_ave = ddetaR_sum / sumPt2;
468 float ddphiR_ave = ddphiR_sum / sumPt2;
469 internalId_.dRMean_ =
sqrt(ddetaR_ave * ddetaR_ave + ddphiR_ave * ddphiR_ave);
481 std::stringstream
out;
483 out << std::setw(15) << it->first << std::setw(3) <<
"=" << std::setw(5) << *it->second.first <<
" ("
484 << std::setw(5) << it->second.second <<
")" << std::endl;
492 *it->second.first = it->second.second;
496 #define INIT_VARIABLE(NAME, TMVANAME, VAL) \
497 internalId_.NAME##_ = VAL; \
498 variables_[#NAME] = std::make_pair(&internalId_.NAME##_, VAL);
constexpr double deltaPhi(double phi1, double phi2)
const edm::EventSetup & c
bool isNonnull() const
Checks for non-null.
std::vector< std::string > tmvaVariables_
float chargedEmEnergy() const
chargedEmEnergy
double pt() const final
transverse momentum
uint16_t *__restrict__ id
Float_t mvacut_[NWPs][NEtas][NPts]
#define INIT_VARIABLE(NAME, TMVANAME, VAL)
virtual double pt() const =0
transverse momentum
Base class for all types of Jets.
std::vector< std::string > tmvaSpectators_
std::vector< Vertex > VertexCollection
collection of Vertex objects
variables_list_t variables_
static constexpr int NPts
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
int chargedMultiplicity() const
chargedMultiplicity
const Point & position() const
position
Jets made from PFObjects.
float neutralEmEnergy() const
neutralEmEnergy
virtual reco::PFCandidatePtr getPFConstituent(unsigned fIndex) const
get specific constituent
edm::RefToBase< reco::Track > TrackBaseRef
persistent reference to a Track, using views
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
bool isNull() const
Checks for null.
PileupJetIdentifier computeIdVariables(const reco::Jet *jet, float jec, const reco::Vertex *, const reco::VertexCollection &, double rho, bool calculateMva=false)
trackRef_iterator tracks_end() const
last iterator over tracks
void SetPtEtaPhi(const reco::Candidate &p, float &pt, float &eta, float &phi)
trackRef_iterator tracks_begin() const
first iterator over tracks
int neutralMultiplicity() const
neutralMultiplicity
bool isNonnull() const
Checks for non-null.
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
static constexpr float d0
T getParameter(std::string const &) const
PileupJetIdentifier internalId_
TMVA::IMethod * loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)
void set(const PileupJetIdentifier &)
double mass() const final
mass
__host__ __device__ constexpr RandomIt lower_bound(RandomIt first, RandomIt last, const T &value, Compare comp={})
std::map< std::string, std::string > tmvaNames_
std::pair< int, int > getJetIdKey(float jetPt, float jetEta)
PileupJetIdentifier computeMva()
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
std::string fullPath() const
float neutralHadronEnergy() const
neutralHadronEnergy
Float_t impactParTkThreshod_
MVAJetPuId(int version=CATEv0, const std::string &tmvaWeight="", const std::string &tmvaMethod="", Float_t impactParTkThreshod_=1., const std::vector< std::string > &tmvaVariables=std::vector< std::string >())
int computeIDflag(float mva, float jetPt, float jetEta)
std::string dumpVariables() const
float chargedHadronEnergy() const
chargedHadronEnergy
void Assign(const std::vector< float > &vec, float &a, float &b, float &c, float &d)
static constexpr int NWPs
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity
double energy() const final
energy
double eta() const final
momentum pseudorapidity