128 if (recoBeamSpotHandle.
isValid()) {
131 edm::LogError(
"UnusableBeamSpot") <<
"No beam spot available from EventSetup";
136 if ((beamVertexState.error().cxx() <= 0.) || (beamVertexState.error().cyy() <= 0.) ||
137 (beamVertexState.error().czz() <= 0.)) {
139 edm::LogError(
"UnusableBeamSpot") <<
"Beamspot with invalid errors " << beamVertexState.error().matrix();
146 for (
auto const& old : oldVertices) {
147 if (!(old.isFake())) {
150 auto result = std::make_unique<reco::VertexCollection>();
165 std::vector<reco::TransientTrack> t_tks;
174 t_tks = (*theB).build(tks,
beamSpot);
177 std::cout <<
"RecoVertex/PrimaryVertexProducer" 178 <<
"Found: " << t_tks.size() <<
" reconstructed tracks" 189 std::cout <<
" clustering returned " <<
clusters.size() <<
" clusters from " << seltks.size()
190 <<
" selected tracks" << std::endl;
195 auto result = std::make_unique<reco::VertexCollection>();
198 std::vector<TransientVertex>
pvs;
199 for (
std::vector<std::vector<reco::TransientTrack>>::const_iterator iclus =
clusters.begin();
205 double meantime = 0.;
208 for (
const auto& tk : *iclus) {
209 const double time = tk.timeExt();
210 const double err = tk.dtErrorExt();
211 const double inverr =
err > 0. ? 1.0 /
err : 0.;
212 const double w = inverr * inverr;
217 meantime = sumwt / sumw;
218 double sumsq = sumwt2 - sumwt * sumwt / sumw;
219 double chisq = iclus->size() > 1 ? sumsq / double(iclus->size() - 1) : sumsq /
double(iclus->size());
220 vartime = chisq / sumw;
225 if (
algorithm->useBeamConstraint && validBS && (iclus->size() > 1)) {
227 }
else if (!(
algorithm->useBeamConstraint) && (iclus->size() > 1)) {
231 std::vector<std::pair<GlobalPoint, GlobalPoint>>
points;
232 if (
algorithm->useBeamConstraint && validBS && (iclus->size() > 1)) {
233 for (
const auto& itrack : *iclus) {
234 GlobalPoint p = itrack.stateAtBeamLine().trackStateAtPCA().position();
235 GlobalPoint err(itrack.stateAtBeamLine().transverseImpactParameter().error(),
236 itrack.stateAtBeamLine().transverseImpactParameter().error(),
237 itrack.track().dzError());
238 std::pair<GlobalPoint, GlobalPoint>
p2(
p,
err);
245 }
else if (!(
algorithm->useBeamConstraint) && (iclus->size() > 1)) {
246 for (
const auto& itrack : *iclus) {
247 GlobalPoint p = itrack.impactPointState().globalPosition();
248 GlobalPoint err(itrack.track().dxyError(), itrack.track().dxyError(), itrack.track().dzError());
249 std::pair<GlobalPoint, GlobalPoint>
p2(
p,
err);
259 "PrimaryVertexProducer: Something went wrong. You are not using the weighted mean fit and no algorithm was " 263 if (
f4D and
v.isValid()) {
264 auto err =
v.positionError().matrix4D();
266 auto trkWeightMap3d =
v.weightMap();
267 v =
TransientVertex(
v.position(), meantime,
err,
v.originalTracks(),
v.totalChiSquared(),
v.degreesOfFreedom());
268 v.weightMap(trkWeightMap3d);
276 std::cout <<
"=" <<
v.position().x() <<
" " <<
v.position().y() <<
" " <<
v.position().z();
279 std::cout <<
" cluster size = " << (*iclus).size() << std::endl;
281 std::cout <<
"Invalid fitted vertex, cluster size=" << (*iclus).size() << std::endl;
287 (!validBS || (*(
algorithm->vertexSelector))(
v, beamVertexState)))
292 std::cout <<
"PrimaryVertexProducerAlgorithm::vertices candidates =" <<
pvs.size() << std::endl;
297 <<
"more than half of candidate vertices lost " <<
pvs.size() <<
' ' <<
clusters.size();
299 if (
pvs.empty() && seltks.size() > 5)
301 <<
"no vertex found with " << seltks.size() <<
" tracks and " <<
clusters.size() <<
" vertex-candidates";
304 if (
pvs.size() > 1) {
309 for (std::vector<TransientVertex>::const_iterator
iv =
pvs.begin();
iv !=
pvs.end();
iv++) {
316 if ((bse.cxx() <= 0.) || (bse.cyy() <= 0.) || (bse.czz() <= 0.)) {
323 std::cout <<
"RecoVertex/PrimaryVertexProducer: " 324 <<
"Beamspot with invalid errors " << bse.matrix() << std::endl;
325 std::cout <<
"Will put Vertex derived from dummy-fake BeamSpot into Event.\n";
330 std::cout <<
"RecoVertex/PrimaryVertexProducer: " 331 <<
" will put Vertex derived from BeamSpot into Event.\n";
338 for (reco::VertexCollection::const_iterator
v = vColl.begin();
v != vColl.end(); ++
v) {
339 std::cout <<
"recvtx " << ivtx++ <<
"#trk " << std::setw(3) <<
v->tracksSize() <<
" chi2 " << std::setw(4)
340 <<
v->chi2() <<
" ndof " << std::setw(3) <<
v->ndof() <<
" x " << std::setw(6) <<
v->position().x()
341 <<
" dx " << std::setw(6) <<
v->xError() <<
" y " << std::setw(6) <<
v->position().y() <<
" dy " 342 << std::setw(6) <<
v->yError() <<
" z " << std::setw(6) <<
v->position().z() <<
" dz " << std::setw(6)
345 std::cout <<
" t " << std::setw(6) <<
v->t() <<
" dt " << std::setw(6) <<
v->tError();
edm::EDGetTokenT< reco::BeamSpot > bsToken
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
const edm::ESGetToken< TransientTrackBuilder, TransientTrackRecord > theTTBToken
std::vector< algo > algorithms
edm::EDGetTokenT< reco::VertexCollection > recoveryVtxToken
T const * product() const
edm::EDGetTokenT< reco::TrackCollection > trkToken
TrackClusterizerInZ * theTrackClusterizer
std::vector< Vertex > VertexCollection
collection of Vertex objects
Log< level::Error, false > LogError
edm::EDGetTokenT< edm::ValueMap< float > > trkTimesToken
edm::EDGetTokenT< edm::ValueMap< float > > trkTimeResosToken
virtual std::vector< std::vector< reco::TransientTrack > > clusterize(const std::vector< reco::TransientTrack > &tracks) const =0
TransientVertex weightedMeanOutlierRejectionBeamSpot(const std::vector< std::pair< GlobalPoint, GlobalPoint >> &points, std::vector< reco::TransientTrack > iclus, const reco::BeamSpot &beamSpot)
TrackFilterForPVFindingBase * theTrackFilter
virtual std::vector< reco::TransientTrack > select(const std::vector< reco::TransientTrack > &tracks) const =0
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
TransientVertex weightedMeanOutlierRejection(const std::vector< std::pair< GlobalPoint, GlobalPoint >> &points, std::vector< reco::TransientTrack > iclus)
constexpr float startError