326 typedef std::map<const Track *, TransientTrack> TransientTrackMap;
350 if( groomedFatJetsHandle->size() > fatJetsHandle->size() )
351 edm::LogError(
"TooManyGroomedJets") <<
"There are more groomed (" << groomedFatJetsHandle->size() <<
") than original fat jets (" << fatJetsHandle->size() <<
"). Please check that the two jet collections belong to each other.";
356 unsigned int bsCovSrc[7] = { 0, };
357 double sigmaZ = 0.0, beamWidth = 0.0;
361 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = bsCovSrc[6] = 1;
362 sigmaZ = beamSpot->
sigmaZ();
368 bsCovSrc[0] = bsCovSrc[1] = 2;
369 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = 1;
370 sigmaZ = beamSpot->
sigmaZ();
374 bsCovSrc[0] = bsCovSrc[1] = bsCovSrc[2] = 2;
387 std::vector<std::vector<int> > clusteredSVs(trackIPTagInfos->size(),std::vector<int>());
391 std::vector<fastjet::PseudoJet> fjInputs;
397 std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
398 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
399 for( m = constituents.begin(); m != constituents.end(); ++
m )
402 if(constit->
pt() == 0)
404 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
407 fjInputs.push_back(fastjet::PseudoJet(constit->
px(),constit->
py(),constit->
pz(),constit->
energy()));
413 for(
typename std::vector<IPTI>::const_iterator it = trackIPTagInfos->begin(); it != trackIPTagInfos->end(); ++it)
415 std::vector<edm::Ptr<reco::Candidate> > constituents = it->jet()->getJetConstituents();
416 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
417 for( m = constituents.begin(); m != constituents.end(); ++
m )
420 if(constit->
pt() == 0)
422 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
425 fjInputs.push_back(fastjet::PseudoJet(constit->
px(),constit->
py(),constit->
pz(),constit->
energy()));
432 const reco::Vertex &
pv = *(trackIPTagInfos->front().primaryVertex());
435 fastjet::PseudoJet
p(dir.
x(),dir.
y(),dir.
z(),dir.
mag());
437 p = fastjet::PseudoJet(it->p4().px(),it->p4().py(),it->p4().pz(),it->p4().energy());
439 p.set_user_info(
new VertexInfo( it - extSecVertex->begin() ));
440 fjInputs.push_back(p);
446 std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(
fjClusterSeq->inclusive_jets(
jetPtMin) );
450 if( inclusiveJets.size() < fatJetsHandle->size() )
451 edm::LogError(
"TooFewReclusteredJets") <<
"There are fewer reclustered (" << inclusiveJets.size() <<
") than original fat jets (" << fatJetsHandle->size() <<
"). Please check that the jet algorithm and jet size match those used for the original jet collection.";
454 std::vector<int> reclusteredIndices;
455 matchReclusteredJets<edm::View<reco::Jet> >(fatJetsHandle,inclusiveJets,reclusteredIndices,
"fat");
458 std::vector<int> groomedIndices;
463 std::vector<std::vector<int> > subjetIndices;
465 matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices);
467 matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices);
470 for(
size_t i=0;
i<fatJetsHandle->size(); ++
i)
472 if( reclusteredIndices.at(
i) < 0 )
continue;
474 if( fatJetsHandle->at(
i).pt() == 0 )
476 edm::LogWarning(
"NullTransverseMomentum") <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
480 if( subjetIndices.at(
i).empty() )
continue;
483 if( (
std::abs( inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - fatJetsHandle->at(
i).pt() ) / fatJetsHandle->at(
i).pt() ) >
relPtTolerance )
485 if( fatJetsHandle->at(
i).pt() < 10. )
486 edm::LogWarning(
"JetPtMismatchAtLowPt") <<
"The reclustered and original fat jet " <<
i <<
" have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << fatJetsHandle->at(
i).pt() <<
" GeV, respectively).\n" 487 <<
"Please check that the jet algorithm and jet size match those used for the original fat jet collection and also make sure the original fat jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n" 488 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n" 489 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
491 edm::LogError(
"JetPtMismatch") <<
"The reclustered and original fat jet " <<
i <<
" have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << fatJetsHandle->at(
i).pt() <<
" GeV, respectively).\n" 492 <<
"Please check that the jet algorithm and jet size match those used for the original fat jet collection and also make sure the original fat jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n" 493 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
497 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
499 std::vector<int> svIndices;
501 for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
503 if( !it->has_user_info() )
continue;
505 svIndices.push_back( it->user_info<VertexInfo>().vertexIndex() );
509 for(
size_t sv=0;
sv<svIndices.size(); ++
sv)
511 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
512 const VTX &extSV = (*extSecVertex)[ svIndices.at(
sv) ];
515 fastjet::PseudoJet
p(dir.
x(),dir.
y(),dir.
z(),dir.
mag());
517 p = fastjet::PseudoJet(extSV.p4().px(),extSV.p4().py(),extSV.p4().pz(),extSV.p4().energy());
519 std::vector<double> dR2toSubjets;
521 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
522 dR2toSubjets.push_back(
Geom::deltaR2( p.rapidity(), p.phi_std(), trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->rapidity(), trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->phi() ) );
525 int closestSubjetIdx =
std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
527 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back( svIndices.at(
sv) );
533 if( inclusiveJets.size() < trackIPTagInfos->size() )
534 edm::LogError(
"TooFewReclusteredJets") <<
"There are fewer reclustered (" << inclusiveJets.size() <<
") than original jets (" << trackIPTagInfos->size() <<
"). Please check that the jet algorithm and jet size match those used for the original jet collection.";
537 std::vector<int> reclusteredIndices;
538 matchReclusteredJets<std::vector<IPTI> >(
trackIPTagInfos,inclusiveJets,reclusteredIndices);
541 for(
size_t i=0;
i<trackIPTagInfos->size(); ++
i)
543 if( reclusteredIndices.at(
i) < 0 )
continue;
545 if( trackIPTagInfos->at(
i).jet()->pt() == 0 )
547 edm::LogWarning(
"NullTransverseMomentum") <<
"The original jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
552 if( (
std::abs( inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - trackIPTagInfos->at(
i).jet()->pt() ) / trackIPTagInfos->at(
i).jet()->pt() ) >
relPtTolerance )
554 if( trackIPTagInfos->at(
i).jet()->pt() < 10. )
555 edm::LogWarning(
"JetPtMismatchAtLowPt") <<
"The reclustered and original jet " <<
i <<
" have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << trackIPTagInfos->at(
i).jet()->pt() <<
" GeV, respectively).\n" 556 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n" 557 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n" 558 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
560 edm::LogError(
"JetPtMismatch") <<
"The reclustered and original jet " <<
i <<
" have different Pt's (" << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << trackIPTagInfos->at(
i).jet()->pt() <<
" GeV, respectively).\n" 561 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection and also make sure the original jets are uncorrected. In addition, make sure you are not using CaloJets which are presently not supported.\n" 562 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine precision in which case make sure the original jet collection is produced and reclustering is performed in the same job.";
566 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
569 for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
571 if( !it->has_user_info() )
continue;
573 clusteredSVs.at(
i).push_back( it->user_info<VertexInfo>().vertexIndex() );
582 std::vector<int> groomedIndices;
587 std::vector<std::vector<int> > subjetIndices;
589 matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices);
591 matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices);
594 for(
size_t i=0;
i<fatJetsHandle->size(); ++
i)
596 if( fatJetsHandle->at(
i).pt() == 0 )
598 edm::LogWarning(
"NullTransverseMomentum") <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
602 if( subjetIndices.at(
i).empty() )
continue;
608 size_t sv = ( it - extSecVertex->begin() );
610 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
611 const VTX &extSV = (*extSecVertex)[
sv];
614 fatJetsHandle->at(
i).py(),
615 fatJetsHandle->at(
i).pz());
621 fastjet::PseudoJet
p(dir.
x(),dir.
y(),dir.
z(),dir.
mag());
623 p = fastjet::PseudoJet(extSV.p4().px(),extSV.p4().py(),extSV.p4().pz(),extSV.p4().energy());
625 std::vector<double> dR2toSubjets;
627 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
628 dR2toSubjets.push_back(
Geom::deltaR2( p.rapidity(), p.phi_std(), trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->rapidity(), trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->phi() ) );
631 int closestSubjetIdx =
std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
633 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back(sv);
639 std::unique_ptr<ConfigurableVertexReconstructor>
vertexReco;
640 std::unique_ptr<GhostTrackVertexFinder> vertexRecoGT;
652 TransientTrackMap primariesMap;
656 auto tagInfos = std::make_unique<Product>();
658 for(
typename std::vector<IPTI>::const_iterator iterJets =
659 trackIPTagInfos->begin(); iterJets != trackIPTagInfos->end();
664 const Vertex &pv = *iterJets->primaryVertex();
666 std::set<TransientTrack> primaries;
670 TransientTrackMap::iterator
pos =
671 primariesMap.lower_bound(iter->get());
673 if (pos != primariesMap.end() &&
674 pos->first == iter->get())
675 primaries.insert(pos->second);
680 primariesMap.insert(pos,
681 std::make_pair(iter->get(),
track));
682 primaries.insert(track);
690 jetRef->momentum().y(),
691 jetRef->momentum().z());
693 std::vector<std::size_t> indices =
698 const std::vector<reco::btag::TrackIPData> &ipData =
699 iterJets->impactParameterData();
703 std::vector<TransientTrack> fitTracks;
704 std::vector<GhostTrackState> gtStates;
705 std::unique_ptr<GhostTrackPrediction> gtPred;
708 *iterJets->ghostTrack()));
710 for(
unsigned int i = 0;
i < indices.size();
i++) {
716 trackData.back().first = indices[
i];
723 trackData.back().second.svStatus =
728 TransientTrackMap::const_iterator pos =
731 if (pos != primariesMap.end()) {
732 primaries.erase(pos->second);
733 fitTrack = pos->second;
735 fitTrack = trackBuilder->
build(trackRef);
736 fitTracks.push_back(fitTrack);
738 trackData.back().second.svStatus =
744 ipData[indices[
i]].closestToGhostTrack;
745 gtState.linearize(*gtPred,
true,
746 gtPred->lambda(pos));
747 gtState.setWeight(ipData[indices[i]].ghostTrackWeight);
748 gtStates.push_back(gtState);
752 std::unique_ptr<GhostTrack> ghostTrack;
759 iterJets->ghostTrack()->px(),
760 iterJets->ghostTrack()->py(),
761 iterJets->ghostTrack()->pz()),
764 iterJets->ghostTrack()->chi2(),
765 iterJets->ghostTrack()->ndof()));
770 std::vector<VTX> extAssoCollection;
771 std::vector<TransientVertex> fittedSVs;
772 std::vector<SecondaryVertex> SVs;
777 fittedSVs = vertexRecoGT->vertices(
780 fittedSVs = vertexReco->vertices(fitTracks);
785 fittedSVs = vertexRecoGT->vertices(
786 pv, *beamSpot, *ghostTrack);
788 fittedSVs = vertexReco->vertices(fitTracks,
796 for(
unsigned int i = 0; i < 7; i++) {
797 unsigned int covSrc = bsCovSrc[
i];
798 for(
unsigned int j = 0; j < 7; j++) {
800 if (!covSrc || bsCovSrc[j] != covSrc)
802 else if (covSrc == 1)
817 fittedSVs = vertexRecoGT->vertices(
818 pv, bs, *ghostTrack);
820 fittedSVs = vertexReco->vertices(fitTracks, bs);
824 std::vector<TransientTrack> primaries_(
825 primaries.begin(), primaries.end());
827 fittedSVs = vertexRecoGT->vertices(
828 pv, *beamSpot, primaries_,
831 fittedSVs = vertexReco->vertices(
832 primaries_, fitTracks,
838 std::remove_copy_if(boost::make_transform_iterator(
839 fittedSVs.begin(), svBuilder),
840 boost::make_transform_iterator(
841 fittedSVs.end(), svBuilder),
842 std::back_inserter(SVs),
847 size_t jetIdx = ( iterJets - trackIPTagInfos->begin() );
849 for(
size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++){
850 const VTX & extVertex = (*extSecVertex)[ clusteredSVs.at(jetIdx).at(iExtSv) ];
851 if( extVertex.p4().M() < 0.3 )
853 extAssoCollection.push_back( extVertex );
857 for(
size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++){
858 const VTX & extVertex = (*extSecVertex)[iExtSv];
861 extAssoCollection.push_back( extVertex );
867 std::remove_copy_if(boost::make_transform_iterator( extAssoCollection.begin(), svBuilder),
868 boost::make_transform_iterator(extAssoCollection.end(), svBuilder),
869 std::back_inserter(SVs),
878 extAssoCollection.clear();
884 std::vector<typename TemplatedSecondaryVertexTagInfo<IPTI,VTX>::VertexData> svData;
886 svData.resize(vtxIndices.size());
887 for(
unsigned int idx = 0;
idx < vtxIndices.size();
idx++) {
890 svData[
idx].vertex =
sv;
903 trackData, svData, SVs.size(),
905 iterJets - trackIPTagInfos->begin())));
math::Error< dimension >::type CovarianceMatrix
reco::Vertex::Point convertPos(const GlobalPoint &p)
T getParameter(std::string const &) const
edm::EDGetTokenT< edm::View< reco::Jet > > token_fatJets
const math::XYZPoint & position(const reco::Vertex &sv)
virtual double pz() const =0
z coordinate of momentum vector
TrackSelector trackSelector
trackRef_iterator tracks_end() const
last iterator over tracks
ClusterSequencePtr fjClusterSeq
edm::EDGetTokenT< reco::BeamSpot > token_BeamSpot
reco::Vertex::Error convertError(const GlobalError &ge)
Measurement1D dist3d() const
reco::TransientTrack build(const reco::Track *p) const
Measurement1D dist2d() const
double covariance(int i, int j) const
(i, j)-th element of error matrix, i, j = 0, ... 2
VertexFilter vertexFilter
reco::btag::SortCriteria sortCriterium
const reco::Track * toTrack(const reco::TrackBaseRef &t)
const Point & position() const
position
static GhostTrackVertexFinder::FitType getGhostTrackFitType(const std::string &name)
ConstraintType constraint
double dydz() const
dydz slope
edm::EDGetTokenT< std::vector< IPTI > > token_trackIPTagInfo
IPTI::input_container::value_type input_item
virtual double energy() const =0
energy
virtual double py() const =0
y coordinate of momentum vector
reco::btag::IndexedTrackData IndexedTrackData
GlobalVector flightDirection(const reco::Vertex &pv, const reco::Vertex &sv)
boost::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
Abs< T >::type abs(const T &t)
Measurement1D dist1d() const
double BeamWidthX() const
beam width X
VertexSorting< SecondaryVertex > vertexSorting
void matchGroomedJets(const edm::Handle< edm::View< reco::Jet > > &jets, const edm::Handle< edm::View< reco::Jet > > &matchedJets, std::vector< int > &matchedIndices)
void matchSubjets(const std::vector< int > &groomedIndices, const edm::Handle< edm::View< reco::Jet > > &groomedJets, const edm::Handle< std::vector< IPTI > > &subjets, std::vector< std::vector< int > > &matchedIndices)
double dxdz() const
dxdz slope
void markUsedTracks(TrackDataVector &trackData, const input_container &trackRefs, const SecondaryVertex &sv, size_t idx)
edm::EDGetTokenT< edm::View< VTX > > token_extSVCollection
Vector3DBase unit() const
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
edm::ParameterSet vtxRecoPSet
edm::EDGetTokenT< edm::View< reco::Jet > > token_groomedFatJets
virtual double pt() const =0
transverse momentum
double sigmaZ() const
sigma z
Error error() const
return SMatrix
double covariance(int i, int j) const
(i,j)-th element of error matrix
std::vector< reco::btag::IndexedTrackData > TrackDataVector
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector<TrackRef>
trackRef_iterator tracks_begin() const
first iterator over tracks
virtual double px() const =0
x coordinate of momentum vector
JetDefPtr fjJetDefinition
std::pair< unsigned int, TrackData > IndexedTrackData
Global3DVector GlobalVector
IPTI::input_container input_container
Public Types inherited from edm::stream::EDProducer<>
1.8.11
