330 typedef std::map<const Track *, TransientTrack> TransientTrackMap;
354 if( groomedFatJetsHandle->size() > fatJetsHandle->size() )
355 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.";
360 unsigned int bsCovSrc[7] = { 0, };
361 double sigmaZ = 0.0, beamWidth = 0.0;
365 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = bsCovSrc[6] = 1;
366 sigmaZ = beamSpot->
sigmaZ();
372 bsCovSrc[0] = bsCovSrc[1] = 2;
373 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = 1;
374 sigmaZ = beamSpot->
sigmaZ();
378 bsCovSrc[0] = bsCovSrc[1] = bsCovSrc[2] = 2;
391 std::vector<std::vector<int> > clusteredSVs(trackIPTagInfos->size(),std::vector<int>());
395 std::vector<fastjet::PseudoJet> fjInputs;
401 std::vector<edm::Ptr<reco::Candidate> > constituents = it->getJetConstituents();
402 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
403 for( m = constituents.begin(); m != constituents.end(); ++
m )
406 if(constit->
pt() == 0)
408 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
411 fjInputs.push_back(fastjet::PseudoJet(constit->
px(),constit->
py(),constit->
pz(),constit->
energy()));
417 for(
typename std::vector<IPTI>::const_iterator it = trackIPTagInfos->begin(); it != trackIPTagInfos->end(); ++it)
419 std::vector<edm::Ptr<reco::Candidate> > constituents = it->jet()->getJetConstituents();
420 std::vector<edm::Ptr<reco::Candidate> >::const_iterator
m;
421 for( m = constituents.begin(); m != constituents.end(); ++
m )
424 if(constit->
pt() == 0)
426 edm::LogWarning(
"NullTransverseMomentum") <<
"dropping input candidate with pt=0";
429 fjInputs.push_back(fastjet::PseudoJet(constit->
px(),constit->
py(),constit->
pz(),constit->
energy()));
436 const reco::Vertex &
pv = *(trackIPTagInfos->front().primaryVertex());
439 fastjet::PseudoJet
p(dir.
x(),dir.
y(),dir.
z(),dir.
mag());
441 p = fastjet::PseudoJet(it->p4().px(),it->p4().py(),it->p4().pz(),it->p4().energy());
443 p.set_user_info(
new VertexInfo( it - extSecVertex->begin() ));
444 fjInputs.push_back(p);
450 std::vector<fastjet::PseudoJet> inclusiveJets = fastjet::sorted_by_pt(
fjClusterSeq->inclusive_jets(
jetPtMin) );
454 if( inclusiveJets.size() < fatJetsHandle->size() )
455 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.";
458 std::vector<int> reclusteredIndices;
459 matchReclusteredJets<edm::View<reco::Jet> >(fatJetsHandle,inclusiveJets,reclusteredIndices,
"fat");
462 std::vector<int> groomedIndices;
467 std::vector<std::vector<int> > subjetIndices;
469 matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices);
471 matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices);
474 for(
size_t i=0;
i<fatJetsHandle->size(); ++
i)
476 if( reclusteredIndices.at(
i) < 0 )
continue;
478 if( fatJetsHandle->at(
i).pt() == 0 )
480 edm::LogWarning(
"NullTransverseMomentum") <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
484 if( subjetIndices.at(
i).size()==0 )
continue;
487 if( (
std::abs( inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - fatJetsHandle->at(
i).pt() ) / fatJetsHandle->at(
i).pt() ) >
relPtTolerance )
489 if( fatJetsHandle->at(
i).pt() < 10. )
490 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" 491 <<
"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" 492 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\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.";
495 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" 496 <<
"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" 497 <<
"\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.";
501 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
503 std::vector<int> svIndices;
505 for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
507 if( !it->has_user_info() )
continue;
509 svIndices.push_back( it->user_info<VertexInfo>().vertexIndex() );
513 for(
size_t sv=0; sv<svIndices.size(); ++sv)
515 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
516 const VTX &extSV = (*extSecVertex)[ svIndices.at(sv) ];
519 fastjet::PseudoJet
p(dir.
x(),dir.
y(),dir.
z(),dir.
mag());
521 p = fastjet::PseudoJet(extSV.p4().px(),extSV.p4().py(),extSV.p4().pz(),extSV.p4().energy());
523 std::vector<double> dR2toSubjets;
525 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
526 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() ) );
529 int closestSubjetIdx =
std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
531 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back( svIndices.at(sv) );
537 if( inclusiveJets.size() < trackIPTagInfos->size() )
538 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.";
541 std::vector<int> reclusteredIndices;
542 matchReclusteredJets<std::vector<IPTI> >(
trackIPTagInfos,inclusiveJets,reclusteredIndices);
545 for(
size_t i=0;
i<trackIPTagInfos->size(); ++
i)
547 if( reclusteredIndices.at(
i) < 0 )
continue;
549 if( trackIPTagInfos->at(
i).jet()->pt() == 0 )
551 edm::LogWarning(
"NullTransverseMomentum") <<
"The original jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
556 if( (
std::abs( inclusiveJets.at(reclusteredIndices.at(
i)).
pt() - trackIPTagInfos->at(
i).jet()->pt() ) / trackIPTagInfos->at(
i).jet()->pt() ) >
relPtTolerance )
558 if( trackIPTagInfos->at(
i).jet()->pt() < 10. )
559 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" 560 <<
"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" 561 <<
"Since the mismatch is at low Pt, it is ignored and only a warning is issued.\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.";
564 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" 565 <<
"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" 566 <<
"\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.";
570 std::vector<fastjet::PseudoJet> constituents = inclusiveJets.at(reclusteredIndices.at(
i)).constituents();
573 for(std::vector<fastjet::PseudoJet>::const_iterator it = constituents.begin(); it != constituents.end(); ++it)
575 if( !it->has_user_info() )
continue;
577 clusteredSVs.at(
i).push_back( it->user_info<VertexInfo>().vertexIndex() );
586 std::vector<int> groomedIndices;
591 std::vector<std::vector<int> > subjetIndices;
593 matchSubjets(groomedIndices,groomedFatJetsHandle,trackIPTagInfos,subjetIndices);
595 matchSubjets(fatJetsHandle,trackIPTagInfos,subjetIndices);
598 for(
size_t i=0;
i<fatJetsHandle->size(); ++
i)
600 if( fatJetsHandle->at(
i).pt() == 0 )
602 edm::LogWarning(
"NullTransverseMomentum") <<
"The original fat jet " <<
i <<
" has Pt=0. This is not expected so the jet will be skipped.";
606 if( subjetIndices.at(
i).size()==0 )
continue;
612 size_t sv = ( it - extSecVertex->begin() );
614 const reco::Vertex &pv = *(trackIPTagInfos->front().primaryVertex());
615 const VTX &extSV = (*extSecVertex)[sv];
618 fatJetsHandle->at(
i).py(),
619 fatJetsHandle->at(
i).pz());
625 fastjet::PseudoJet
p(dir.
x(),dir.
y(),dir.
z(),dir.
mag());
627 p = fastjet::PseudoJet(extSV.p4().px(),extSV.p4().py(),extSV.p4().pz(),extSV.p4().energy());
629 std::vector<double> dR2toSubjets;
631 for(
size_t sj=0; sj<subjetIndices.at(
i).size(); ++sj)
632 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() ) );
635 int closestSubjetIdx =
std::distance( dR2toSubjets.begin(), std::min_element(dR2toSubjets.begin(), dR2toSubjets.end()) );
637 clusteredSVs.at(subjetIndices.at(
i).at(closestSubjetIdx)).push_back(sv);
643 std::auto_ptr<ConfigurableVertexReconstructor>
vertexReco;
644 std::auto_ptr<GhostTrackVertexFinder> vertexRecoGT;
656 TransientTrackMap primariesMap;
660 auto tagInfos = std::make_unique<Product>();
662 for(
typename std::vector<IPTI>::const_iterator iterJets =
663 trackIPTagInfos->begin(); iterJets != trackIPTagInfos->end();
668 const Vertex &pv = *iterJets->primaryVertex();
670 std::set<TransientTrack> primaries;
674 TransientTrackMap::iterator
pos =
675 primariesMap.lower_bound(iter->get());
677 if (pos != primariesMap.end() &&
678 pos->first == iter->get())
679 primaries.insert(pos->second);
684 primariesMap.insert(pos,
685 std::make_pair(iter->get(),
track));
686 primaries.insert(track);
694 jetRef->momentum().y(),
695 jetRef->momentum().z());
697 std::vector<std::size_t> indices =
702 const std::vector<reco::btag::TrackIPData> &ipData =
703 iterJets->impactParameterData();
707 std::vector<TransientTrack> fitTracks;
708 std::vector<GhostTrackState> gtStates;
709 std::auto_ptr<GhostTrackPrediction> gtPred;
712 *iterJets->ghostTrack()));
714 for(
unsigned int i = 0;
i < indices.size();
i++) {
720 trackData.back().first = indices[
i];
727 trackData.back().second.svStatus =
732 TransientTrackMap::const_iterator pos =
735 if (pos != primariesMap.end()) {
736 primaries.erase(pos->second);
737 fitTrack = pos->second;
739 fitTrack = trackBuilder->
build(trackRef);
740 fitTracks.push_back(fitTrack);
742 trackData.back().second.svStatus =
748 ipData[indices[
i]].closestToGhostTrack;
749 gtState.linearize(*gtPred,
true,
750 gtPred->lambda(pos));
751 gtState.setWeight(ipData[indices[i]].ghostTrackWeight);
752 gtStates.push_back(gtState);
756 std::auto_ptr<GhostTrack> ghostTrack;
763 iterJets->ghostTrack()->px(),
764 iterJets->ghostTrack()->py(),
765 iterJets->ghostTrack()->pz()),
768 iterJets->ghostTrack()->chi2(),
769 iterJets->ghostTrack()->ndof()));
774 std::vector<VTX> extAssoCollection;
775 std::vector<TransientVertex> fittedSVs;
776 std::vector<SecondaryVertex> SVs;
781 fittedSVs = vertexRecoGT->vertices(
784 fittedSVs = vertexReco->vertices(fitTracks);
789 fittedSVs = vertexRecoGT->vertices(
790 pv, *beamSpot, *ghostTrack);
792 fittedSVs = vertexReco->vertices(fitTracks,
800 for(
unsigned int i = 0; i < 7; i++) {
801 unsigned int covSrc = bsCovSrc[
i];
802 for(
unsigned int j = 0; j < 7; j++) {
804 if (!covSrc || bsCovSrc[j] != covSrc)
806 else if (covSrc == 1)
821 fittedSVs = vertexRecoGT->vertices(
822 pv, bs, *ghostTrack);
824 fittedSVs = vertexReco->vertices(fitTracks, bs);
828 std::vector<TransientTrack> primaries_(
829 primaries.begin(), primaries.end());
831 fittedSVs = vertexRecoGT->vertices(
832 pv, *beamSpot, primaries_,
835 fittedSVs = vertexReco->vertices(
836 primaries_, fitTracks,
842 std::remove_copy_if(boost::make_transform_iterator(
843 fittedSVs.begin(), svBuilder),
844 boost::make_transform_iterator(
845 fittedSVs.end(), svBuilder),
846 std::back_inserter(SVs),
851 size_t jetIdx = ( iterJets - trackIPTagInfos->begin() );
853 for(
size_t iExtSv = 0; iExtSv < clusteredSVs.at(jetIdx).size(); iExtSv++){
854 const VTX & extVertex = (*extSecVertex)[ clusteredSVs.at(jetIdx).at(iExtSv) ];
855 if( extVertex.p4().M() < 0.3 )
857 extAssoCollection.push_back( extVertex );
861 for(
size_t iExtSv = 0; iExtSv < extSecVertex->size(); iExtSv++){
862 const VTX & extVertex = (*extSecVertex)[iExtSv];
865 extAssoCollection.push_back( extVertex );
870 std::remove_copy_if(boost::make_transform_iterator( extAssoCollection.begin(), svBuilder),
871 boost::make_transform_iterator(extAssoCollection.end(), svBuilder),
872 std::back_inserter(SVs),
881 extAssoCollection.clear();
887 std::vector<typename TemplatedSecondaryVertexTagInfo<IPTI,VTX>::VertexData> svData;
889 svData.resize(vtxIndices.size());
890 for(
unsigned int idx = 0;
idx < vtxIndices.size();
idx++) {
893 svData[
idx].vertex = sv;
906 trackData, svData, SVs.size(),
908 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
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
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
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 Member Functions inherited from edm::stream::EDProducer<>
1.8.11
