AdaptiveVertexFit/src/AdaptiveVertexFitter.cc

0001 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0002 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0003 #include "RecoVertex/VertexTools/interface/AnnealingSchedule.h"
0004 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
0005 #include "RecoVertex/VertexTools/interface/VertexTrackFactory.h"
0006 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include <algorithm>
0010
0011 using namespace edm;
0012
0013 // #define STORE_WEIGHTS
0014 #ifdef STORE_WEIGHTS
0015 #include <dataharvester/Writer.h>
0016 #endif
0017
0018 using namespace std;
0019
0020 namespace {
0021   void sortTracksByPt(std::vector<reco::TransientTrack>& cont) {
0022     auto s = cont.size();
0023     float pt2[s];
0024     int ind[s];
0025     int i = 0;
0026     for (auto const& tk : cont) {
0027       ind[i] = i;
0028       pt2[i++] = tk.impactPointState().globalMomentum().perp2();
0029     }
0030     //clang can not handle lambdas with variable length arrays
0031     auto* p_pt2 = pt2;
0032     std::sort(ind, ind + s, [p_pt2](int i, int j) { return p_pt2[i] > p_pt2[j]; });
0033     std::vector<reco::TransientTrack> tmp;
0034     tmp.reserve(s);
0035     for (auto i = 0U; i < s; ++i)
0036       tmp.emplace_back(std::move(cont[ind[i]]));
0037     cont.swap(tmp);
0038   }
0039
0040   // typedef ReferenceCountingPointer<VertexTrack<5> > RefCountedVertexTrack;
0041   typedef AdaptiveVertexFitter::RefCountedVertexTrack RefCountedVertexTrack;
0042
0043   AlgebraicSymMatrix33 initFitError() {
0044     // that's how we model the lin pt error for the initial seed!
0045     const float initialError = 10000;
0046     AlgebraicSymMatrix33 ret;
0047     ret(0, 0) = initialError;
0048     ret(1, 1) = initialError;
0049     ret(2, 2) = initialError;
0050     return ret;
0051   }
0052
0053   GlobalError const fitError = initFitError();
0054
0055   AlgebraicSymMatrix33 initLinePointError() {
0056     // that's how we model the error of the linearization point.
0057     // for track weighting!
0058     AlgebraicSymMatrix33 ret;
0059     ret(0, 0) = .3;
0060     ret(1, 1) = .3;
0061     ret(2, 2) = 3.;
0062     // ret(0,0)=1e-7; ret(1,1)=1e-7; ret(2,2)=1e-7;
0063     return ret;
0064   }
0065
0066   GlobalError const linPointError = initLinePointError();
0067
0068   void sortByDistanceToRefPoint(std::vector<RefCountedVertexTrack>& cont, const GlobalPoint ref) {
0069     auto s = cont.size();
0070     float d2[s];
0071     int ind[s];
0072     int i = 0;
0073     for (auto const& tk : cont) {
0074       ind[i] = i;
0075       d2[i++] = (tk->linearizedTrack()->track().initialFreeState().position() - ref).mag2();
0076     }
0077     //clang can not handle lambdas with variable length arrays
0078     auto* p_d2 = d2;
0079     std::sort(ind, ind + s, [p_d2](int i, int j) { return p_d2[i] < p_d2[j]; });
0080     std::vector<RefCountedVertexTrack> tmp;
0081     tmp.reserve(s);
0082     for (auto i = 0U; i < s; ++i)
0083       tmp.emplace_back(std::move(cont[ind[i]]));
0084     cont.swap(tmp);
0085   }
0086
0087 #ifdef STORE_WEIGHTS
0088   //NOTE: This is not thread safe
0089   map<RefCountedLinearizedTrackState, int> ids;
0090   int iter = 0;
0091
0092   int getId(const RefCountedLinearizedTrackState& r) {
0093     static int ctr = 1;
0094     if (ids.count(r) == 0) {
0095       ids[r] = ctr++;
0096     }
0097     return ids[r];
0098   }
0099 #endif
0100 }  // namespace
0101
0102 AdaptiveVertexFitter::AdaptiveVertexFitter(const AnnealingSchedule& ann,
0103                                            const LinearizationPointFinder& linP,
0104                                            const VertexUpdator<5>& updator,
0105                                            const VertexTrackCompatibilityEstimator<5>& crit,
0106                                            const VertexSmoother<5>& smoother,
0107                                            const AbstractLTSFactory<5>& ltsf)
0108     : theNr(0),
0109       theLinP(linP.clone()),
0110       theUpdator(updator.clone()),
0111       theSmoother(smoother.clone()),
0112       theAssProbComputer(ann.clone()),
0113       theComp(crit.clone()),
0114       theLinTrkFactory(ltsf.clone()),
0115       gsfIntermediarySmoothing_(false) {
0116   setParameters();
0117 }
0118
0119 void AdaptiveVertexFitter::setWeightThreshold(float w) { theWeightThreshold = w; }
0120
0121 AdaptiveVertexFitter::AdaptiveVertexFitter(const AdaptiveVertexFitter& o)
0122     : theMaxShift(o.theMaxShift),
0123       theMaxLPShift(o.theMaxLPShift),
0124       theMaxStep(o.theMaxStep),
0125       theWeightThreshold(o.theWeightThreshold),
0126       theNr(o.theNr),
0127       theLinP(o.theLinP->clone()),
0128       theUpdator(o.theUpdator->clone()),
0129       theSmoother(o.theSmoother->clone()),
0130       theAssProbComputer(o.theAssProbComputer->clone()),
0131       theComp(o.theComp->clone()),
0132       theLinTrkFactory(o.theLinTrkFactory->clone()),
0133       gsfIntermediarySmoothing_(o.gsfIntermediarySmoothing_) {}
0134
0135 AdaptiveVertexFitter::~AdaptiveVertexFitter() {
0136   delete theLinP;
0137   delete theUpdator;
0138   delete theSmoother;
0139   delete theAssProbComputer;
0140   delete theComp;
0141   delete theLinTrkFactory;
0142 }
0143
0144 void AdaptiveVertexFitter::setParameters(double maxshift, double maxlpshift, unsigned maxstep, double weightthreshold) {
0145   theMaxShift = maxshift;
0146   theMaxLPShift = maxlpshift;
0147   theMaxStep = maxstep;
0148   theWeightThreshold = weightthreshold;
0149 }
0150
0151 void AdaptiveVertexFitter::setParameters(const edm::ParameterSet& s) {
0152   setParameters(s.getParameter<double>("maxshift"),
0153                 s.getParameter<double>("maxlpshift"),
0154                 s.getParameter<int>("maxstep"),
0155                 s.getParameter<double>("weightthreshold"));
0156 }
0157
0158 CachingVertex<5> AdaptiveVertexFitter::vertex(const vector<reco::TransientTrack>& unstracks) const {
0159   if (unstracks.size() < 2) {
0160     LogError("RecoVertex|AdaptiveVertexFitter") << "Supplied fewer than two tracks. Vertex is invalid.";
0161     return CachingVertex<5>();  // return invalid vertex
0162   };
0163   vector<reco::TransientTrack> tracks = unstracks;
0164   sortTracksByPt(tracks);
0165   // Linearization Point
0166   GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
0167   // Initial vertex seed, with a very large error matrix
0168   VertexState lseed(linP, linPointError);
0169   vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, lseed);
0170
0171   VertexState seed(linP, fitError);
0172   return fit(vtContainer, seed, false);
0173 }
0174
0175 CachingVertex<5> AdaptiveVertexFitter::vertex(const vector<RefCountedVertexTrack>& tracks) const {
0176   if (tracks.size() < 2) {
0177     LogError("RecoVertex|AdaptiveVertexFitter") << "Supplied fewer than two tracks. Vertex is invalid.";
0178     return CachingVertex<5>();  // return invalid vertex
0179   };
0180   // Initial vertex seed, with a very small weight matrix
0181   GlobalPoint linP = tracks[0]->linearizedTrack()->linearizationPoint();
0182   VertexState seed(linP, fitError);
0183   return fit(tracks, seed, false);
0184 }
0185
0186 CachingVertex<5> AdaptiveVertexFitter::vertex(const vector<RefCountedVertexTrack>& tracks,
0187                                               const reco::BeamSpot& spot) const {
0188   if (tracks.empty()) {
0189     LogError("RecoVertex|AdaptiveVertexFitter") << "Supplied no tracks. Vertex is invalid.";
0190     return CachingVertex<5>();  // return invalid vertex
0191   };
0192   VertexState beamSpotState(spot);
0193   return fit(tracks, beamSpotState, true);
0194 }
0195
0196 /** Fit vertex out of a set of reco::TransientTracks.
0197  *  Uses the specified linearization point.
0198  */
0199 CachingVertex<5> AdaptiveVertexFitter::vertex(const vector<reco::TransientTrack>& tracks,
0200                                               const GlobalPoint& linPoint) const {
0201   if (tracks.size() < 2) {
0202     LogError("RecoVertex|AdaptiveVertexFitter") << "Supplied fewer than two tracks. Vertex is invalid.";
0203     return CachingVertex<5>();  // return invalid vertex
0204   };
0205   // Initial vertex seed, with a very large error matrix
0206   VertexState seed(linPoint, linPointError);
0207   vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, seed);
0208   VertexState fitseed(linPoint, fitError);
0209   return fit(vtContainer, fitseed, false);
0210 }
0211
0212 /** Fit vertex out of a set of TransientTracks.
0213  *  The specified BeamSpot will be used as priot, but NOT for the linearization.
0214  *  The specified LinearizationPointFinder will be used to find the linearization point.
0215  */
0216 CachingVertex<5> AdaptiveVertexFitter::vertex(const vector<reco::TransientTrack>& unstracks,
0217                                               const reco::BeamSpot& beamSpot) const {
0218   if (unstracks.empty()) {
0219     LogError("RecoVertex|AdaptiveVertexFitter") << "Supplied no tracks. Vertex is invalid.";
0220     return CachingVertex<5>();  // return invalid vertex
0221   };
0222
0223   VertexState beamSpotState(beamSpot);
0224   vector<RefCountedVertexTrack> vtContainer;
0225
0226   vector<reco::TransientTrack> tracks = unstracks;
0227   sortTracksByPt(tracks);
0228
0229   if (tracks.size() > 1) {
0230     // Linearization Point search if there are more than 1 track
0231     GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
0232     VertexState lpState(linP, linPointError);
0233     vtContainer = linearizeTracks(tracks, lpState);
0234   } else {
0235     // otherwise take the beamspot position.
0236     vtContainer = linearizeTracks(tracks, beamSpotState);
0237   }
0238
0239   return fit(vtContainer, beamSpotState, true);
0240 }
0241
0242 /** Fit vertex out of a set of reco::TransientTracks.
0243  *   Uses the position as both the linearization point AND as prior
0244  *   estimate of the vertex position. The error is used for the
0245  *   weight of the prior estimate.
0246  */
0247 CachingVertex<5> AdaptiveVertexFitter::vertex(const vector<reco::TransientTrack>& tracks,
0248                                               const GlobalPoint& priorPos,
0249                                               const GlobalError& priorError) const
0250
0251 {
0252   if (tracks.empty()) {
0253     LogError("RecoVertex|AdaptiveVertexFitter") << "Supplied no tracks. Vertex is invalid.";
0254     return CachingVertex<5>();  // return invalid vertex
0255   };
0256   VertexState seed(priorPos, priorError);
0257   vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, seed);
0258   return fit(vtContainer, seed, true);
0259 }
0260
0261 /** Fit vertex out of a set of VertexTracks
0262  *   Uses the position and error for the prior estimate of the vertex.
0263  *   This position is not used to relinearize the tracks.
0264  */
0265 CachingVertex<5> AdaptiveVertexFitter::vertex(const vector<RefCountedVertexTrack>& tracks,
0266                                               const GlobalPoint& priorPos,
0267                                               const GlobalError& priorError) const {
0268   if (tracks.empty()) {
0269     LogError("RecoVertex|AdaptiveVertexFitter") << "Supplied no tracks. Vertex is invalid.";
0270     return CachingVertex<5>();  // return invalid vertex
0271   };
0272   VertexState seed(priorPos, priorError);
0273   return fit(tracks, seed, true);
0274 }
0275
0276 /**
0277  * Construct a container of VertexTrack from a set of reco::TransientTracks.
0278  * As this is the first iteration of the adaptive fit, the initial error
0279  * does not enter in the computation of the weights.
0280  * This is to avoid that all tracks get the same weight when
0281  * using a very large initial error matrix.
0282  */
0283 vector<AdaptiveVertexFitter::RefCountedVertexTrack> AdaptiveVertexFitter::linearizeTracks(
0284     const vector<reco::TransientTrack>& tracks, const VertexState& seed) const {
0285   const GlobalPoint& linP(seed.position());
0286   vector<RefCountedLinearizedTrackState> lTracks;
0287   for (vector<reco::TransientTrack>::const_iterator i = tracks.begin(); i != tracks.end(); ++i) {
0288     try {
0289       RefCountedLinearizedTrackState lTrData = theLinTrkFactory->linearizedTrackState(linP, *i);
0290       lTracks.push_back(lTrData);
0291     } catch (exception& e) {
0292       LogInfo("RecoVertex/AdaptiveVertexFitter") << "Exception " << e.what() << " in ::linearizeTracks."
0293                                                  << "Your future vertex has just lost a track.";
0294     };
0295   }
0296   return weightTracks(lTracks, seed);
0297 }
0298
0299 /**
0300  * Construct new a container of VertexTrack with a new linearization point
0301  * and vertex seed, from an existing set of VertexTrack, from which only the
0302  * recTracks will be used.
0303  */
0304 vector<AdaptiveVertexFitter::RefCountedVertexTrack> AdaptiveVertexFitter::reLinearizeTracks(
0305     const vector<RefCountedVertexTrack>& tracks, const CachingVertex<5>& vertex) const {
0306   const VertexState& seed = vertex.vertexState();
0307   GlobalPoint linP = seed.position();
0308   vector<RefCountedLinearizedTrackState> lTracks;
0309   for (vector<RefCountedVertexTrack>::const_iterator i = tracks.begin(); i != tracks.end(); i++) {
0310     try {
0311       RefCountedLinearizedTrackState lTrData =
0312           theLinTrkFactory->linearizedTrackState(linP, (**i).linearizedTrack()->track());
0313       /*
0314       RefCountedLinearizedTrackState lTrData =
0315               (**i).linearizedTrack()->stateWithNewLinearizationPoint(linP);
0316               */
0317       lTracks.push_back(lTrData);
0318     } catch (exception& e) {
0319       LogInfo("RecoVertex/AdaptiveVertexFitter") << "Exception " << e.what() << " in ::relinearizeTracks. "
0320                                                  << "Will not relinearize this track.";
0321       lTracks.push_back((**i).linearizedTrack());
0322     };
0323   };
0324   return reWeightTracks(lTracks, vertex);
0325 }
0326
0327 AdaptiveVertexFitter* AdaptiveVertexFitter::clone() const { return new AdaptiveVertexFitter(*this); }
0328
0329 double AdaptiveVertexFitter::getWeight(float chi2) const {
0330   double weight = theAssProbComputer->weight(chi2);
0331
0332   if (weight > 1.0) {
0333     LogInfo("RecoVertex/AdaptiveVertexFitter") << "Weight " << weight << " > 1.0!";
0334     weight = 1.0;
0335   };
0336
0337   if (weight < 1e-20) {
0338     // LogInfo("RecoVertex/AdaptiveVertexFitter") << "Weight " << weight << " < 0.0!";
0339     weight = 1e-20;
0340   };
0341   return weight;
0342 }
0343
0344 vector<AdaptiveVertexFitter::RefCountedVertexTrack> AdaptiveVertexFitter::reWeightTracks(
0345     const vector<RefCountedLinearizedTrackState>& lTracks, const CachingVertex<5>& vertex) const {
0346   const VertexState& seed = vertex.vertexState();
0347   // cout << "[AdaptiveVertexFitter] now reweight around " << seed.position() << endl;
0348   theNr++;
0349   // GlobalPoint pos = seed.position();
0350
0351   vector<RefCountedVertexTrack> finalTracks;
0352   VertexTrackFactory<5> vTrackFactory;
0353 #ifdef STORE_WEIGHTS
0354   iter++;
0355 #endif
0356   for (vector<RefCountedLinearizedTrackState>::const_iterator i = lTracks.begin(); i != lTracks.end(); i++) {
0357     double weight = 0.;
0358     // cout << "[AdaptiveVertexFitter] estimate " << endl;
0359     pair<bool, double> chi2Res(false, 0.);
0360     try {
0361       chi2Res = theComp->estimate(vertex, *i, std::distance(lTracks.begin(), i));
0362     } catch (exception const& e) {
0363     };
0364     // cout << "[AdaptiveVertexFitter] /estimate " << endl;
0365     if (!chi2Res.first) {
0366       // cout << "[AdaptiveVertexFitter] aie... vertex candidate is at  " << vertex.position() << endl;
0367       LogInfo("AdaptiveVertexFitter") << "When reweighting, chi2<0. Will add this track with w=0.";
0368       // edm::LogInfo("AdaptiveVertexFitter" ) << "pt=" << (**i).track().pt();
0369     } else {
0370       weight = getWeight(chi2Res.second);
0371     }
0372
0373     RefCountedVertexTrack vTrData = vTrackFactory.vertexTrack(*i, seed, weight);
0374
0375 #ifdef STORE_WEIGHTS
0376     map<string, dataharvester::MultiType> m;
0377     m["chi2"] = chi2;
0378     m["w"] = theAssProbComputer->weight(chi2);
0379     m["T"] = theAssProbComputer->currentTemp();
0380     m["n"] = iter;
0381     m["pos"] = "reweight";
0382     m["id"] = getId(*i);
0383     dataharvester::Writer::file("w.txt").save(m);
0384 #endif
0385
0386     finalTracks.push_back(vTrData);
0387   }
0388   sortByDistanceToRefPoint(finalTracks, vertex.position());
0389   // cout << "[AdaptiveVertexFitter] /now reweight" << endl;
0390   return finalTracks;
0391 }
0392
0393 vector<AdaptiveVertexFitter::RefCountedVertexTrack> AdaptiveVertexFitter::weightTracks(
0394     const vector<RefCountedLinearizedTrackState>& lTracks, const VertexState& seed) const {
0395   theNr++;
0396   CachingVertex<5> seedvtx(seed, vector<RefCountedVertexTrack>(), 0.);
0397   /** track weighting, as opposed to re-weighting, must always
0398    * be done with a reset annealer! */
0399   theAssProbComputer->resetAnnealing();
0400
0401   vector<RefCountedVertexTrack> finalTracks;
0402   VertexTrackFactory<5> vTrackFactory;
0403 #ifdef STORE_WEIGHTS
0404   iter++;
0405 #endif
0406   for (vector<RefCountedLinearizedTrackState>::const_iterator i = lTracks.begin(); i != lTracks.end(); i++) {
0407     double weight = 0.;
0408     pair<bool, double> chi2Res = theComp->estimate(seedvtx, *i, std::distance(lTracks.begin(), i));
0409     if (!chi2Res.first) {
0410       // cout << "[AdaptiveVertexFitter] Aiee! " << endl;
0411       LogInfo("AdaptiveVertexFitter") << "When weighting a track, chi2 calculation failed;"
0412                                       << " will add with w=0.";
0413     } else {
0414       weight = getWeight(chi2Res.second);
0415     }
0416     RefCountedVertexTrack vTrData = vTrackFactory.vertexTrack(*i, seed, weight);
0417 #ifdef STORE_WEIGHTS
0418     map<string, dataharvester::MultiType> m;
0419     m["chi2"] = chi2;
0420     m["w"] = theAssProbComputer->weight(chi2);
0421     m["T"] = theAssProbComputer->currentTemp();
0422     m["n"] = iter;
0423     m["id"] = getId(*i);
0424     m["pos"] = "weight";
0425     dataharvester::Writer::file("w.txt").save(m);
0426 #endif
0427     finalTracks.push_back(vTrData);
0428   }
0429   return finalTracks;
0430 }
0431
0432 /**
0433  * Construct new a container of VertexTrack with new weights
0434  * accounting for vertex error, from an existing set of VertexTracks.
0435  * From these the LinearizedTracks will be reused.
0436  */
0437 vector<AdaptiveVertexFitter::RefCountedVertexTrack> AdaptiveVertexFitter::reWeightTracks(
0438     const vector<RefCountedVertexTrack>& tracks, const CachingVertex<5>& seed) const {
0439   vector<RefCountedLinearizedTrackState> lTracks;
0440   for (vector<RefCountedVertexTrack>::const_iterator i = tracks.begin(); i != tracks.end(); i++) {
0441     lTracks.push_back((**i).linearizedTrack());
0442   }
0443
0444   return reWeightTracks(lTracks, seed);
0445 }
0446
0447 /*
0448  * The method where the vertex fit is actually done!
0449  */
0450
0451 CachingVertex<5> AdaptiveVertexFitter::fit(const vector<RefCountedVertexTrack>& tracks,
0452                                            const VertexState& priorSeed,
0453                                            bool withPrior) const {
0454   // cout << "[AdaptiveVertexFit] fit with " << tracks.size() << endl;
0455   theAssProbComputer->resetAnnealing();
0456
0457   vector<RefCountedVertexTrack> initialTracks;
0458   GlobalPoint priorVertexPosition = priorSeed.position();
0459   GlobalError priorVertexError = priorSeed.error();
0460
0461   CachingVertex<5> returnVertex(priorVertexPosition, priorVertexError, initialTracks, 0);
0462   if (withPrior) {
0463     returnVertex = CachingVertex<5>(
0464         priorVertexPosition, priorVertexError, priorVertexPosition, priorVertexError, initialTracks, 0);
0465   }
0466
0467   std::vector<RefCountedVertexTrack> globalVTracks = tracks;
0468   // sort the tracks, according to distance to seed!
0469   sortByDistanceToRefPoint(globalVTracks, priorSeed.position());
0470
0471   // main loop through all the VTracks
0472   int step = 0;
0473
0474   CachingVertex<5> initialVertex = returnVertex;
0475
0476   GlobalPoint newPosition = priorVertexPosition;
0477   GlobalPoint previousPosition = newPosition;
0478
0479   int ns_trks = 0;  // number of significant tracks.
0480   // If we have only two significant tracks, we return an invalid vertex
0481
0482   // cout << "[AdaptiveVertexFit] start " << tracks.size() << endl;
0483   /*
0484   for ( vector< RefCountedVertexTrack >::const_iterator
0485         i=globalVTracks.begin(); i!=globalVTracks.end() ; ++i )
0486   {
0487     cout << "  " << (**i).linearizedTrack()->track().initialFreeState().momentum() << endl;
0488   }*/
0489   do {
0490     ns_trks = 0;
0491     CachingVertex<5> fVertex = initialVertex;
0492     // cout << "[AdaptiveVertexFit] step " << step << " at " << fVertex.position() << endl;
0493     if ((previousPosition - newPosition).transverse() > theMaxLPShift) {
0494       // relinearize and reweight.
0495       // (reLinearizeTracks also reweights tracks)
0496       // cout << "[AdaptiveVertexFit] relinearize at " << returnVertex.position() << endl;
0497       if (gsfIntermediarySmoothing_)
0498         returnVertex = theSmoother->smooth(returnVertex);
0499       globalVTracks = reLinearizeTracks(globalVTracks, returnVertex);
0500     } else if (step) {
0501       // reweight, if it is not the first step
0502       // cout << "[AdaptiveVertexFit] reweight at " << returnVertex.position() << endl;
0503       if (gsfIntermediarySmoothing_)
0504         returnVertex = theSmoother->smooth(returnVertex);
0505       globalVTracks = reWeightTracks(globalVTracks, returnVertex);
0506     }
0507     // cout << "[AdaptiveVertexFit] relinarized, reweighted" << endl;
0508     // update sequentially the vertex estimate
0509     CachingVertex<5> nVertex;
0510     for (vector<RefCountedVertexTrack>::const_iterator i = globalVTracks.begin(); i != globalVTracks.end(); i++) {
0511       if ((**i).weight() > 0.)
0512         nVertex = theUpdator->add(fVertex, *i);
0513       else
0514         nVertex = fVertex;
0515       if (nVertex.isValid()) {
0516         if ((**i).weight() >= theWeightThreshold)
0517           ns_trks++;
0518
0519         if (fabs(nVertex.position().z()) > 10000. || nVertex.position().perp() > 120.) {
0520           // were more than 100 m off!!
0521           LogInfo("AdaptiveVertexFitter")
0522               << "Vertex candidate just took off to " << nVertex.position() << "! Will discard this update!";
0523           //        //<< "track pt was " << (**i).linearizedTrack()->track().pt()
0524           //                         << "track momentum was " << (**i).linearizedTrack()->track().initialFreeState().momentum()
0525           //                         << "track position was " << (**i).linearizedTrack()->track().initialFreeState().position()
0526           //                         << "track chi2 was " << (**i).linearizedTrack()->track().chi2()
0527           //                         << "track ndof was " << (**i).linearizedTrack()->track().ndof()
0528           //                         << "track w was " << (**i).weight()
0529           //                         << "track schi2 was " << (**i).smoothedChi2();
0530         } else {
0531           fVertex = nVertex;
0532         }
0533       } else {
0534         LogInfo("RecoVertex/AdaptiveVertexFitter")
0535             << "The updator returned an invalid vertex when adding track " << i - globalVTracks.begin()
0536             << ".\n Your vertex might just have lost one good track.";
0537       }
0538     }
0539     previousPosition = newPosition;
0540     newPosition = fVertex.position();
0541     returnVertex = fVertex;
0542     theAssProbComputer->anneal();
0543     step++;
0544     if (step >= theMaxStep)
0545       break;
0546
0547   } while (
0548       // repeat as long as
0549       // - vertex moved too much or
0550       // - we're not yet annealed
0551       (((previousPosition - newPosition).mag() > theMaxShift) || (!(theAssProbComputer->isAnnealed()))));
0552
0553   if (theWeightThreshold > 0. && ns_trks < 2 && !withPrior) {
0554     LogDebug("AdaptiveVertexFitter") << "fewer than two significant tracks (w>" << theWeightThreshold << ")."
0555                                      << " Fitted vertex is invalid.";
0556     return CachingVertex<5>();  // return invalid vertex
0557   }
0558
0559 #ifdef STORE_WEIGHTS
0560   map<string, dataharvester::MultiType> m;
0561   m["chi2"] = chi2;
0562   m["w"] = theAssProbComputer->weight(chi2);
0563   m["T"] = theAssProbComputer->currentTemp();
0564   m["n"] = iter;
0565   m["id"] = getId(*i);
0566   m["pos"] = "final";
0567   dataharvester::Writer::file("w.txt").save(m);
0568 #endif
0569   // cout << "[AdaptiveVertexFit] /fit" << endl;
0570   return theSmoother->smooth(returnVertex);
0571 }