10 #include "TMatrixDSym.h"
11 #include "TMatrixDSymEigen.h"
28 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++)
mvacut_[i0][0][i2] = pt010 [i2];
38 for(
int i2 = 0; i2 <
NPts; i2++)
mvacut_[i0][1][i2] = pt1020[i2];
39 for(
int i2 = 0; i2 <
NPts; i2++)
mvacut_[i0][2][i2] = pt2030[i2];
40 for(
int i2 = 0; i2 <
NPts; i2++)
mvacut_[i0][3][i2] = pt3050[i2];
54 Float_t impactParTkThreshod,
122 void Assign(
const std::vector<float> & vec,
float &
a,
float &
b,
float &
c,
float &
d )
124 size_t sz = vec.size();
125 a = ( sz > 0 ? vec[0] : 0. );
126 b = ( sz > 1 ? vec[1] : 0. );
127 c = ( sz > 2 ? vec[2] : 0. );
128 d = ( sz > 3 ? vec[3] : 0. );
142 reader_ =
new TMVA::Reader(
"!Color:Silent");
177 if(jetPt > 10 && jetPt < 20) ptId = 1;
178 if(jetPt >= 20 && jetPt < 30) ptId = 2;
179 if(jetPt >= 30 ) ptId = 3;
182 if(fabs(jetEta) > 2.5 && fabs(jetEta) < 2.75) etaId = 1;
183 if(fabs(jetEta) >= 2.75 && fabs(jetEta) < 3.0 ) etaId = 2;
184 if(fabs(jetEta) >= 3.0 && fabs(jetEta) < 5.0 ) etaId = 3;
185 return std::pair<int,int>(ptId,etaId);
191 std::pair<int,int> jetIdKey =
getJetIdKey(jetPt,jetEta);
215 typedef std::vector <reco::PFCandidatePtr> constituents_type;
216 typedef std::vector <reco::PFCandidatePtr>::iterator constituents_iterator;
229 std::vector<float>
frac, fracCh, fracEm, fracNeut;
231 std::array<float, ncones> cones{ {0.1,0.2,0.3,0.4}};
233 TMatrixDSym covMatrix(2); covMatrix = 0.;
237 float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0;
float sum_deta =0 ;
float sum_dphi =0 ;
float sum_deta2 =0 ;
float sum_detadphi =0 ;
float sum_dphi2=0;
241 for(constituents_iterator it=constituents.begin(); it!=constituents.end(); ++it) {
243 float candPt = icand->pt();
244 float candPtFrac = candPt/
jetPt;
246 float candDeta = fabs( (*it)->eta() - jet->
eta() );
248 float candPtDr = candPt * candDr;
249 size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
250 float weight2 = candPt * candPt;
253 if( lLead.
isNull() || candPt > lLead->pt() ) {
256 }
else if( (lSecond.
isNull() || candPt > lSecond->pt()) && (candPt < lLead->
pt()) ) {
265 sumPt2 += candPt*candPt;
266 sum_deta += candDeta*weight2;
267 sum_dphi += candDphi*weight2;
268 sum_deta2 += candDeta*candDeta*weight2;
269 sum_detadphi += candDeta*candDphi*weight2;
270 sum_dphi2 += candDphi*candDphi*weight2;
276 frac.push_back(candPtFrac);
277 if( icone < ncones ) { *coneFracs[icone] += candPt; }
281 if (lLeadNeut.
isNull() || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
283 fracNeut.push_back(candPtFrac);
288 if(lLeadEm.
isNull() || candPt > lLeadEm->pt()) { lLeadEm = icand; }
290 fracEm.push_back(candPtFrac);
295 if (lLeadCh.
isNull() || candPt > lLeadCh->pt()) { lLeadCh = icand; }
297 fracCh.push_back(candPtFrac);
302 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 ) {
311 if( iv.
isFake() || iv.
ndof() < 4 ) {
continue; }
315 if( ! isVtx0 && ! inAnyOther ) {
321 if( inVtx0 && ! inAnyOther ) {
323 }
else if( ! inVtx0 && inAnyOther ) {
329 }
else if( dZ < 0.2 ) {
334 if( lTrail.
isNull() || candPt < lTrail->pt() ) {
340 if ( lSecond.
isNull() ) { lSecond = lTrail; }
341 if ( lLeadNeut.
isNull() ) { lLeadNeut = lTrail; }
342 if ( lLeadEm.
isNull() ) { lLeadEm = lTrail; }
343 if ( lLeadCh.
isNull() ) { lLeadCh = lTrail; }
344 impactTrack = lLeadCh->trackRef();
363 std::sort(frac.begin(),frac.end(),std::greater<float>());
364 std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
365 std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
366 std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
386 for(
size_t ic=0; ic<ncones; ++ic){
387 *coneFracs[ic] /=
jetPt;
392 for(
unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
411 float ave_deta = sum_deta/sumPt2;
412 float ave_dphi = sum_dphi/sumPt2;
413 float ave_deta2 = sum_deta2/sumPt2;
414 float ave_dphi2 = sum_dphi2/sumPt2;
415 float a = ave_deta2-ave_deta*ave_deta;
416 float b = ave_dphi2-ave_dphi*ave_dphi;
417 float c = -(sum_detadphi/sumPt2-ave_deta*ave_dphi);
418 float axis1=0;
float axis2=0;
419 if((((a-b)*(a-b)+4*c*c))>0) {
422 axis1 =
sqrt(0.5*(a+b+delta));
425 axis2 =
sqrt(0.5*(a+b-delta));
436 float ddetaR_sum(0.0), ddphiR_sum(0.0);
439 float weight =part->pt()*part->pt() ;
440 float deta = part->eta() - jet->
eta();
442 float ddeta, ddphi, ddR;
443 ddeta = deta - ave_deta ;
445 ddR =
sqrt(ddeta*ddeta + ddphi*ddphi);
446 ddetaR_sum += ddR*ddeta*
weight;
447 ddphiR_sum += ddR*ddphi*
weight;
449 float ddetaR_ave = ddetaR_sum/sumPt2;
450 float ddphiR_ave = ddphiR_sum/sumPt2;
451 internalId_.dRMean_ =
sqrt(ddetaR_ave*ddetaR_ave+ddphiR_ave*ddphiR_ave);
467 std::stringstream
out;
468 for(variables_list_t::const_iterator it=
variables_.begin();
470 out << std::setw(15) << it->first << std::setw(3) <<
"="
471 << std::setw(5) << *it->second.first
472 <<
" (" << std::setw(5) << it->second.second <<
")" << std::endl;
480 for(variables_list_t::iterator it=
variables_.begin();
482 *it->second.first = it->second.second;
486 #define INIT_VARIABLE(NAME,TMVANAME,VAL) \
487 internalId_.NAME ## _ = VAL; \
488 variables_[ # NAME ] = std::make_pair(& internalId_.NAME ## _, VAL);
T getParameter(std::string const &) const
bool isNonnull() const
Checks for non-null.
std::vector< std::string > tmvaVariables_
float chargedEmEnergy() const
chargedEmEnergy
trackRef_iterator tracks_end() const
last iterator over tracks
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
double deltaR(const T1 &t1, const T2 &t2)
variables_list_t variables_
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.
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
float neutralEmEnergy() const
neutralEmEnergy
virtual reco::PFCandidatePtr getPFConstituent(unsigned fIndex) const
get specific constituent
std::map< std::string, std::string > tmvaNames_
virtual double mass() const
mass
virtual double energy() const
energy
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)
void SetPtEtaPhi(const reco::Candidate &p, float &pt, float &eta, float &phi)
int neutralMultiplicity() const
neutralMultiplicity
bool isNonnull() const
Checks for non-null.
double deltaPhi(double phi1, double phi2)
unsigned int nCharged(const GenJet &jet)
PileupJetIdentifier internalId_
TMVA::IMethod * loadTMVAWeights(TMVA::Reader *reader, const std::string &method, const std::string &weightFile, bool verbose=false)
Geom::Phi< T > phi() const
void set(const PileupJetIdentifier &)
std::pair< int, int > getJetIdKey(float jetPt, float jetEta)
PileupJetIdentifier computeMva()
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
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 >())
std::string fullPath() const
trackRef_iterator tracks_begin() const
first iterator over tracks
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)
virtual double phi() const =0
momentum azimuthal angle
virtual double eta() const =0
momentum pseudorapidity