306 typedef std::map<const Track *, TransientTrack> TransientTrackMap;
327 if (groomedFatJetsHandle->size() > fatJetsHandle->size())
329 <<
"There are more groomed (" << groomedFatJetsHandle->size() <<
") than original fat jets (" 330 << fatJetsHandle->size() <<
"). Please check that the two jet collections belong to each other.";
335 unsigned int bsCovSrc[7] = {
338 double sigmaZ = 0.0, beamWidth = 0.0;
342 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = bsCovSrc[6] = 1;
343 sigmaZ = beamSpot->
sigmaZ();
349 bsCovSrc[0] = bsCovSrc[1] = 2;
350 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = 1;
351 sigmaZ = beamSpot->
sigmaZ();
355 bsCovSrc[0] = bsCovSrc[1] = bsCovSrc[2] = 2;
368 std::vector<std::vector<int> > clusteredSVs(trackIPTagInfos->size(), std::vector<int>());
371 std::vector<fastjet::PseudoJet> fjInputs;
375 std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
376 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
377 for (m = constituents.begin(); m != constituents.end(); ++
m) {
379 if (constit->
pt() == 0) {
380 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
383 fjInputs.push_back(fastjet::PseudoJet(constit->
px(), constit->
py(), constit->
pz(), constit->
energy()));
387 for (
typename std::vector<IPTI>::const_iterator it = trackIPTagInfos->begin(); it != trackIPTagInfos->end();
389 std::vector<edm::Ptr<reco::Candidate> > constituents = it->jet()->getJetConstituents();
390 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
391 for (m = constituents.begin(); m != constituents.end(); ++
m) {
393 if (constit->
pt() == 0) {
394 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
397 fjInputs.push_back(fastjet::PseudoJet(constit->
px(), constit->
py(), constit->
pz(), constit->
energy()));
403 const reco::Vertex &
pv = *(trackIPTagInfos->front().primaryVertex());
406 fastjet::PseudoJet
p(
407 dir.
x(), dir.
y(), dir.
z(), dir.
mag());
409 p = fastjet::PseudoJet(it->p4().px(), it->p4().py(), it->p4().pz(), it->p4().energy());
411 p.set_user_info(
new VertexInfo(it - extSecVertex->begin()));
412 fjInputs.push_back(p);
418 std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(
fjClusterSeq->inclusive_jets(
jetPtMin));
421 if (inclusiveJets.size() < fatJetsHandle->size())
423 <<
"There are fewer reclustered (" << inclusiveJets.size() <<
") than original fat jets (" 424 << fatJetsHandle->size()
425 <<
"). Please check that the jet algorithm and jet size match those used for the original jet collection.";
428 std::vector<int> reclusteredIndices;
429 matchReclusteredJets<edm::View<reco::Jet> >(fatJetsHandle, inclusiveJets, reclusteredIndices,
"fat");
432 std::vector<int> groomedIndices;
437 std::vector<std::vector<int> > subjetIndices;
439 matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
441 matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
444 for (
size_t i = 0;
i < fatJetsHandle->size(); ++
i) {
445 if (reclusteredIndices.at(
i) < 0)
448 if (fatJetsHandle->at(
i).pt() == 0)
451 <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
455 if (subjetIndices.at(
i).empty())
459 if ((
std::abs(inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - fatJetsHandle->at(
i).pt()) /
461 if (fatJetsHandle->at(
i).pt() < 10.)
463 <<
"The reclustered and original fat jet " <<
i <<
" have different Pt's (" 464 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << fatJetsHandle->at(
i).pt()
465 <<
" GeV, respectively).\n" 466 <<
"Please check that the jet algorithm and jet size match those used for the original fat jet " 467 "collection and also make sure the original fat jets are uncorrected. In addition, make sure you " 468 "are not using CaloJets which are presently not supported.\n" 469 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n" 470 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine " 471 "precision in which case make sure the original jet collection is produced and reclustering is " 472 "performed in the same job.";
475 <<
"The reclustered and original fat jet " <<
i <<
" have different Pt's (" 476 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << fatJetsHandle->at(
i).pt()
477 <<
" GeV, respectively).\n" 478 <<
"Please check that the jet algorithm and jet size match those used for the original fat jet " 479 "collection and also make sure the original fat jets are uncorrected. In addition, make sure you " 480 "are not using CaloJets which are presently not supported.\n" 481 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine " 482 "precision in which case make sure the original jet collection is produced and reclustering is " 483 "performed in the same job.";
487 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
489 std::vector<int> svIndices;
491 for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
493 if (!it->has_user_info())
496 svIndices.push_back(it->user_info<VertexInfo>().vertexIndex());
500 for (
size_t sv = 0;
sv < svIndices.size(); ++
sv) {
501 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
502 const VTX &extSV = (*extSecVertex)[svIndices.at(
sv)];
505 fastjet::PseudoJet
p(
506 dir.
x(), dir.
y(), dir.
z(), dir.
mag());
508 p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
510 std::vector<double> dR2toSubjets;
512 for (
size_t sj = 0; sj < subjetIndices.at(
i).size(); ++sj)
515 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->rapidity(),
516 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->phi()));
519 int closestSubjetIdx =
520 std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
522 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back(svIndices.at(
sv));
526 if (inclusiveJets.size() < trackIPTagInfos->size())
528 <<
"There are fewer reclustered (" << inclusiveJets.size() <<
") than original jets (" 529 << trackIPTagInfos->size()
530 <<
"). Please check that the jet algorithm and jet size match those used for the original jet collection.";
533 std::vector<int> reclusteredIndices;
534 matchReclusteredJets<std::vector<IPTI> >(
trackIPTagInfos, inclusiveJets, reclusteredIndices);
537 for (
size_t i = 0;
i < trackIPTagInfos->size(); ++
i) {
538 if (reclusteredIndices.at(
i) < 0)
541 if (trackIPTagInfos->at(
i).jet()->pt() == 0)
544 <<
"The original jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
549 if ((
std::abs(inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - trackIPTagInfos->at(
i).jet()->pt()) /
551 if (trackIPTagInfos->at(
i).jet()->pt() < 10.)
553 <<
"The reclustered and original jet " <<
i <<
" have different Pt's (" 554 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << trackIPTagInfos->at(
i).jet()->pt()
555 <<
" GeV, respectively).\n" 556 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection " 557 "and also make sure the original jets are uncorrected. In addition, make sure you are not using " 558 "CaloJets which are presently not supported.\n" 559 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n" 560 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine " 561 "precision in which case make sure the original jet collection is produced and reclustering is " 562 "performed in the same job.";
565 <<
"The reclustered and original jet " <<
i <<
" have different Pt's (" 566 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << trackIPTagInfos->at(
i).jet()->pt()
567 <<
" GeV, respectively).\n" 568 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection " 569 "and also make sure the original jets are uncorrected. In addition, make sure you are not using " 570 "CaloJets which are presently not supported.\n" 571 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine " 572 "precision in which case make sure the original jet collection is produced and reclustering is " 573 "performed in the same job.";
577 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
580 for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
582 if (!it->has_user_info())
585 clusteredSVs.at(
i).push_back(it->user_info<VertexInfo>().vertexIndex());
593 std::vector<int> groomedIndices;
598 std::vector<std::vector<int> > subjetIndices;
600 matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
602 matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
605 for (
size_t i = 0;
i < fatJetsHandle->size(); ++
i) {
606 if (fatJetsHandle->at(
i).pt() == 0)
609 <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
613 if (subjetIndices.at(
i).empty())
619 size_t sv = (it - extSecVertex->begin());
621 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
622 const VTX &extSV = (*extSecVertex)[
sv];
624 GlobalVector jetDir(fatJetsHandle->at(
i).px(), fatJetsHandle->at(
i).py(), fatJetsHandle->at(
i).pz());
630 fastjet::PseudoJet
p(
631 dir.
x(), dir.
y(), dir.
z(), dir.
mag());
633 p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
635 std::vector<double> dR2toSubjets;
637 for (
size_t sj = 0; sj < subjetIndices.at(
i).size(); ++sj)
640 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->rapidity(),
641 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->phi()));
644 int closestSubjetIdx =
645 std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
647 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back(sv);
653 std::unique_ptr<ConfigurableVertexReconstructor>
vertexReco;
654 std::unique_ptr<GhostTrackVertexFinder> vertexRecoGT;
665 TransientTrackMap primariesMap;
669 auto tagInfos = std::make_unique<Product>();
671 for (
typename std::vector<IPTI>::const_iterator iterJets = trackIPTagInfos->begin();
672 iterJets != trackIPTagInfos->end();
677 const Vertex &pv = *iterJets->primaryVertex();
679 std::set<TransientTrack> primaries;
682 TransientTrackMap::iterator
pos = primariesMap.lower_bound(iter->get());
684 if (pos != primariesMap.end() && pos->first == iter->get())
685 primaries.insert(pos->second);
688 primariesMap.insert(pos, std::make_pair(iter->get(),
track));
689 primaries.insert(track);
696 GlobalVector jetDir(jetRef->momentum().x(), jetRef->momentum().y(), jetRef->momentum().z());
702 const std::vector<reco::btag::TrackIPData> &ipData = iterJets->impactParameterData();
706 std::vector<TransientTrack> fitTracks;
707 std::vector<GhostTrackState> gtStates;
708 std::unique_ptr<GhostTrackPrediction> gtPred;
712 for (
unsigned int i = 0;
i < indices.size();
i++) {
718 trackData.back().first = indices[
i];
728 TransientTrackMap::const_iterator pos = primariesMap.find(
reco::btag::toTrack((trackRef)));
730 if (pos != primariesMap.end()) {
731 primaries.erase(pos->second);
732 fitTrack = pos->second;
734 fitTrack = trackBuilder->
build(trackRef);
735 fitTracks.push_back(fitTrack);
741 GlobalPoint pos = ipData[indices[
i]].closestToGhostTrack;
742 gtState.linearize(*gtPred,
true, gtPred->lambda(pos));
743 gtState.setWeight(ipData[indices[i]].ghostTrackWeight);
744 gtStates.push_back(gtState);
748 std::unique_ptr<GhostTrack> ghostTrack;
754 GlobalVector(iterJets->ghostTrack()->px(), iterJets->ghostTrack()->py(), iterJets->ghostTrack()->pz()),
758 iterJets->ghostTrack()->chi2(),
759 iterJets->ghostTrack()->ndof()));
763 std::vector<VTX> extAssoCollection;
764 std::vector<TransientVertex> fittedSVs;
765 std::vector<SecondaryVertex> SVs;
770 fittedSVs = vertexRecoGT->vertices(pv, *ghostTrack);
772 fittedSVs = vertexReco->vertices(fitTracks);
777 fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, *ghostTrack);
779 fittedSVs = vertexReco->vertices(fitTracks, *beamSpot);
786 for (
unsigned int i = 0; i < 7; i++) {
787 unsigned int covSrc = bsCovSrc[
i];
788 for (
unsigned int j = 0;
j < 7;
j++) {
790 if (!covSrc || bsCovSrc[
j] != covSrc)
792 else if (covSrc == 1)
794 else if (
j < 3 && i < 3)
809 fittedSVs = vertexRecoGT->vertices(pv, bs, *ghostTrack);
811 fittedSVs = vertexReco->vertices(fitTracks, bs);
815 std::vector<TransientTrack> primaries_(primaries.begin(), primaries.end());
817 fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, primaries_, *ghostTrack);
819 fittedSVs = vertexReco->vertices(primaries_, fitTracks, *beamSpot);
824 std::remove_copy_if(boost::make_transform_iterator(fittedSVs.begin(), svBuilder),
825 boost::make_transform_iterator(fittedSVs.end(), svBuilder),
826 std::back_inserter(SVs),
831 size_t jetIdx = (iterJets - trackIPTagInfos->begin());
833 for (
size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++) {
834 const VTX &extVertex = (*extSecVertex)[clusteredSVs.at(jetIdx).at(iExtSv)];
835 if (extVertex.p4().M() < 0.3)
837 extAssoCollection.push_back(extVertex);
840 for (
size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++) {
841 const VTX &extVertex = (*extSecVertex)[iExtSv];
844 extVertex.p4().M() < 0.3)
846 extAssoCollection.push_back(extVertex);
851 std::remove_copy_if(boost::make_transform_iterator(extAssoCollection.begin(), svBuilder),
852 boost::make_transform_iterator(extAssoCollection.end(), svBuilder),
853 std::back_inserter(SVs),
862 extAssoCollection.clear();
868 std::vector<typename TemplatedSecondaryVertexTagInfo<IPTI, VTX>::VertexData> svData;
870 svData.resize(vtxIndices.size());
871 for (
unsigned int idx = 0;
idx < vtxIndices.size();
idx++) {
874 svData[
idx].vertex =
sv;
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)
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
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
trackRef_iterator tracks_begin() const
first iterator over tracks
virtual double px() const =0
x coordinate of momentum vector
JetDefPtr fjJetDefinition
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector<TrackRef>
std::pair< unsigned int, TrackData > IndexedTrackData
Global3DVector GlobalVector
IPTI::input_container input_container
Public Types inherited from edm::stream::EDProducer<>
1.8.11
