320 typedef std::map<const Track *, TransientTrack> TransientTrackMap;
339 if (groomedFatJetsHandle->size() > fatJetsHandle->size())
341 <<
"There are more groomed (" << groomedFatJetsHandle->size() <<
") than original fat jets ("
342 << fatJetsHandle->size() <<
"). Please check that the two jet collections belong to each other.";
350 unsigned int bsCovSrc[7] = {
353 double sigmaZ = 0.0, beamWidth = 0.0;
357 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = bsCovSrc[6] = 1;
358 sigmaZ = beamSpot->sigmaZ();
359 beamWidth = beamSpot->BeamWidthX();
364 bsCovSrc[0] = bsCovSrc[1] = 2;
365 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = 1;
366 sigmaZ = beamSpot->sigmaZ();
370 bsCovSrc[0] = bsCovSrc[1] = bsCovSrc[2] = 2;
383 std::vector<std::vector<int> > clusteredSVs(trackIPTagInfos->size(), std::vector<int>());
386 std::vector<fastjet::PseudoJet> fjInputs;
390 std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
391 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
392 for (m = constituents.begin(); m != constituents.end(); ++
m) {
394 if (constit->pt() == 0) {
395 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
398 if (it->isWeighted()) {
401 <<
"TemplatedSecondaryVertexProducer: No weights (e.g. PUPPI) given for weighted jet collection"
403 float w = (*weightsHandle)[constit];
405 fastjet::PseudoJet(constit->px() *
w, constit->py() *
w, constit->pz() *
w, constit->energy() *
w));
407 fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
412 for (
typename std::vector<IPTI>::const_iterator it = trackIPTagInfos->begin(); it != trackIPTagInfos->end();
414 std::vector<edm::Ptr<reco::Candidate> > constituents = it->jet()->getJetConstituents();
415 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
416 for (m = constituents.begin(); m != constituents.end(); ++
m) {
418 if (constit->pt() == 0) {
419 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
422 if (it->jet()->isWeighted()) {
425 <<
"TemplatedSecondaryVertexProducer: No weights (e.g. PUPPI) given for weighted jet collection"
427 float w = (*weightsHandle)[constit];
429 fastjet::PseudoJet(constit->px() *
w, constit->py() *
w, constit->pz() *
w, constit->energy() *
w));
431 fjInputs.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
438 const reco::Vertex &
pv = *(trackIPTagInfos->front().primaryVertex());
441 fastjet::PseudoJet
p(
442 dir.
x(), dir.
y(), dir.
z(), dir.
mag());
444 p = fastjet::PseudoJet(it->p4().px(), it->p4().py(), it->p4().pz(), it->p4().energy());
446 p.set_user_info(
new VertexInfo(it - extSecVertex->begin()));
447 fjInputs.push_back(p);
453 std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(
fjClusterSeq->inclusive_jets(
jetPtMin));
456 if (inclusiveJets.size() < fatJetsHandle->size())
458 <<
"There are fewer reclustered (" << inclusiveJets.size() <<
") than original fat jets ("
459 << fatJetsHandle->size()
460 <<
"). Please check that the jet algorithm and jet size match those used for the original jet collection.";
463 std::vector<int> reclusteredIndices;
464 matchReclusteredJets<edm::View<reco::Jet> >(fatJetsHandle, inclusiveJets, reclusteredIndices,
"fat");
467 std::vector<int> groomedIndices;
472 std::vector<std::vector<int> > subjetIndices;
474 matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
476 matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
479 for (
size_t i = 0;
i < fatJetsHandle->size(); ++
i) {
480 if (reclusteredIndices.at(
i) < 0)
483 if (fatJetsHandle->at(
i).pt() == 0)
486 <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
490 if (subjetIndices.at(
i).empty())
494 if ((
std::abs(inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - fatJetsHandle->at(
i).pt()) /
496 if (fatJetsHandle->at(
i).pt() < 10.)
498 <<
"The reclustered and original fat jet " <<
i <<
" have different Pt's ("
499 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << fatJetsHandle->at(
i).pt()
500 <<
" GeV, respectively).\n"
501 <<
"Please check that the jet algorithm and jet size match those used for the original fat jet "
502 "collection and also make sure the original fat jets are uncorrected. In addition, make sure you "
503 "are not using CaloJets which are presently not supported.\n"
504 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
505 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
506 "precision in which case make sure the original jet collection is produced and reclustering is "
507 "performed in the same job.";
510 <<
"The reclustered and original fat jet " <<
i <<
" have different Pt's ("
511 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << fatJetsHandle->at(
i).pt()
512 <<
" GeV, respectively).\n"
513 <<
"Please check that the jet algorithm and jet size match those used for the original fat jet "
514 "collection and also make sure the original fat jets are uncorrected. In addition, make sure you "
515 "are not using CaloJets which are presently not supported.\n"
516 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
517 "precision in which case make sure the original jet collection is produced and reclustering is "
518 "performed in the same job.";
522 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
524 std::vector<int> svIndices;
526 for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
528 if (!it->has_user_info())
531 svIndices.push_back(it->user_info<VertexInfo>().vertexIndex());
535 for (
size_t sv = 0; sv < svIndices.size(); ++sv) {
536 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
537 const VTX &extSV = (*extSecVertex)[svIndices.at(sv)];
540 fastjet::PseudoJet
p(
541 dir.
x(), dir.
y(), dir.
z(), dir.
mag());
543 p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
545 std::vector<double> dR2toSubjets;
547 for (
size_t sj = 0; sj < subjetIndices.at(
i).size(); ++sj)
550 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->rapidity(),
551 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->phi()));
554 int closestSubjetIdx =
555 std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
557 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back(svIndices.at(sv));
561 if (inclusiveJets.size() < trackIPTagInfos->size())
563 <<
"There are fewer reclustered (" << inclusiveJets.size() <<
") than original jets ("
564 << trackIPTagInfos->size()
565 <<
"). Please check that the jet algorithm and jet size match those used for the original jet collection.";
568 std::vector<int> reclusteredIndices;
569 matchReclusteredJets<std::vector<IPTI> >(
trackIPTagInfos, inclusiveJets, reclusteredIndices);
572 for (
size_t i = 0;
i < trackIPTagInfos->size(); ++
i) {
573 if (reclusteredIndices.at(
i) < 0)
576 if (trackIPTagInfos->at(
i).jet()->pt() == 0)
579 <<
"The original jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
584 if ((
std::abs(inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - trackIPTagInfos->at(
i).jet()->pt()) /
586 if (trackIPTagInfos->at(
i).jet()->pt() < 10.)
588 <<
"The reclustered and original jet " <<
i <<
" have different Pt's ("
589 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << trackIPTagInfos->at(
i).jet()->pt()
590 <<
" GeV, respectively).\n"
591 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection "
592 "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
593 "CaloJets which are presently not supported.\n"
594 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\n"
595 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
596 "precision in which case make sure the original jet collection is produced and reclustering is "
597 "performed in the same job.";
600 <<
"The reclustered and original jet " <<
i <<
" have different Pt's ("
601 << inclusiveJets.at(reclusteredIndices.at(
i)).
pt() <<
" vs " << trackIPTagInfos->at(
i).jet()->pt()
602 <<
" GeV, respectively).\n"
603 <<
"Please check that the jet algorithm and jet size match those used for the original jet collection "
604 "and also make sure the original jets are uncorrected. In addition, make sure you are not using "
605 "CaloJets which are presently not supported.\n"
606 <<
"\nIn extremely rare instances the mismatch could be caused by a difference in the machine "
607 "precision in which case make sure the original jet collection is produced and reclustering is "
608 "performed in the same job.";
612 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
615 for (std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end();
617 if (!it->has_user_info())
620 clusteredSVs.at(
i).push_back(it->user_info<VertexInfo>().vertexIndex());
628 std::vector<int> groomedIndices;
633 std::vector<std::vector<int> > subjetIndices;
635 matchSubjets(groomedIndices, groomedFatJetsHandle, trackIPTagInfos, subjetIndices);
637 matchSubjets(fatJetsHandle, trackIPTagInfos, subjetIndices);
640 for (
size_t i = 0;
i < fatJetsHandle->size(); ++
i) {
641 if (fatJetsHandle->at(
i).pt() == 0)
644 <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
648 if (subjetIndices.at(
i).empty())
654 size_t sv = (it - extSecVertex->begin());
656 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
657 const VTX &extSV = (*extSecVertex)[sv];
659 GlobalVector jetDir(fatJetsHandle->at(
i).px(), fatJetsHandle->at(
i).py(), fatJetsHandle->at(
i).pz());
665 fastjet::PseudoJet
p(
666 dir.
x(), dir.
y(), dir.
z(), dir.
mag());
668 p = fastjet::PseudoJet(extSV.p4().px(), extSV.p4().py(), extSV.p4().pz(), extSV.p4().energy());
670 std::vector<double> dR2toSubjets;
672 for (
size_t sj = 0; sj < subjetIndices.at(
i).size(); ++sj)
675 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->rapidity(),
676 trackIPTagInfos->at(subjetIndices.at(
i).at(sj)).
jet()->phi()));
679 int closestSubjetIdx =
680 std::distance(dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()));
682 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back(sv);
688 std::unique_ptr<ConfigurableVertexReconstructor>
vertexReco;
689 std::unique_ptr<GhostTrackVertexFinder> vertexRecoGT;
691 vertexRecoGT = std::make_unique<GhostTrackVertexFinder>(
698 vertexReco = std::make_unique<ConfigurableVertexReconstructor>(
vtxRecoPSet);
700 TransientTrackMap primariesMap;
704 auto tagInfos = std::make_unique<Product>();
706 for (
typename std::vector<IPTI>::const_iterator iterJets = trackIPTagInfos->begin();
707 iterJets != trackIPTagInfos->end();
712 const Vertex &pv = *iterJets->primaryVertex();
714 std::set<TransientTrack> primaries;
717 TransientTrackMap::iterator pos = primariesMap.lower_bound(iter->get());
719 if (pos != primariesMap.end() && pos->first == iter->get())
720 primaries.insert(pos->second);
723 primariesMap.insert(pos, std::make_pair(iter->get(),
track));
724 primaries.insert(track);
731 GlobalVector jetDir(jetRef->momentum().x(), jetRef->momentum().y(), jetRef->momentum().z());
737 const std::vector<reco::btag::TrackIPData> &ipData = iterJets->impactParameterData();
741 std::vector<TransientTrack> fitTracks;
742 std::vector<GhostTrackState> gtStates;
743 std::unique_ptr<GhostTrackPrediction> gtPred;
745 gtPred = std::make_unique<GhostTrackPrediction>(*iterJets->ghostTrack());
747 for (
unsigned int i = 0;
i < indices.size();
i++) {
753 trackData.back().first = indices[
i];
763 TransientTrackMap::const_iterator pos = primariesMap.find(
reco::btag::toTrack((trackRef)));
765 if (pos != primariesMap.end()) {
766 primaries.erase(pos->second);
767 fitTrack = pos->second;
769 fitTrack = trackBuilder->build(trackRef);
770 fitTracks.push_back(fitTrack);
776 GlobalPoint pos = ipData[indices[
i]].closestToGhostTrack;
777 gtState.linearize(*gtPred,
true, gtPred->lambda(pos));
778 gtState.setWeight(ipData[indices[i]].ghostTrackWeight);
779 gtStates.push_back(gtState);
783 std::unique_ptr<GhostTrack> ghostTrack;
785 ghostTrack = std::make_unique<GhostTrack>(
789 GlobalVector(iterJets->ghostTrack()->px(), iterJets->ghostTrack()->py(), iterJets->ghostTrack()->pz()),
793 iterJets->ghostTrack()->chi2(),
794 iterJets->ghostTrack()->ndof());
798 std::vector<VTX> extAssoCollection;
799 std::vector<TransientVertex> fittedSVs;
800 std::vector<SecondaryVertex> SVs;
805 fittedSVs = vertexRecoGT->vertices(pv, *ghostTrack);
807 fittedSVs = vertexReco->vertices(fitTracks);
812 fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, *ghostTrack);
814 fittedSVs = vertexReco->vertices(fitTracks, *beamSpot);
821 for (
unsigned int i = 0; i < 7; i++) {
822 unsigned int covSrc = bsCovSrc[
i];
823 for (
unsigned int j = 0;
j < 7;
j++) {
825 if (!covSrc || bsCovSrc[
j] != covSrc)
827 else if (covSrc == 1)
828 v = beamSpot->covariance(i,
j);
829 else if (
j < 3 && i < 3)
837 beamSpot.
isValid() ? beamSpot->dxdz() : 0.,
838 beamSpot.
isValid() ? beamSpot->dydz() : 0.,
844 fittedSVs = vertexRecoGT->vertices(pv,
bs, *ghostTrack);
846 fittedSVs = vertexReco->vertices(fitTracks,
bs);
850 std::vector<TransientTrack> primaries_(primaries.begin(), primaries.end());
852 fittedSVs = vertexRecoGT->vertices(pv, *beamSpot, primaries_, *ghostTrack);
854 fittedSVs = vertexReco->vertices(primaries_, fitTracks, *beamSpot);
859 std::remove_copy_if(boost::make_transform_iterator(fittedSVs.begin(), svBuilder),
860 boost::make_transform_iterator(fittedSVs.end(), svBuilder),
861 std::back_inserter(SVs),
866 size_t jetIdx = (iterJets - trackIPTagInfos->begin());
868 for (
size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++) {
869 const VTX &extVertex = (*extSecVertex)[clusteredSVs.at(jetIdx).at(iExtSv)];
870 if (extVertex.p4().M() < 0.3)
872 extAssoCollection.push_back(extVertex);
875 for (
size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++) {
876 const VTX &extVertex = (*extSecVertex)[iExtSv];
879 extVertex.p4().M() < 0.3)
881 extAssoCollection.push_back(extVertex);
886 std::remove_copy_if(boost::make_transform_iterator(extAssoCollection.begin(), svBuilder),
887 boost::make_transform_iterator(extAssoCollection.end(), svBuilder),
888 std::back_inserter(SVs),
897 extAssoCollection.clear();
903 std::vector<typename TemplatedSecondaryVertexTagInfo<IPTI, VTX>::VertexData> svData;
905 svData.resize(vtxIndices.size());
906 for (
unsigned int idx = 0; idx < vtxIndices.size(); idx++) {
909 svData[idx].vertex = sv;
911 svData[idx].dist2d = sv.
dist2d();
912 svData[idx].dist3d = sv.
dist3d();
math::Error< dimension >::type CovarianceMatrix
reco::Vertex::Point convertPos(const GlobalPoint &p)
edm::EDGetTokenT< edm::View< reco::Jet > > token_fatJets
edm::ESGetToken< TransientTrackBuilder, TransientTrackRecord > token_trackBuilder
TrackSelector trackSelector
bool getByToken(EDGetToken token, Handle< PROD > &result) const
ClusterSequencePtr fjClusterSeq
edm::EDGetTokenT< reco::BeamSpot > token_BeamSpot
reco::Vertex::Error convertError(const GlobalError &ge)
constexpr bool isUninitialized() const noexcept
Measurement1D dist3d() 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
Log< level::Error, false > LogError
const reco::Track * toTrack(const reco::TrackBaseRef &t)
const Point & position() const
position
static GhostTrackVertexFinder::FitType getGhostTrackFitType(const std::string &name)
ConstraintType constraint
edm::EDGetTokenT< std::vector< IPTI > > token_trackIPTagInfo
IPTI::input_container::value_type input_item
reco::btag::IndexedTrackData IndexedTrackData
GlobalVector flightDirection(const reco::Vertex &pv, const reco::Vertex &sv)
edm::EDGetTokenT< edm::ValueMap< float > > token_weights
Abs< T >::type abs(const T &t)
trackRef_iterator tracks_end() const
last iterator over tracks
Measurement1D dist1d() const
VertexSorting< SecondaryVertex > vertexSorting
trackRef_iterator tracks_begin() const
first iterator over tracks
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)
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
T getParameter(std::string const &) const
Error error() const
return SMatrix
std::vector< reco::btag::IndexedTrackData > TrackDataVector
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
static int position[264][3]
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Log< level::Warning, false > LogWarning
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.5
