GhostTrackFitter/src/SequentialGhostTrackFitter.cc

0001 #include <vector>
0002
0003 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
0004 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
0005
0006 #include "RecoVertex/GhostTrackFitter/interface/SequentialGhostTrackFitter.h"
0007
0008 using namespace reco;
0009
0010 namespace {
0011   inline double sqr(double arg) { return arg * arg; }
0012 }  // namespace
0013
0014 SequentialGhostTrackFitter::SequentialGhostTrackFitter()
0015     : maxIteration(15), minDeltaR(0.0015), minDistance(0.002), weightThreshold(0.001) {}
0016
0017 bool SequentialGhostTrackFitter::stable(const GhostTrackPrediction &before, const GhostTrackPrediction &after) const {
0018   return (sqr(after.sz() - before.sz()) + sqr(after.ip() - before.ip()) < sqr(minDistance) &&
0019           sqr(after.eta() - before.eta()) + sqr(after.phi() - before.phi()) < sqr(minDeltaR));
0020 }
0021
0022 GhostTrackPrediction SequentialGhostTrackFitter::fit(const GhostTrackFitter::PredictionUpdater &updater,
0023                                                      const GhostTrackPrediction &prior,
0024                                                      std::vector<GhostTrackState> &states,
0025                                                      double &ndof,
0026                                                      double &chi2) {
0027   GhostTrackPrediction pred, lastPred = prior;
0028
0029   reset();
0030
0031   ndof = 0.;
0032   chi2 = 0.;
0033
0034   unsigned int iteration = 0;
0035   for (;;) {
0036     pred = prior;
0037
0038     if (states.begin() == states.end())
0039       break;
0040
0041     if (iteration > 0) {
0042       for (unsigned int i = 0; i < states.size(); i++) {
0043         GhostTrackState &state = states[i];
0044         state.linearize(lastPred);
0045       }
0046     }
0047
0048     ndof = 0.;  // prior gives us an initial ndof
0049     chi2 = 0.;
0050
0051     for (std::vector<GhostTrackState>::const_iterator state = states.begin(); state != states.end(); ++state) {
0052       if (state->isValid() && state->weight() > weightThreshold)
0053         pred = updater.update(pred, *state, ndof, chi2);
0054     }
0055
0056     if (++iteration >= maxIteration || stable(lastPred, pred))
0057       break;
0058
0059     postFit(updater, pred, states);
0060
0061     lastPred = pred;
0062   }
0063
0064   return pred;
0065 }