Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:29:08

0001 #include <algorithm>
0002 #include <iterator>
0003 #include <memory>
0004 
0005 #include <map>
0006 #include <set>
0007 #include <vector>
0008 
0009 #include <Math/SVector.h>
0010 #include <Math/SMatrix.h>
0011 #include <Math/MatrixFunctions.h>
0012 
0013 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0014 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0015 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0016 
0017 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0018 
0019 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0021 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0022 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0023 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0024 #include "TrackingTools/TransientTrack/interface/TransientTrackFromFTSFactory.h"
0025 
0026 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0027 #include "RecoVertex/VertexPrimitives/interface/VertexFitter.h"
0028 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0029 #include "RecoVertex/VertexPrimitives/interface/CachingVertex.h"
0030 #include "RecoVertex/VertexPrimitives/interface/LinearizedTrackState.h"
0031 #include "RecoVertex/VertexPrimitives/interface/ConvertError.h"
0032 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0033 #include "RecoVertex/VertexTools/interface/LinearizedTrackStateFactory.h"
0034 #include "RecoVertex/VertexTools/interface/VertexTrackFactory.h"
0035 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0036 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
0037 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0038 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0039 
0040 #include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"
0041 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
0042 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackFitter.h"
0043 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
0044 #include "RecoVertex/GhostTrackFitter/interface/SequentialGhostTrackFitter.h"
0045 #include "RecoVertex/GhostTrackFitter/interface/KalmanGhostTrackUpdater.h"
0046 
0047 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackVertexFinder.h"
0048 
0049 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0050 
0051 // #define DEBUG
0052 
0053 #ifdef DEBUG
0054 #include "RecoBTag/SecondaryVertex/interface/TrackKinematics.h"
0055 #endif
0056 
0057 using namespace reco;
0058 
0059 namespace {
0060   using namespace ROOT::Math;
0061 
0062   typedef SVector<double, 3> Vector3;
0063 
0064   typedef SMatrix<double, 3, 3, MatRepSym<double, 3> > Matrix3S;
0065   typedef SMatrix<double, 2, 2, MatRepSym<double, 2> > Matrix2S;
0066   typedef SMatrix<double, 2, 3> Matrix23;
0067 
0068   typedef ReferenceCountingPointer<VertexTrack<5> > RefCountedVertexTrack;
0069   typedef ReferenceCountingPointer<LinearizedTrackState<5> > RefCountedLinearizedTrackState;
0070 
0071   struct VtxTrackIs {
0072     VtxTrackIs(const TransientTrack &track) : track(track) {}
0073     bool operator()(const RefCountedVertexTrack &vtxTrack) const {
0074       return vtxTrack->linearizedTrack()->track() == track;
0075     }
0076 
0077     const TransientTrack &track;
0078   };
0079 
0080   inline Vector3 conv(const GlobalVector &vec) {
0081     Vector3 result;
0082     result[0] = vec.x();
0083     result[1] = vec.y();
0084     result[2] = vec.z();
0085     return result;
0086   }
0087 
0088   inline double sqr(double arg) { return arg * arg; }
0089 }  // namespace
0090 
0091 struct GhostTrackVertexFinder::FinderInfo {
0092   FinderInfo(const CachingVertex<5> &primaryVertex,
0093              const GhostTrack &ghostTrack,
0094              const BeamSpot &beamSpot,
0095              bool hasBeamSpot,
0096              bool hasPrimaries)
0097       : primaryVertex(primaryVertex),
0098         pred(ghostTrack.prediction()),
0099         prior(ghostTrack.prior()),
0100         states(ghostTrack.states()),
0101         beamSpot(beamSpot),
0102         hasBeamSpot(hasBeamSpot),
0103         hasPrimaries(hasPrimaries),
0104         field(nullptr) {}
0105 
0106   const CachingVertex<5> &primaryVertex;
0107   GhostTrackPrediction pred;
0108   GhostTrackPrediction prior;
0109   std::vector<GhostTrackState> states;
0110   VertexState beamSpot;
0111   bool hasBeamSpot;
0112   bool hasPrimaries;
0113   const MagneticField *field;
0114   TransientTrack ghostTrack;
0115 };
0116 
0117 static TransientTrack transientGhostTrack(const GhostTrackPrediction &pred, const MagneticField *field) {
0118   return TransientTrackFromFTSFactory().build(pred.fts(field));
0119 }
0120 
0121 static double vtxErrorLong(const GlobalError &error, const GlobalVector &dir) {
0122   return ROOT::Math::Similarity(conv(dir), error.matrix());
0123 }
0124 
0125 static GlobalPoint vtxMean(const GlobalPoint &p1, const GlobalError &e1, const GlobalPoint &p2, const GlobalError &e2) {
0126   GlobalVector diff = p2 - p1;
0127 
0128   if UNLIKELY (p1 == p2) {
0129     edm::LogWarning("GhostTrackVertexFinder") << "Averaging with self at " << p1;
0130     return p1;
0131   }
0132 
0133   double err1 = vtxErrorLong(e1, diff);
0134   double err2 = vtxErrorLong(e2, diff);
0135 
0136   double weight = err1 / (err1 + err2);
0137 
0138   return p1 + weight * diff;
0139 }
0140 
0141 static VertexState stateMean(const VertexState &v1, const VertexState &v2) {
0142   return VertexState(vtxMean(v1.position(), v1.error(), v2.position(), v2.error()), v1.error() + v2.error());
0143 }
0144 
0145 static bool covarianceUpdate(
0146     Matrix3S &cov, const Vector3 &residual, const Matrix3S &error, double &chi2, double theta, double phi) {
0147   using namespace ROOT::Math;
0148 
0149   Matrix23 jacobian;
0150   jacobian(0, 0) = std::cos(phi) * std::cos(theta);
0151   jacobian(0, 1) = std::sin(phi) * std::cos(theta);
0152   jacobian(0, 2) = -std::sin(theta);
0153   jacobian(1, 0) = -std::sin(phi);
0154   jacobian(1, 1) = std::cos(phi);
0155 
0156   Matrix2S measErr = Similarity(jacobian, error);
0157   Matrix2S combErr = Similarity(jacobian, cov) + measErr;
0158   if (!measErr.Invert() || !combErr.Invert())
0159     return false;
0160 
0161   cov -= SimilarityT(jacobian * cov, combErr);
0162   chi2 += Similarity(jacobian * residual, measErr);
0163 
0164   return true;
0165 }
0166 
0167 static CachingVertex<5> vertexAtState(const TransientTrack &ghostTrack,
0168                                       const GhostTrackPrediction &pred,
0169                                       const GhostTrackState &state) {
0170   LinearizedTrackStateFactory linTrackFactory;
0171   VertexTrackFactory<5> vertexTrackFactory;
0172 
0173   if (!state.isValid())
0174     return CachingVertex<5>();
0175 
0176   GlobalPoint pca1 = pred.position(state.lambda());
0177   GlobalError err1 = pred.positionError(state.lambda());
0178 
0179   GlobalPoint pca2 = state.globalPosition();
0180   GlobalError err2 = state.cartesianError();
0181 
0182   GlobalPoint point = vtxMean(pca1, err1, pca2, err2);
0183 
0184   const TransientTrack &recTrack = state.track();
0185 
0186   RefCountedLinearizedTrackState linState[2] = {linTrackFactory.linearizedTrackState(point, ghostTrack),
0187                                                 linTrackFactory.linearizedTrackState(point, recTrack)};
0188   if (!linState[0]->isValid() || !linState[1]->isValid())
0189     return CachingVertex<5>();
0190 
0191   Matrix3S cov = SMatrixIdentity();
0192   cov *= 10000;
0193   double chi2 = 0.;
0194   if (!covarianceUpdate(cov,
0195                         conv(pca1 - point),
0196                         err1.matrix(),
0197                         chi2,
0198                         linState[0]->predictedStateParameters()[1],
0199                         linState[0]->predictedStateParameters()[2]) ||
0200       !covarianceUpdate(cov,
0201                         conv(pca2 - point),
0202                         err2.matrix(),
0203                         chi2,
0204                         linState[1]->predictedStateParameters()[1],
0205                         linState[1]->predictedStateParameters()[2]))
0206     return CachingVertex<5>();
0207 
0208   GlobalError error(cov);
0209   VertexState vtxState(point, error);
0210 
0211   std::vector<RefCountedVertexTrack> linTracks(2);
0212   linTracks[0] = vertexTrackFactory.vertexTrack(linState[0], vtxState);
0213   linTracks[1] = vertexTrackFactory.vertexTrack(linState[1], vtxState);
0214 
0215   return CachingVertex<5>(point, error, linTracks, chi2);
0216 }
0217 
0218 static RefCountedVertexTrack relinearizeTrack(const RefCountedVertexTrack &track,
0219                                               const VertexState &state,
0220                                               const VertexTrackFactory<5> factory) {
0221   RefCountedLinearizedTrackState linTrack = track->linearizedTrack();
0222   linTrack = linTrack->stateWithNewLinearizationPoint(state.position());
0223   return factory.vertexTrack(linTrack, state);
0224 }
0225 
0226 static std::vector<RefCountedVertexTrack> relinearizeTracks(const std::vector<RefCountedVertexTrack> &tracks,
0227                                                             const VertexState &state) {
0228   VertexTrackFactory<5> vertexTrackFactory;
0229 
0230   std::vector<RefCountedVertexTrack> finalTracks;
0231   finalTracks.reserve(tracks.size());
0232 
0233   for (std::vector<RefCountedVertexTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter)
0234     finalTracks.push_back(relinearizeTrack(*iter, state, vertexTrackFactory));
0235 
0236   return finalTracks;
0237 }
0238 
0239 double GhostTrackVertexFinder::vertexCompat(
0240     const CachingVertex<5> &vtx1, const CachingVertex<5> &vtx2, const FinderInfo &info, double scale1, double scale2) {
0241   Vector3 diff = conv(vtx2.position() - vtx1.position());
0242   Matrix3S cov = scale1 * vtx1.error().matrix() + scale2 * vtx2.error().matrix();
0243 
0244   return sqr(ROOT::Math::Dot(diff, diff)) / ROOT::Math::Similarity(cov, diff);
0245 }
0246 
0247 static double trackVertexCompat(const CachingVertex<5> &vtx, const RefCountedVertexTrack &vertexTrack) {
0248   using namespace ROOT::Math;
0249 
0250   TransientTrack track = vertexTrack->linearizedTrack()->track();
0251   GlobalPoint point = vtx.position();
0252   AnalyticalImpactPointExtrapolator extrap(track.field());
0253   TrajectoryStateOnSurface tsos = extrap.extrapolate(track.impactPointState(), point);
0254 
0255   if (!tsos.isValid())
0256     return 1.0e6;
0257 
0258   GlobalPoint point1 = vtx.position();
0259   GlobalPoint point2 = tsos.globalPosition();
0260   Vector3 dir = conv(point2 - point1);
0261   Matrix3S error = vtx.error().matrix() + tsos.cartesianError().matrix().Sub<Matrix3S>(0, 0);
0262   if (!error.Invert())
0263     return 1.0e6;
0264 
0265   return ROOT::Math::Similarity(error, dir);
0266 }
0267 
0268 GhostTrackVertexFinder::GhostTrackVertexFinder()
0269     : maxFitChi2_(10.0), mergeThreshold_(3.0), primcut_(2.0), seccut_(5.0), fitType_(kAlwaysWithGhostTrack) {}
0270 
0271 GhostTrackVertexFinder::GhostTrackVertexFinder(
0272     double maxFitChi2, double mergeThreshold, double primcut, double seccut, FitType fitType)
0273     : maxFitChi2_(maxFitChi2), mergeThreshold_(mergeThreshold), primcut_(primcut), seccut_(seccut), fitType_(fitType) {}
0274 
0275 GhostTrackVertexFinder::~GhostTrackVertexFinder() {}
0276 
0277 GhostTrackFitter &GhostTrackVertexFinder::ghostTrackFitter() const {
0278   if (!ghostTrackFitter_.get())
0279     ghostTrackFitter_ = std::make_unique<GhostTrackFitter>();
0280 
0281   return *ghostTrackFitter_;
0282 }
0283 
0284 VertexFitter<5> &GhostTrackVertexFinder::vertexFitter(bool primary) const {
0285   std::unique_ptr<VertexFitter<5> > *ptr = primary ? &primVertexFitter_ : &secVertexFitter_;
0286   if (!ptr->get())
0287     *ptr = std::make_unique<AdaptiveVertexFitter>(GeometricAnnealing(primary ? primcut_ : seccut_));
0288 
0289   return **ptr;
0290 }
0291 
0292 #ifdef DEBUG
0293 static void debugTrack(const TransientTrack &track, const TransientVertex *vtx = 0) {
0294   const FreeTrajectoryState &fts = track.initialFreeState();
0295   GlobalVector mom = fts.momentum();
0296   std::cout << "\tpt = " << mom.perp() << ", eta = " << mom.eta() << ", phi = " << mom.phi()
0297             << ", charge = " << fts.charge();
0298   if (vtx && vtx->hasTrackWeight())
0299     std::cout << ", w = " << vtx->trackWeight(track) << std::endl;
0300   else
0301     std::cout << std::endl;
0302 }
0303 
0304 static void debugTracks(const std::vector<TransientTrack> &tracks, const TransientVertex *vtx = 0) {
0305   TrackKinematics kin;
0306   for (std::vector<TransientTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter) {
0307     debugTrack(*iter, vtx);
0308     try {
0309       TrackBaseRef track = iter->trackBaseRef();
0310       kin.add(*track);
0311     } catch (exception const &e) {
0312       // ignore, yeah this is debugging, so shut up ^^
0313     }
0314   }
0315 
0316   std::cout << "mass = " << kin.vectorSum().M() << std::endl;
0317 }
0318 
0319 static void debugTracks(const std::vector<RefCountedVertexTrack> &tracks) {
0320   std::vector<TransientTrack> tracks_;
0321   for (std::vector<RefCountedVertexTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter) {
0322     std::cout << "+ " << (*iter)->linearizedTrack()->linearizationPoint() << std::endl;
0323     tracks_.push_back((*iter)->linearizedTrack()->track());
0324   }
0325   debugTracks(tracks_);
0326 }
0327 
0328 static void debugVertex(const TransientVertex &vtx, const GhostTrackPrediction &pred) {
0329   double origin = 0.;
0330   std::cout << vtx.position() << "-> " << vtx.totalChiSquared() << " with " << vtx.degreesOfFreedom() << std::endl;
0331 
0332   double err = std::sqrt(vtxErrorLong(vtx.positionError(), pred.direction()) / pred.rho2());
0333   std::cout << "* " << (pred.lambda(vtx.position()) * pred.rho() - origin) << " +-" << err << " => "
0334             << ((pred.lambda(vtx.position()) * pred.rho() - origin) / err) << std::endl;
0335 
0336   std::vector<TransientTrack> tracks = vtx.originalTracks();
0337   debugTracks(tracks, &vtx);
0338 }
0339 #endif
0340 
0341 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0342                                                               const GlobalVector &direction,
0343                                                               double coneRadius,
0344                                                               const std::vector<TransientTrack> &tracks) const {
0345   return vertices(RecoVertex::convertPos(primaryVertex.position()),
0346                   RecoVertex::convertError(primaryVertex.error()),
0347                   direction,
0348                   coneRadius,
0349                   tracks);
0350 }
0351 
0352 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0353                                                               const GlobalVector &direction,
0354                                                               double coneRadius,
0355                                                               const reco::BeamSpot &beamSpot,
0356                                                               const std::vector<TransientTrack> &tracks) const {
0357   return vertices(RecoVertex::convertPos(primaryVertex.position()),
0358                   RecoVertex::convertError(primaryVertex.error()),
0359                   direction,
0360                   coneRadius,
0361                   beamSpot,
0362                   tracks);
0363 }
0364 
0365 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0366                                                               const GlobalVector &direction,
0367                                                               double coneRadius,
0368                                                               const reco::BeamSpot &beamSpot,
0369                                                               const std::vector<TransientTrack> &primaries,
0370                                                               const std::vector<TransientTrack> &tracks) const {
0371   return vertices(RecoVertex::convertPos(primaryVertex.position()),
0372                   RecoVertex::convertError(primaryVertex.error()),
0373                   direction,
0374                   coneRadius,
0375                   beamSpot,
0376                   primaries,
0377                   tracks);
0378 }
0379 
0380 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0381                                                               const GlobalError &primaryError,
0382                                                               const GlobalVector &direction,
0383                                                               double coneRadius,
0384                                                               const std::vector<TransientTrack> &tracks) const {
0385   GhostTrack ghostTrack = ghostTrackFitter().fit(primaryPosition, primaryError, direction, coneRadius, tracks);
0386 
0387   CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0388 
0389   std::vector<TransientVertex> result = vertices(ghostTrack, primary);
0390 
0391 #ifdef DEBUG
0392   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0393     debugVertex(*iter, ghostTrack.prediction());
0394 #endif
0395 
0396   return result;
0397 }
0398 
0399 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0400                                                               const GlobalError &primaryError,
0401                                                               const GlobalVector &direction,
0402                                                               double coneRadius,
0403                                                               const BeamSpot &beamSpot,
0404                                                               const std::vector<TransientTrack> &tracks) const {
0405   GhostTrack ghostTrack = ghostTrackFitter().fit(primaryPosition, primaryError, direction, coneRadius, tracks);
0406 
0407   CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0408 
0409   std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true);
0410 
0411 #ifdef DEBUG
0412   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0413     debugVertex(*iter, ghostTrack.prediction());
0414 #endif
0415 
0416   return result;
0417 }
0418 
0419 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0420                                                               const GlobalError &primaryError,
0421                                                               const GlobalVector &direction,
0422                                                               double coneRadius,
0423                                                               const BeamSpot &beamSpot,
0424                                                               const std::vector<TransientTrack> &primaries,
0425                                                               const std::vector<TransientTrack> &tracks) const {
0426   GhostTrack ghostTrack = ghostTrackFitter().fit(primaryPosition, primaryError, direction, coneRadius, tracks);
0427 
0428   std::vector<RefCountedVertexTrack> primaryVertexTracks;
0429   if (!primaries.empty()) {
0430     LinearizedTrackStateFactory linTrackFactory;
0431     VertexTrackFactory<5> vertexTrackFactory;
0432 
0433     VertexState state(primaryPosition, primaryError);
0434 
0435     for (std::vector<TransientTrack>::const_iterator iter = primaries.begin(); iter != primaries.end(); ++iter) {
0436       RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(primaryPosition, *iter);
0437 
0438       primaryVertexTracks.push_back(vertexTrackFactory.vertexTrack(linState, state));
0439     }
0440   }
0441 
0442   CachingVertex<5> primary(primaryPosition, primaryError, primaryVertexTracks, 0.);
0443 
0444   std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true, true);
0445 
0446 #ifdef DEBUG
0447   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0448     debugVertex(*iter, ghostTrack.prediction());
0449 #endif
0450 
0451   return result;
0452 }
0453 
0454 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0455                                                               const GlobalError &primaryError,
0456                                                               const GhostTrack &ghostTrack) const {
0457   CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0458 
0459   std::vector<TransientVertex> result = vertices(ghostTrack, primary);
0460 
0461 #ifdef DEBUG
0462   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0463     debugVertex(*iter, ghostTrack.prediction());
0464 #endif
0465 
0466   return result;
0467 }
0468 
0469 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0470                                                               const GlobalError &primaryError,
0471                                                               const BeamSpot &beamSpot,
0472                                                               const GhostTrack &ghostTrack) const {
0473   CachingVertex<5> primary(primaryPosition, primaryError, std::vector<RefCountedVertexTrack>(), 0.);
0474 
0475   std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true);
0476 
0477 #ifdef DEBUG
0478   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0479     debugVertex(*iter, ghostTrack.prediction());
0480 #endif
0481 
0482   return result;
0483 }
0484 
0485 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GlobalPoint &primaryPosition,
0486                                                               const GlobalError &primaryError,
0487                                                               const BeamSpot &beamSpot,
0488                                                               const std::vector<TransientTrack> &primaries,
0489                                                               const GhostTrack &ghostTrack) const {
0490   std::vector<RefCountedVertexTrack> primaryVertexTracks;
0491   if (!primaries.empty()) {
0492     LinearizedTrackStateFactory linTrackFactory;
0493     VertexTrackFactory<5> vertexTrackFactory;
0494 
0495     VertexState state(primaryPosition, primaryError);
0496 
0497     for (std::vector<TransientTrack>::const_iterator iter = primaries.begin(); iter != primaries.end(); ++iter) {
0498       RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(primaryPosition, *iter);
0499 
0500       primaryVertexTracks.push_back(vertexTrackFactory.vertexTrack(linState, state));
0501     }
0502   }
0503 
0504   CachingVertex<5> primary(primaryPosition, primaryError, primaryVertexTracks, 0.);
0505 
0506   std::vector<TransientVertex> result = vertices(ghostTrack, primary, beamSpot, true, true);
0507 
0508 #ifdef DEBUG
0509   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0510     debugVertex(*iter, ghostTrack.prediction());
0511 #endif
0512 
0513   return result;
0514 }
0515 
0516 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0517                                                               const GhostTrack &ghostTrack) const {
0518   return vertices(
0519       RecoVertex::convertPos(primaryVertex.position()), RecoVertex::convertError(primaryVertex.error()), ghostTrack);
0520 }
0521 
0522 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0523                                                               const BeamSpot &beamSpot,
0524                                                               const GhostTrack &ghostTrack) const {
0525   return vertices(RecoVertex::convertPos(primaryVertex.position()),
0526                   RecoVertex::convertError(primaryVertex.error()),
0527                   beamSpot,
0528                   ghostTrack);
0529 }
0530 
0531 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0532                                                               const BeamSpot &beamSpot,
0533                                                               const std::vector<TransientTrack> &primaries,
0534                                                               const GhostTrack &ghostTrack) const {
0535   return vertices(RecoVertex::convertPos(primaryVertex.position()),
0536                   RecoVertex::convertError(primaryVertex.error()),
0537                   beamSpot,
0538                   primaries,
0539                   ghostTrack);
0540 }
0541 
0542 static GhostTrackPrediction dummyPrediction(const Vertex &primaryVertex, const Track &ghostTrack) {
0543   return GhostTrackPrediction(RecoVertex::convertPos(primaryVertex.position()),
0544                               RecoVertex::convertError(primaryVertex.error()),
0545                               GlobalVector(ghostTrack.px(), ghostTrack.py(), ghostTrack.pz()),
0546                               0.05);
0547 }
0548 
0549 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0550                                                               const Track &ghostTrack,
0551                                                               const std::vector<TransientTrack> &tracks,
0552                                                               const std::vector<float> &weights) const {
0553   GhostTrack ghostTrack_(ghostTrack,
0554                          tracks,
0555                          weights,
0556                          dummyPrediction(primaryVertex, ghostTrack),
0557                          RecoVertex::convertPos(primaryVertex.position()),
0558                          true);
0559 
0560   CachingVertex<5> primary(RecoVertex::convertPos(primaryVertex.position()),
0561                            RecoVertex::convertError(primaryVertex.error()),
0562                            std::vector<RefCountedVertexTrack>(),
0563                            0.);
0564 
0565   std::vector<TransientVertex> result = vertices(ghostTrack_, primary);
0566 
0567 #ifdef DEBUG
0568   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0569     debugVertex(*iter, ghostTrack_.prediction());
0570 #endif
0571 
0572   return result;
0573 }
0574 
0575 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0576                                                               const Track &ghostTrack,
0577                                                               const BeamSpot &beamSpot,
0578                                                               const std::vector<TransientTrack> &tracks,
0579                                                               const std::vector<float> &weights) const {
0580   GhostTrack ghostTrack_(ghostTrack,
0581                          tracks,
0582                          weights,
0583                          dummyPrediction(primaryVertex, ghostTrack),
0584                          RecoVertex::convertPos(primaryVertex.position()),
0585                          true);
0586 
0587   CachingVertex<5> primary(RecoVertex::convertPos(primaryVertex.position()),
0588                            RecoVertex::convertError(primaryVertex.error()),
0589                            std::vector<RefCountedVertexTrack>(),
0590                            0.);
0591 
0592   std::vector<TransientVertex> result = vertices(ghostTrack_, primary, beamSpot, true);
0593 
0594 #ifdef DEBUG
0595   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0596     debugVertex(*iter, ghostTrack_.prediction());
0597 #endif
0598 
0599   return result;
0600 }
0601 
0602 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const Vertex &primaryVertex,
0603                                                               const Track &ghostTrack,
0604                                                               const BeamSpot &beamSpot,
0605                                                               const std::vector<TransientTrack> &primaries,
0606                                                               const std::vector<TransientTrack> &tracks,
0607                                                               const std::vector<float> &weights) const {
0608   GhostTrack ghostTrack_(ghostTrack,
0609                          tracks,
0610                          weights,
0611                          dummyPrediction(primaryVertex, ghostTrack),
0612                          RecoVertex::convertPos(primaryVertex.position()),
0613                          true);
0614 
0615   std::vector<RefCountedVertexTrack> primaryVertexTracks;
0616   if (!primaries.empty()) {
0617     LinearizedTrackStateFactory linTrackFactory;
0618     VertexTrackFactory<5> vertexTrackFactory;
0619 
0620     VertexState state(RecoVertex::convertPos(primaryVertex.position()),
0621                       RecoVertex::convertError(primaryVertex.error()));
0622 
0623     for (std::vector<TransientTrack>::const_iterator iter = primaries.begin(); iter != primaries.end(); ++iter) {
0624       RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(state.position(), *iter);
0625 
0626       primaryVertexTracks.push_back(vertexTrackFactory.vertexTrack(linState, state));
0627     }
0628   }
0629 
0630   CachingVertex<5> primary(RecoVertex::convertPos(primaryVertex.position()),
0631                            RecoVertex::convertError(primaryVertex.error()),
0632                            primaryVertexTracks,
0633                            0.);
0634 
0635   std::vector<TransientVertex> result = vertices(ghostTrack_, primary, beamSpot, true, true);
0636 
0637 #ifdef DEBUG
0638   for (std::vector<TransientVertex>::const_iterator iter = result.begin(); iter != result.end(); ++iter)
0639     debugVertex(*iter, ghostTrack_.prediction());
0640 #endif
0641 
0642   return result;
0643 }
0644 
0645 static double fitChi2(const CachingVertex<5> &vtx) { return vtx.totalChiSquared() / vtx.degreesOfFreedom(); }
0646 
0647 std::vector<CachingVertex<5> > GhostTrackVertexFinder::initialVertices(const FinderInfo &info) const {
0648   std::vector<CachingVertex<5> > vertices;
0649   for (std::vector<GhostTrackState>::const_iterator iter = info.states.begin(); iter != info.states.end(); ++iter) {
0650     if (!iter->isValid())
0651       continue;
0652 
0653     GhostTrackState state(*iter);
0654     state.linearize(info.pred);
0655 
0656     if (!state.isValid())
0657       continue;
0658 
0659     CachingVertex<5> vtx = vertexAtState(info.ghostTrack, info.pred, state);
0660 
0661     if (vtx.isValid())  // && fitChi2(vtx) < maxFitChi2_)
0662       vertices.push_back(vtx);
0663   }
0664 
0665   return vertices;
0666 }
0667 
0668 static void mergeTrackHelper(const std::vector<RefCountedVertexTrack> &tracks,
0669                              std::vector<RefCountedVertexTrack> &newTracks,
0670                              const VertexState &state,
0671                              const VtxTrackIs &ghostTrackFinder,
0672                              RefCountedVertexTrack &ghostTrack,
0673                              const VertexTrackFactory<5> &factory) {
0674   for (std::vector<RefCountedVertexTrack>::const_iterator iter = tracks.begin(); iter != tracks.end(); ++iter) {
0675     bool gt = ghostTrackFinder(*iter);
0676     if (gt && ghostTrack)
0677       continue;
0678 
0679     RefCountedVertexTrack track = relinearizeTrack(*iter, state, factory);
0680 
0681     if (gt)
0682       ghostTrack = *iter;
0683     else
0684       newTracks.push_back(*iter);
0685   }
0686 }
0687 
0688 CachingVertex<5> GhostTrackVertexFinder::mergeVertices(const CachingVertex<5> &vertex1,
0689                                                        const CachingVertex<5> &vertex2,
0690                                                        const FinderInfo &info,
0691                                                        bool isPrimary) const {
0692   VertexTrackFactory<5> vertexTrackFactory;
0693 
0694   VertexState state = stateMean(vertex1.vertexState(), vertex2.vertexState());
0695   std::vector<RefCountedVertexTrack> linTracks;
0696   VtxTrackIs isGhostTrack(info.ghostTrack);
0697   RefCountedVertexTrack vtxGhostTrack;
0698 
0699   mergeTrackHelper(vertex1.tracks(), linTracks, state, isGhostTrack, vtxGhostTrack, vertexTrackFactory);
0700   mergeTrackHelper(vertex2.tracks(), linTracks, state, isGhostTrack, vtxGhostTrack, vertexTrackFactory);
0701 
0702   if (vtxGhostTrack && (fitType_ == kAlwaysWithGhostTrack || linTracks.size() < 2))
0703     linTracks.push_back(vertexTrackFactory.vertexTrack(vtxGhostTrack->linearizedTrack(), vtxGhostTrack->vertexState()));
0704 
0705   try {
0706     CachingVertex<5> vtx;
0707     if (info.hasBeamSpot && isPrimary)
0708       vtx = vertexFitter(true).vertex(linTracks, info.beamSpot.position(), info.beamSpot.error());
0709     else
0710       vtx = vertexFitter(isPrimary).vertex(linTracks);
0711     if (vtx.isValid())
0712       return vtx;
0713   } catch (const VertexException &e) {
0714     // fit failed
0715   }
0716 
0717   return CachingVertex<5>();
0718 }
0719 
0720 bool GhostTrackVertexFinder::recursiveMerge(std::vector<CachingVertex<5> > &vertices, const FinderInfo &info) const {
0721   typedef std::pair<unsigned int, unsigned int> Indices;
0722 
0723   std::multimap<float, Indices> compatMap;
0724   unsigned int n = vertices.size();
0725   for (unsigned int i = 0; i < n; i++) {
0726     const CachingVertex<5> &v1 = vertices[i];
0727     for (unsigned int j = i + 1; j < n; j++) {
0728       const CachingVertex<5> &v2 = vertices[j];
0729 
0730       float compat = vertexCompat(v1, v2, info, i == 0 ? primcut_ : seccut_, seccut_);
0731 
0732       if (compat > mergeThreshold_)
0733         continue;
0734 
0735       compatMap.insert(std::make_pair(compat, Indices(i, j)));
0736     }
0737   }
0738 
0739   bool changesMade = false;
0740   bool repeat = true;
0741   while (repeat) {
0742     repeat = false;
0743     for (std::multimap<float, Indices>::const_iterator iter = compatMap.begin(); iter != compatMap.end(); ++iter) {
0744       unsigned int v1 = iter->second.first;
0745       unsigned int v2 = iter->second.second;
0746 
0747       CachingVertex<5> newVtx = mergeVertices(vertices[v1], vertices[v2], info, v1 == 0);
0748       if (!newVtx.isValid() || (v1 != 0 && fitChi2(newVtx) > maxFitChi2_))
0749         continue;
0750 
0751       std::swap(vertices[v1], newVtx);
0752       vertices.erase(vertices.begin() + v2);
0753       n--;
0754 
0755       std::multimap<float, Indices> newCompatMap;
0756       for (++iter; iter != compatMap.end(); ++iter) {
0757         if (iter->second.first == v1 || iter->second.first == v2 || iter->second.second == v1 ||
0758             iter->second.second == v2)
0759           continue;
0760 
0761         Indices indices = iter->second;
0762         indices.first -= indices.first > v2;
0763         indices.second -= indices.second > v2;
0764 
0765         newCompatMap.insert(std::make_pair(iter->first, indices));
0766       }
0767       std::swap(compatMap, newCompatMap);
0768 
0769       for (unsigned int i = 0; i < n; i++) {
0770         if (i == v1)
0771           continue;
0772 
0773         const CachingVertex<5> &other = vertices[i];
0774         float compat =
0775             vertexCompat(vertices[v1], other, info, v1 == 0 ? primcut_ : seccut_, i == 0 ? primcut_ : seccut_);
0776 
0777         if (compat > mergeThreshold_)
0778           continue;
0779 
0780         compatMap.insert(std::make_pair(compat, Indices(std::min(i, v1), std::max(i, v1))));
0781       }
0782 
0783       changesMade = true;
0784       repeat = true;
0785       break;
0786     }
0787   }
0788 
0789   return changesMade;
0790 }
0791 
0792 bool GhostTrackVertexFinder::reassignTracks(std::vector<CachingVertex<5> > &vertices_, const FinderInfo &info) const {
0793   std::vector<std::pair<RefCountedVertexTrack, std::vector<RefCountedVertexTrack> > > trackBundles(vertices_.size());
0794 
0795   VtxTrackIs isGhostTrack(info.ghostTrack);
0796 
0797   bool assignmentChanged = false;
0798   for (std::vector<CachingVertex<5> >::const_iterator iter = vertices_.begin(); iter != vertices_.end(); ++iter) {
0799     std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
0800 
0801     if (vtxTracks.empty()) {
0802       LinearizedTrackStateFactory linTrackFactory;
0803       VertexTrackFactory<5> vertexTrackFactory;
0804 
0805       RefCountedLinearizedTrackState linState = linTrackFactory.linearizedTrackState(iter->position(), info.ghostTrack);
0806 
0807       trackBundles[iter - vertices_.begin()].first = vertexTrackFactory.vertexTrack(linState, iter->vertexState());
0808     }
0809 
0810     for (std::vector<RefCountedVertexTrack>::const_iterator track = vtxTracks.begin(); track != vtxTracks.end();
0811          ++track) {
0812       if (isGhostTrack(*track)) {
0813         trackBundles[iter - vertices_.begin()].first = *track;
0814         continue;
0815       }
0816 
0817       if ((*track)->weight() < 1e-3) {
0818         trackBundles[iter - vertices_.begin()].second.push_back(*track);
0819         continue;
0820       }
0821 
0822       unsigned int idx = iter - vertices_.begin();
0823       double best = 1.0e9;
0824       for (std::vector<CachingVertex<5> >::const_iterator vtx = vertices_.begin(); vtx != vertices_.end(); ++vtx) {
0825         if (info.pred.lambda(vtx->position()) < info.pred.lambda(vertices_[0].position()))
0826           continue;
0827 
0828         double compat = sqr(trackVertexCompat(*vtx, *track));
0829 
0830         compat /= (vtx == vertices_.begin()) ? primcut_ : seccut_;
0831 
0832         if (compat < best) {
0833           best = compat;
0834           idx = vtx - vertices_.begin();
0835         }
0836       }
0837 
0838       if ((int)idx != iter - vertices_.begin())
0839         assignmentChanged = true;
0840 
0841       trackBundles[idx].second.push_back(*track);
0842     }
0843   }
0844 
0845   if (!assignmentChanged)
0846     return false;
0847 
0848   VertexTrackFactory<5> vertexTrackFactory;
0849   std::vector<CachingVertex<5> > vertices;
0850   vertices.reserve(vertices_.size());
0851 
0852   for (std::vector<CachingVertex<5> >::const_iterator iter = vertices_.begin(); iter != vertices_.end(); ++iter) {
0853     const std::vector<RefCountedVertexTrack> &tracks = trackBundles[iter - vertices_.begin()].second;
0854     if (tracks.empty())
0855       continue;
0856 
0857     CachingVertex<5> vtx;
0858 
0859     if (tracks.size() == 1) {
0860       const TransientTrack &track = tracks[0]->linearizedTrack()->track();
0861 
0862       int idx = -1;
0863       for (std::vector<GhostTrackState>::const_iterator iter = info.states.begin(); iter != info.states.end(); ++iter) {
0864         if (iter->isTrack() && iter->track() == track) {
0865           idx = iter - info.states.begin();
0866           break;
0867         }
0868       }
0869       if (idx < 0)
0870         continue;
0871 
0872       vtx = vertexAtState(info.ghostTrack, info.pred, info.states[idx]);
0873       if (!vtx.isValid())
0874         continue;
0875     } else {
0876       std::vector<RefCountedVertexTrack> linTracks = relinearizeTracks(tracks, iter->vertexState());
0877 
0878       if (fitType_ == kAlwaysWithGhostTrack)
0879         linTracks.push_back(
0880             relinearizeTrack(trackBundles[iter - vertices_.begin()].first, iter->vertexState(), vertexTrackFactory));
0881 
0882       bool primary = iter == vertices_.begin();
0883       try {
0884         if (primary && info.hasBeamSpot)
0885           vtx = vertexFitter(true).vertex(linTracks, info.beamSpot.position(), info.beamSpot.error());
0886         else
0887           vtx = vertexFitter(primary).vertex(linTracks);
0888       } catch (const VertexException &e) {
0889         // fit failed;
0890       }
0891       if (!vtx.isValid())
0892         return false;
0893     }
0894 
0895     vertices.push_back(vtx);
0896   };
0897 
0898   std::swap(vertices_, vertices);
0899   return true;
0900 }
0901 
0902 void GhostTrackVertexFinder::refitGhostTrack(std::vector<CachingVertex<5> > &vertices, FinderInfo &info) const {
0903   VtxTrackIs isGhostTrack(info.ghostTrack);
0904 
0905   std::vector<GhostTrackState> states;
0906   std::vector<unsigned int> oldStates;
0907   oldStates.reserve(info.states.size());
0908 
0909   for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter) {
0910     std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
0911 
0912     oldStates.clear();
0913     for (std::vector<RefCountedVertexTrack>::const_iterator track = vtxTracks.begin(); track != vtxTracks.end();
0914          ++track) {
0915       if (isGhostTrack(*track) || (*track)->weight() < 1e-3)
0916         continue;
0917 
0918       const TransientTrack &tt = (*track)->linearizedTrack()->track();
0919 
0920       int idx = -1;
0921       for (std::vector<GhostTrackState>::const_iterator iter = info.states.begin(); iter != info.states.end(); ++iter) {
0922         if (iter->isTrack() && iter->track() == tt) {
0923           idx = iter - info.states.begin();
0924           break;
0925         }
0926       }
0927 
0928       if (idx >= 0)
0929         oldStates.push_back(idx);
0930     }
0931 
0932     if (oldStates.size() == 1)
0933       states.push_back(info.states[oldStates[0]]);
0934     else
0935       states.push_back(GhostTrackState(iter->vertexState()));
0936   }
0937 
0938   KalmanGhostTrackUpdater updater;
0939   SequentialGhostTrackFitter fitter;
0940   double ndof, chi2;
0941   info.pred = fitter.fit(updater, info.prior, states, ndof, chi2);
0942   TransientTrack ghostTrack = transientGhostTrack(info.pred, info.field);
0943 
0944   std::swap(info.states, states);
0945   states.clear();
0946 
0947   std::vector<CachingVertex<5> > newVertices;
0948   for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter) {
0949     std::vector<RefCountedVertexTrack> vtxTracks = iter->tracks();
0950 
0951     int idx = -1;
0952     bool redo = false;
0953     for (std::vector<RefCountedVertexTrack>::iterator track = vtxTracks.begin(); track != vtxTracks.end(); ++track) {
0954       if (isGhostTrack(*track)) {
0955         LinearizedTrackStateFactory linTrackFactory;
0956         VertexTrackFactory<5> vertexTrackFactory;
0957 
0958         *track = vertexTrackFactory.vertexTrack(linTrackFactory.linearizedTrackState(iter->position(), ghostTrack),
0959                                                 iter->vertexState());
0960         redo = true;
0961         continue;
0962       }
0963 
0964       const TransientTrack &tt = (*track)->linearizedTrack()->track();
0965 
0966       if (idx >= 0) {
0967         idx = -1;
0968         break;
0969       }
0970 
0971       for (std::vector<GhostTrackState>::const_iterator it = info.states.begin(); it != info.states.end(); ++it) {
0972         if (!it->isTrack())
0973           continue;
0974 
0975         if (it->track() == tt) {
0976           idx = it - info.states.begin();
0977           break;
0978         }
0979       }
0980     }
0981 
0982     if (idx >= 0) {
0983       CachingVertex<5> vtx = vertexAtState(ghostTrack, info.pred, info.states[idx]);
0984       if (vtx.isValid())
0985         newVertices.push_back(vtx);
0986     } else if (redo) {
0987       bool primary = iter == vertices.begin();
0988       CachingVertex<5> vtx;
0989       if (primary && info.hasBeamSpot)
0990         vtx = vertexFitter(true).vertex(vtxTracks, info.beamSpot.position(), info.beamSpot.error());
0991       else
0992         vtx = vertexFitter(primary).vertex(vtxTracks);
0993       if (vtx.isValid())
0994         newVertices.push_back(vtx);
0995     } else
0996       newVertices.push_back(*iter);
0997   }
0998 
0999   std::swap(newVertices, vertices);
1000   info.ghostTrack = ghostTrack;
1001 }
1002 
1003 // implementation
1004 
1005 std::vector<TransientVertex> GhostTrackVertexFinder::vertices(const GhostTrack &ghostTrack,
1006                                                               const CachingVertex<5> &primary,
1007                                                               const BeamSpot &beamSpot,
1008                                                               bool hasBeamSpot,
1009                                                               bool hasPrimaries) const {
1010   FinderInfo info(primary, ghostTrack, beamSpot, hasBeamSpot, hasPrimaries);
1011 
1012   std::vector<TransientVertex> result;
1013   if (info.states.empty())
1014     return result;
1015 
1016   info.field = info.states[0].track().field();
1017   info.ghostTrack = transientGhostTrack(info.pred, info.field);
1018 
1019   std::vector<CachingVertex<5> > vertices = initialVertices(info);
1020   if (primary.isValid()) {
1021     vertices.push_back(primary);
1022     if (vertices.size() > 1)
1023       std::swap(vertices.front(), vertices.back());
1024   }
1025 
1026   unsigned int reassigned = 0;
1027   while (reassigned < 3) {
1028     if (vertices.size() < 2)
1029       break;
1030 
1031 #ifdef DEBUG
1032     for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter)
1033 
1034       debugVertex(*iter, ghostTrack.prediction());
1035 
1036     std::cout << "----- recursive merging: ---------" << std::endl;
1037 #endif
1038 
1039     bool changed = recursiveMerge(vertices, info);
1040     if ((!changed && !reassigned) || vertices.size() < 2)
1041       break;
1042     if (changed)
1043       reassigned = 0;
1044 
1045     if (fitType_ == kRefitGhostTrackWithVertices) {
1046       refitGhostTrack(vertices, info);
1047       changed = true;
1048     }
1049 
1050     try {
1051 #ifdef DEBUG
1052       for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter)
1053         debugVertex(*iter, ghostTrack.prediction());
1054       std::cout << "----- reassignment: ---------" << std::endl;
1055 #endif
1056       if (reassignTracks(vertices, info)) {
1057         reassigned++;
1058         changed = true;
1059       } else
1060         reassigned = 0;
1061     } catch (const VertexException &e) {
1062       // just keep vertices as they are
1063     }
1064 
1065     if (!changed)
1066       break;
1067   }
1068 
1069   for (std::vector<CachingVertex<5> >::const_iterator iter = vertices.begin(); iter != vertices.end(); ++iter) {
1070     std::vector<RefCountedVertexTrack> tracks = iter->tracks();
1071     std::vector<RefCountedVertexTrack> newTracks;
1072     newTracks.reserve(tracks.size());
1073 
1074     std::remove_copy_if(tracks.begin(), tracks.end(), std::back_inserter(newTracks), VtxTrackIs(info.ghostTrack));
1075 
1076     if (newTracks.empty())
1077       continue;
1078 
1079     CachingVertex<5> vtx(iter->vertexState(), newTracks, iter->totalChiSquared());
1080     result.push_back(vtx);
1081   }
1082 
1083   return result;
1084 }