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,
207 void assign(
const std::vector<float> & vec,
float &
a,
float &
b,
float &
c,
float &
d )
209 size_t sz = vec.size();
210 a = ( sz > 0 ? vec[0] : 0. );
211 b = ( sz > 1 ? vec[1] : 0. );
212 c = ( sz > 2 ? vec[2] : 0. );
213 d = ( sz > 3 ? vec[3] : 0. );
227 reader_jteta_0_2_ = std::unique_ptr<TMVA::Reader>(
new TMVA::Reader(
"!Color:Silent"));
230 reader_jteta_3_5_ = std::unique_ptr<TMVA::Reader>(
new TMVA::Reader(
"!Color:Silent"));
232 reader_ = std::unique_ptr<TMVA::Reader>(
new TMVA::Reader(
"!Color:Silent"));
312 if(jetPt >= 10 && jetPt < 20) ptId = 1;
313 if(jetPt >= 20 && jetPt < 30) ptId = 2;
314 if(jetPt >= 30 ) ptId = 3;
321 return std::pair<int,int>(ptId,etaId);
326 std::pair<int,int> jetIdKey =
getJetIdKey(jetPt,jetEta);
327 float betaStarModified = betaStarClassic/
log(nvtx-0.64);
345 std::pair<int,int> jetIdKey =
getJetIdKey(jetPt,jetEta);
372 static std::atomic<int> printWarning{10};
380 assert( patjet !=
nullptr || pfjet !=
nullptr );
381 if( patjet !=
nullptr && jec == 0. ) {
388 const reco::Candidate* lLead =
nullptr, *lSecond =
nullptr, *lLeadNeut =
nullptr, *lLeadEm =
nullptr, *lLeadCh =
nullptr, *lTrail =
nullptr;
389 std::vector<float>
frac, fracCh, fracEm, fracNeut;
390 float cones[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
391 size_t ncones =
sizeof(cones)/
sizeof(
float);
400 TMatrixDSym covMatrix(2); covMatrix = 0.;
402 float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
416 bool isPacked =
true;
417 if (lPack ==
nullptr){
421 float candPt = icand->
pt();
422 float candPtFrac = candPt/
jetPt;
426 float candPtDr = candPt * candDr;
427 size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
429 if(candDr < dRmin) dRmin = candDr;
432 if( lLead ==
nullptr || candPt > lLead->
pt() ) {
435 }
else if( (lSecond ==
nullptr || candPt > lSecond->pt()) && (candPt < lLead->
pt()) ) {
442 covMatrix(0,0) += candPt*candPt*candDeta*candDeta;
443 covMatrix(0,1) += candPt*candPt*candDeta*candDphi;
444 covMatrix(1,1) += candPt*candPt*candDphi*candDphi;
447 sumPt2 += candPt*candPt;
450 frac.push_back(candPtFrac);
452 if( icone < ncones ) { *coneFracs[icone] += candPt; }
456 if (lLeadNeut ==
nullptr || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
458 fracNeut.push_back(candPtFrac);
459 if( icone < ncones ) { *coneNeutFracs[icone] += candPt; }
464 if( icand->
pdgId() == 22 ) {
465 if(lLeadEm ==
nullptr || candPt > lLeadEm->pt()) { lLeadEm = icand; }
467 fracEm.push_back(candPtFrac);
468 if( icone < ncones ) { *coneEmFracs[icone] += candPt; }
473 if( icand->
charge() != 0 ) {
474 if (lLeadCh ==
nullptr || candPt > lLeadCh->pt()) {
478 if (lPF &&
std::abs(icand->
pdgId()) == 13 && pfTrk ==
nullptr){
483 edm::LogWarning(
"BadMuon")<<
"Found a PFCandidate muon without a trackRef: falling back to Muon::bestTrack ";
502 fracCh.push_back(candPtFrac);
503 if( icone < ncones ) { *coneChFracs[icone] += candPt; }
507 if( icand->
charge() != 0 ) {
515 bool inAnyOther =
false;
519 for(reco::VertexCollection::const_iterator vi=allvtx.begin(); vi!=allvtx.end(); ++vi ) {
521 if( iv.
isFake() || iv.
ndof() < 4 ) {
continue; }
526 if( ! isVtx0 && ! inAnyOther ) {
533 if( inVtx0 && ! inAnyOther ) {
535 }
else if( ! inVtx0 && inAnyOther ) {
541 }
else if( dZ < 0.2 ) {
551 bool inVtxOther =
false;
553 if (lPack->
fromPV() == 0) inVtxOther =
true;
565 if( lTrail ==
nullptr || candPt < lTrail->
pt() ) {
572 assert( !(lLead ==
nullptr) );
574 if ( lSecond ==
nullptr ) { lSecond = lTrail; }
575 if ( lLeadNeut ==
nullptr ) { lLeadNeut = lTrail; }
576 if ( lLeadEm ==
nullptr ) { lLeadEm = lTrail; }
577 if ( lLeadCh ==
nullptr ) { lLeadCh = lTrail; }
593 std::sort(frac.begin(),frac.end(),std::greater<float>());
594 std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
595 std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
596 std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
602 covMatrix(0,0) /= sumPt2;
603 covMatrix(0,1) /= sumPt2;
604 covMatrix(1,1) /= sumPt2;
605 covMatrix(1,0) = covMatrix(0,1);
609 TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
622 for(
size_t ic=0; ic<ncones; ++ic){
623 *coneFracs[ic] /=
jetPt;
624 *coneEmFracs[ic] /=
jetPt;
625 *coneNeutFracs[ic] /=
jetPt;
626 *coneChFracs[ic] /=
jetPt;
631 for(
unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
649 if( sumTkPt != 0. ) {
670 std::stringstream
out;
671 for(variables_list_t::const_iterator it=
variables_.begin();
673 out << std::setw(15) << it->first << std::setw(3) <<
"="
674 << std::setw(5) << *it->second.first
675 <<
" (" << std::setw(5) << it->second.second <<
")" << std::endl;
684 for(variables_list_t::iterator it=
variables_.begin();
686 *it->second.first = it->second.second;
691 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
692 internalId_.NAME ## _ = VAL; \
693 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 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
double deltaR(const T1 &t1, const T2 &t2)
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.
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
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 double mass() const
mass
virtual double energy() const
energy
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)
float trackWeight(const TrackBaseRef &r) const
returns the weight with which a Track has contributed to the vertex-fit.
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...
void assign(const std::vector< float > &vec, float &a, float &b, float &c, float &d)
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)
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)
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_
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