11 #include "TMatrixDSym.h"
12 #include "TMatrixDSymEigen.h"
50 for(
int i0 = 0; i0 < 3; i0++) {
56 for(
int i1 = 0; i1 < nCut; i1++) {
58 if(
cutBased_ && i1 == 0) lFullCutType =
"BetaStar"+ lCutType;
59 if(
cutBased_ && i1 == 1) lFullCutType =
"RMS" + lCutType;
60 std::vector<double> pt010 = jetConfig.
getParameter<std::vector<double> >((
"Pt010_" +lFullCutType).c_str());
61 std::vector<double> pt1020 = jetConfig.
getParameter<std::vector<double> >((
"Pt1020_"+lFullCutType).c_str());
62 std::vector<double> pt2030 = jetConfig.
getParameter<std::vector<double> >((
"Pt2030_"+lFullCutType).c_str());
63 std::vector<double> pt3050 = jetConfig.
getParameter<std::vector<double> >((
"Pt3050_"+lFullCutType).c_str());
65 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][0][i2] = pt010 [i2];
66 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][1][i2] = pt1020[i2];
67 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][2][i2] = pt2030[i2];
68 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][3][i2] = pt3050[i2];
71 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][0][i2] = pt010 [i2];
72 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][1][i2] = pt1020[i2];
73 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][2][i2] = pt2030[i2];
74 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][3][i2] = pt3050[i2];
77 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][0][i2] = pt010 [i2];
78 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][1][i2] = pt1020[i2];
79 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][2][i2] = pt2030[i2];
80 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][3][i2] = pt3050[i2];
92 Float_t impactParTkThreshod,
125 void assign(
const std::vector<float> & vec,
float &
a,
float &
b,
float &
c,
float &
d )
127 size_t sz = vec.size();
128 a = ( sz > 0 ? vec[0] : 0. );
129 b = ( sz > 1 ? vec[1] : 0. );
130 c = ( sz > 2 ? vec[2] : 0. );
131 d = ( sz > 3 ? vec[3] : 0. );
145 reader_jteta_0_2_ = std::unique_ptr<TMVA::Reader>(
new TMVA::Reader(
"!Color:Silent"));
148 reader_jteta_3_5_ = std::unique_ptr<TMVA::Reader>(
new TMVA::Reader(
"!Color:Silent"));
150 reader_ = std::unique_ptr<TMVA::Reader>(
new TMVA::Reader(
"!Color:Silent"));
230 if(jetPt >= 10 && jetPt < 20) ptId = 1;
231 if(jetPt >= 20 && jetPt < 30) ptId = 2;
232 if(jetPt >= 30 ) ptId = 3;
239 return std::pair<int,int>(ptId,etaId);
244 std::pair<int,int> jetIdKey =
getJetIdKey(jetPt,jetEta);
245 float betaStarModified = betaStarClassic/
log(nvtx-0.64);
263 std::pair<int,int> jetIdKey =
getJetIdKey(jetPt,jetEta);
290 static std::atomic<int> printWarning{10};
298 assert( patjet !=
nullptr || pfjet !=
nullptr );
299 if( patjet !=
nullptr && jec == 0. ) {
306 const reco::Candidate* lLead =
nullptr, *lSecond =
nullptr, *lLeadNeut =
nullptr, *lLeadEm =
nullptr, *lLeadCh =
nullptr, *lTrail =
nullptr;
307 std::vector<float>
frac, fracCh, fracEm, fracNeut;
308 float cones[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
309 size_t ncones =
sizeof(cones)/
sizeof(
float);
318 TMatrixDSym covMatrix(2); covMatrix = 0.;
320 float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
334 bool isPacked =
true;
335 if (lPack ==
nullptr){
339 float candPt = icand->
pt();
340 float candPtFrac = candPt/
jetPt;
344 float candPtDr = candPt * candDr;
345 size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
347 if(candDr < dRmin) dRmin = candDr;
350 if( lLead ==
nullptr || candPt > lLead->
pt() ) {
353 }
else if( (lSecond ==
nullptr || candPt > lSecond->pt()) && (candPt < lLead->
pt()) ) {
360 covMatrix(0,0) += candPt*candPt*candDeta*candDeta;
361 covMatrix(0,1) += candPt*candPt*candDeta*candDphi;
362 covMatrix(1,1) += candPt*candPt*candDphi*candDphi;
365 sumPt2 += candPt*candPt;
368 frac.push_back(candPtFrac);
370 if( icone < ncones ) { *coneFracs[icone] += candPt; }
374 if (lLeadNeut ==
nullptr || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
376 fracNeut.push_back(candPtFrac);
377 if( icone < ncones ) { *coneNeutFracs[icone] += candPt; }
382 if( icand->
pdgId() == 22 ) {
383 if(lLeadEm ==
nullptr || candPt > lLeadEm->pt()) { lLeadEm = icand; }
385 fracEm.push_back(candPtFrac);
386 if( icone < ncones ) { *coneEmFracs[icone] += candPt; }
391 if( icand->
charge() != 0 ) {
392 if (lLeadCh ==
nullptr || candPt > lLeadCh->pt()) {
396 if (lPF &&
std::abs(icand->
pdgId()) == 13 && pfTrk ==
nullptr){
401 edm::LogWarning(
"BadMuon")<<
"Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack ";
420 fracCh.push_back(candPtFrac);
421 if( icone < ncones ) { *coneChFracs[icone] += candPt; }
425 if( icand->
charge() != 0 ) {
433 bool inAnyOther =
false;
437 for(reco::VertexCollection::const_iterator vi=allvtx.begin(); vi!=allvtx.end(); ++vi ) {
439 if( iv.
isFake() || iv.
ndof() < 4 ) {
continue; }
444 if( ! isVtx0 && ! inAnyOther ) {
451 if( inVtx0 && ! inAnyOther ) {
453 }
else if( ! inVtx0 && inAnyOther ) {
459 }
else if( dZ < 0.2 ) {
469 bool inVtxOther =
false;
471 if (lPack->
fromPV() == 0) inVtxOther =
true;
483 if( lTrail ==
nullptr || candPt < lTrail->
pt() ) {
490 assert( !(lLead ==
nullptr) );
492 if ( lSecond ==
nullptr ) { lSecond = lTrail; }
493 if ( lLeadNeut ==
nullptr ) { lLeadNeut = lTrail; }
494 if ( lLeadEm ==
nullptr ) { lLeadEm = lTrail; }
495 if ( lLeadCh ==
nullptr ) { lLeadCh = lTrail; }
507 float sumW(0.0), sumW2(0.0);
508 float sum_deta(0.0),sum_dphi(0.0);
509 float ave_deta(0.0), ave_dphi(0.0);
512 if (!(
part.isAvailable() &&
part.isNonnull()) ){
517 float weight2 = weight *
weight;
519 float deta =
part->eta() - jet->
eta();
521 sum_deta += deta*weight2;
522 sum_dphi += dphi*weight2;
524 ave_deta = sum_deta/sumW2;
525 ave_dphi = sum_dphi/sumW2;
529 float ddetaR_sum(0.0), ddphiR_sum(0.0), pull_tmp(0.0);
532 if (!(
part.isAvailable() &&
part.isNonnull()) ){
536 float deta =
part->eta() - jet->
eta();
538 float ddeta, ddphi, ddR;
539 ddeta = deta - ave_deta ;
540 ddphi = dphi-ave_dphi;
541 ddR =
sqrt(ddeta*ddeta + ddphi*ddphi);
542 ddetaR_sum += ddR*ddeta*
weight;
543 ddphiR_sum += ddR*ddphi*
weight;
546 float ddetaR_ave = ddetaR_sum/sumW2;
547 float ddphiR_ave = ddphiR_sum/sumW2;
548 pull_tmp =
sqrt(ddetaR_ave*ddetaR_ave+ddphiR_ave*ddphiR_ave);
560 std::sort(frac.begin(),frac.end(),std::greater<float>());
561 std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
562 std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
563 std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
569 covMatrix(0,0) /= sumPt2;
570 covMatrix(0,1) /= sumPt2;
571 covMatrix(1,1) /= sumPt2;
572 covMatrix(1,0) = covMatrix(0,1);
576 TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
589 for(
size_t ic=0; ic<ncones; ++ic){
590 *coneFracs[ic] /=
jetPt;
591 *coneEmFracs[ic] /=
jetPt;
592 *coneNeutFracs[ic] /=
jetPt;
593 *coneChFracs[ic] /=
jetPt;
598 for(
unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
616 if( sumTkPt != 0. ) {
637 std::stringstream
out;
638 for(variables_list_t::const_iterator it=
variables_.begin();
640 out << std::setw(15) << it->first << std::setw(3) <<
"="
641 << std::setw(5) << *it->second.first
642 <<
" (" << std::setw(5) << it->second.second <<
")" << std::endl;
651 for(variables_list_t::iterator it=
variables_.begin();
653 *it->second.first = it->second.second;
658 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
659 internalId_.NAME ## _ = VAL; \
660 variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
T getParameter(std::string const &) const
void set(const PileupJetIdentifier &)
std::string tmvaWeights_jteta_2p5_3_
Float_t impactParTkThreshod_
std::vector< std::string > tmvaVariables_
bool isNonnull() const
Checks for non-null.
#define INIT_VARIABLE(NAME, TMVANAME, VAL)
std::string tmvaWeights_jteta_2_2p5_
float chargedEmEnergy() const
chargedEmEnergy
virtual double energy() const final
energy
CandidatePtr daughterPtr(size_type i) const
reference to daughter at given position
virtual const Track * bestTrack() const
virtual double pt() const =0
transverse momentum
T const * get() const
Returns C++ pointer to the item.
Base class for all types of Jets.
std::unique_ptr< TMVA::Reader > reader_jteta_2p5_3_
std::pair< int, int > getJetIdKey(float jetPt, float jetEta)
std::map< std::string, std::string > tmvaNames_
std::unique_ptr< TMVA::Reader > reader_jteta_0_2_
std::vector< Vertex > VertexCollection
collection of Vertex objects
int computeCutIDflag(float betaStarClassic, float dR2Mean, float nvtx, float jetPt, float jetEta)
Float_t betaStarCut_[3][4][4]
int chargedMultiplicity() const
chargedMultiplicity
const Point & position() const
position
std::unique_ptr< TMVA::Reader > reader_jteta_3_5_
Jets made from PFObjects.
float neutralEmEnergy() const
neutralEmEnergy
std::vector< std::string > tmvaVariables_jteta_0_3_
reco::TrackRef trackRef() const
virtual const Track * bestTrack() const
best track pointer
virtual size_t numberOfDaughters() const
number of daughters
virtual CandidatePtr sourceCandidatePtr(size_type i) const
PileupJetIdentifier internalId_
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
const PVAssoc fromPV(size_t ipv=0) const
std::string dumpVariables() const
std::vector< std::string > tmvaVariables_jteta_3_5_
Abs< T >::type abs(const T &t)
virtual const reco::Track & pseudoTrack() const
Return reference to a pseudo track made with candidate kinematics, parameterized error for eta...
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
void assign(const std::vector< float > &vec, float &a, float &b, float &c, float &d)
float trackWeight(const TREF &r) const
returns the weight with which a Track has contributed to the vertex-fit.
virtual int charge() const =0
electric charge
virtual size_type numberOfSourceCandidatePtrs() const
T const * get() const
Returns C++ pointer to the item.
PileupJetIdAlgo(int version=PHILv0, const std::string &tmvaWeight="", const std::string &tmvaMethod="", Float_t impactParTkThreshod_=1., const std::vector< std::string > &tmvaVariables=std::vector< std::string >(), bool runMvas=true)
virtual double mass() const final
mass
int computeIDflag(float mva, float jetPt, float jetEta)
int neutralMultiplicity() const
neutralMultiplicity
reco::MuonRef muonRef() const
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
virtual int pdgId() const =0
PDG identifier.
void setPtEtaPhi(const reco::Candidate &p, float &pt, float &eta, float &phi)
double deltaPhi(double phi1, double phi2)
std::string tmvaWeights_jteta_3_5_
PileupJetIdentifier computeIdVariables(const reco::Jet *jet, float jec, const reco::Vertex *, const reco::VertexCollection &, double rho)
unsigned int nCharged(const GenJet &jet)
PileupJetIdentifier computeMva()
Analysis-level calorimeter jet class.
std::vector< std::string > tmvaSpectators_
TMVA::IMethod * loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)
Geom::Phi< T > phi() const
Particle reconstructed by the particle flow algorithm.
std::unique_ptr< TMVA::Reader > reader_
reco::GsfTrackRef gsfTrackRef() const
variables_list_t variables_
float neutralHadronEnergy() const
neutralHadronEnergy
std::unique_ptr< TMVA::Reader > reader_jteta_2_2p5_
virtual double eta() const final
momentum pseudorapidity
std::string fullPath() const
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Jet correctedJet(const std::string &level, const std::string &flavor="none", const std::string &set="") const
float chargedHadronEnergy() const
chargedHadronEnergy
std::string tmvaWeights_jteta_0_2_
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity
virtual double pt() const final
transverse momentum