11 #include "TMatrixDSym.h"
12 #include "TMatrixDSymEigen.h"
35 for(
int i0 = 0; i0 < 3; i0++) {
41 for(
int i1 = 0; i1 < nCut; i1++) {
43 if(
cutBased_ && i1 == 0) lFullCutType =
"BetaStar"+ lCutType;
44 if(
cutBased_ && i1 == 1) lFullCutType =
"RMS" + lCutType;
45 std::vector<double> pt010 = jetConfig.
getParameter<std::vector<double> >((
"Pt010_" +lFullCutType).c_str());
46 std::vector<double> pt1020 = jetConfig.
getParameter<std::vector<double> >((
"Pt1020_"+lFullCutType).c_str());
47 std::vector<double> pt2030 = jetConfig.
getParameter<std::vector<double> >((
"Pt2030_"+lFullCutType).c_str());
48 std::vector<double> pt3050 = jetConfig.
getParameter<std::vector<double> >((
"Pt3050_"+lFullCutType).c_str());
50 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][0][i2] = pt010 [i2];
51 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][1][i2] = pt1020[i2];
52 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][2][i2] = pt2030[i2];
53 for(
int i2 = 0; i2 < 4; i2++)
mvacut_[i0][3][i2] = pt3050[i2];
56 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][0][i2] = pt010 [i2];
57 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][1][i2] = pt1020[i2];
58 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][2][i2] = pt2030[i2];
59 for(
int i2 = 0; i2 < 4; i2++)
betaStarCut_[i0][3][i2] = pt3050[i2];
62 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][0][i2] = pt010 [i2];
63 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][1][i2] = pt1020[i2];
64 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][2][i2] = pt2030[i2];
65 for(
int i2 = 0; i2 < 4; i2++)
rmsCut_[i0][3][i2] = pt3050[i2];
77 Float_t impactParTkThreshod,
175 void assign(
const std::vector<float> & vec,
float &
a,
float &
b,
float &
c,
float &
d )
177 size_t sz = vec.size();
178 a = ( sz > 0 ? vec[0] : 0. );
179 b = ( sz > 1 ? vec[1] : 0. );
180 c = ( sz > 2 ? vec[2] : 0. );
181 d = ( sz > 3 ? vec[3] : 0. );
194 reader_ =
new TMVA::Reader(
"!Color:Silent");
234 if(jetPt >= 10 && jetPt < 20) ptId = 1;
235 if(jetPt >= 20 && jetPt < 30) ptId = 2;
236 if(jetPt >= 30 ) ptId = 3;
243 return std::pair<int,int>(ptId,etaId);
248 std::pair<int,int> jetIdKey =
getJetIdKey(jetPt,jetEta);
249 float betaStarModified = betaStarClassic/
log(nvtx-0.64);
267 std::pair<int,int> jetIdKey =
getJetIdKey(jetPt,jetEta);
295 static std::atomic<int> printWarning{10};
303 assert( patjet !=
nullptr || pfjet !=
nullptr );
304 if( patjet !=
nullptr && jec == 0. ) {
311 const reco::Candidate* lLead =
nullptr, *lSecond =
nullptr, *lLeadNeut =
nullptr, *lLeadEm =
nullptr, *lLeadCh =
nullptr, *lTrail =
nullptr;
312 std::vector<float>
frac, fracCh, fracEm, fracNeut;
313 float cones[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
314 size_t ncones =
sizeof(cones)/
sizeof(
float);
323 TMatrixDSym covMatrix(2); covMatrix = 0.;
325 float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
336 bool isPacked =
true;
337 if (lPack ==
nullptr){
341 float candPt = icand->
pt();
342 float candPtFrac = candPt/
jetPt;
346 float candPtDr = candPt * candDr;
347 size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
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; }
511 std::sort(frac.begin(),frac.end(),std::greater<float>());
512 std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
513 std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
514 std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
520 covMatrix(0,0) /= sumPt2;
521 covMatrix(0,1) /= sumPt2;
522 covMatrix(1,1) /= sumPt2;
523 covMatrix(1,0) = covMatrix(0,1);
527 TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
540 for(
size_t ic=0; ic<ncones; ++ic){
541 *coneFracs[ic] /=
jetPt;
542 *coneEmFracs[ic] /=
jetPt;
543 *coneNeutFracs[ic] /=
jetPt;
544 *coneChFracs[ic] /=
jetPt;
549 for(
unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
563 if( sumTkPt != 0. ) {
584 std::stringstream
out;
585 for(variables_list_t::const_iterator it=
variables_.begin();
587 out << std::setw(15) << it->first << std::setw(3) <<
"="
588 << std::setw(5) << *it->second.first
589 <<
" (" << std::setw(5) << it->second.second <<
")" << std::endl;
598 for(variables_list_t::iterator it=
variables_.begin();
600 *it->second.first = it->second.second;
605 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
606 internalId_.NAME ## _ = VAL; \
607 variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
T getParameter(std::string const &) const
void set(const PileupJetIdentifier &)
Float_t impactParTkThreshod_
std::vector< std::string > tmvaVariables_
bool isNonnull() const
Checks for non-null.
#define INIT_VARIABLE(NAME, TMVANAME, VAL)
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::pair< int, int > getJetIdKey(float jetPt, float jetEta)
std::map< std::string, std::string > tmvaNames_
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
Jets made from PFObjects.
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
float neutralEmEnergy() const
neutralEmEnergy
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
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.
PileupJetIdentifier computeIdVariables(const reco::Jet *jet, float jec, const reco::Vertex *, const reco::VertexCollection &, bool calculateMva=false)
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)
unsigned int nCharged(const GenJet &jet)
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 >())
PileupJetIdentifier computeMva()
Analysis-level calorimeter jet class.
std::vector< std::string > tmvaSpectators_
Particle reconstructed by the particle flow algorithm.
reco::GsfTrackRef gsfTrackRef() const
variables_list_t variables_
float neutralHadronEnergy() const
neutralHadronEnergy
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
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity
void loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)