AdaptiveVertexFinder/src/AdaptiveVertexReconstructor.cc

0001 #include "RecoVertex/AdaptiveVertexFinder/interface/AdaptiveVertexReconstructor.h"
0002 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
0003 #include "FWCore/Utilities/interface/EDMException.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "RecoVertex/VertexTools/interface/DummyVertexSmoother.h"
0006 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
0007 #include <algorithm>
0008
0009 using namespace std;
0010
0011 TransientVertex AdaptiveVertexReconstructor::cleanUp(const TransientVertex& old) const {
0012   vector<reco::TransientTrack> trks = old.originalTracks();
0013   vector<reco::TransientTrack> newtrks;
0014   TransientVertex::TransientTrackToFloatMap mp;
0015   static const float minweight = 1.e-8;  // discard all tracks with lower weight
0016   for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0017     if (old.trackWeight(*i) > minweight) {
0018       newtrks.push_back(*i);
0019       mp[*i] = old.trackWeight(*i);
0020     }
0021   }
0022
0023   TransientVertex ret;
0024
0025   if (old.hasPrior()) {
0026     VertexState priorstate(old.priorPosition(), old.priorError());
0027     ret = TransientVertex(priorstate, old.vertexState(), newtrks, old.totalChiSquared(), old.degreesOfFreedom());
0028   } else {
0029     ret = TransientVertex(old.vertexState(), newtrks, old.totalChiSquared(), old.degreesOfFreedom());
0030   }
0031   ret.weightMap(mp);  // set weight map
0032
0033   if (old.hasRefittedTracks()) {
0034     // we have refitted tracks -- copy them!
0035     vector<reco::TransientTrack> newrfs;
0036     vector<reco::TransientTrack> oldrfs = old.refittedTracks();
0037     vector<reco::TransientTrack>::const_iterator origtrkiter = trks.begin();
0038     for (vector<reco::TransientTrack>::const_iterator i = oldrfs.begin(); i != oldrfs.end(); ++i) {
0039       if (old.trackWeight(*origtrkiter) > minweight) {
0040         newrfs.push_back(*i);
0041       }
0042       origtrkiter++;
0043     }
0044     if (!newrfs.empty())
0045       ret.refittedTracks(newrfs);  // copy refitted tracks
0046   }
0047
0048   if (ret.refittedTracks().size() > ret.originalTracks().size()) {
0049     edm::LogError("AdaptiveVertexReconstructor") << "More refitted tracks than original tracks!";
0050   }
0051
0052   return ret;
0053 }
0054
0055 void AdaptiveVertexReconstructor::erase(const TransientVertex& newvtx,
0056                                         set<reco::TransientTrack>& remainingtrks,
0057                                         float w) const {
0058   /*
0059    * Erase tracks that are in newvtx from remainingtrks
0060    * But erase only if trackweight > w
0061    */
0062   const vector<reco::TransientTrack>& origtrks = newvtx.originalTracks();
0063   for (vector<reco::TransientTrack>::const_iterator i = origtrks.begin(); i != origtrks.end(); ++i) {
0064     double weight = newvtx.trackWeight(*i);
0065     if (weight > w) {
0066       remainingtrks.erase(*i);
0067     };
0068   };
0069 }
0070
0071 AdaptiveVertexReconstructor::AdaptiveVertexReconstructor(float primcut, float seccut, float min_weight, bool smoothing)
0072     : thePrimaryFitter(nullptr), theSecondaryFitter(nullptr), theMinWeight(min_weight), theWeightThreshold(0.001) {
0073   setupFitters(primcut, 256., 0.25, seccut, 256., 0.25, smoothing);
0074 }
0075
0076 AdaptiveVertexReconstructor::~AdaptiveVertexReconstructor() {
0077   if (thePrimaryFitter)
0078     delete thePrimaryFitter;
0079   if (theSecondaryFitter)
0080     delete theSecondaryFitter;
0081 }
0082
0083 void AdaptiveVertexReconstructor::setupFitters(
0084     float primcut, float primT, float primr, float seccut, float secT, float secr, bool smoothing) {
0085   VertexSmoother<5>* smoother;
0086   if (smoothing) {
0087     smoother = new KalmanVertexSmoother();
0088   } else {
0089     smoother = new DummyVertexSmoother<5>();
0090   }
0091
0092   if (thePrimaryFitter)
0093     delete thePrimaryFitter;
0094   if (theSecondaryFitter)
0095     delete theSecondaryFitter;
0096
0097   /*
0098   edm::LogError ("AdaptiveVertexReconstructor" )
0099     << "Tini and r are hardcoded now!";
0100     */
0101   thePrimaryFitter = new AdaptiveVertexFitter(GeometricAnnealing(primcut, primT, primr),
0102                                               DefaultLinearizationPointFinder(),
0103                                               KalmanVertexUpdator<5>(),
0104                                               KalmanVertexTrackCompatibilityEstimator<5>(),
0105                                               *smoother);
0106   thePrimaryFitter->setWeightThreshold(theWeightThreshold);
0107   // if the primary fails, sth is wrong, so here we set a threshold on the weight.
0108   theSecondaryFitter = new AdaptiveVertexFitter(GeometricAnnealing(seccut, secT, secr),
0109                                                 DefaultLinearizationPointFinder(),
0110                                                 KalmanVertexUpdator<5>(),
0111                                                 KalmanVertexTrackCompatibilityEstimator<5>(),
0112                                                 *smoother);
0113   theSecondaryFitter->setWeightThreshold(0.);
0114   // need to set it or else we have
0115   // unwanted exceptions to deal with.
0116   // cleanup can come later!
0117   delete smoother;
0118 }
0119
0120 AdaptiveVertexReconstructor::AdaptiveVertexReconstructor(const edm::ParameterSet& m)
0121     : thePrimaryFitter(nullptr), theSecondaryFitter(nullptr), theMinWeight(0.5), theWeightThreshold(0.001) {
0122   float primcut = 2.0;
0123   float seccut = 6.0;
0124   bool smoothing = false;
0125   // float primT = 4096.;
0126   // float primr = 0.125;
0127   float primT = 256.;
0128   float primr = 0.25;
0129   float secT = 256.;
0130   float secr = 0.25;
0131
0132   try {
0133     primcut = m.getParameter<double>("primcut");
0134     primT = m.getParameter<double>("primT");
0135     primr = m.getParameter<double>("primr");
0136     seccut = m.getParameter<double>("seccut");
0137     secT = m.getParameter<double>("secT");
0138     secr = m.getParameter<double>("secr");
0139     theMinWeight = m.getParameter<double>("minweight");
0140     theWeightThreshold = m.getParameter<double>("weightthreshold");
0141     smoothing = m.getParameter<bool>("smoothing");
0142   } catch (edm::Exception& e) {
0143     edm::LogError("AdaptiveVertexReconstructor") << e.what();
0144   }
0145
0146   setupFitters(primcut, primT, primr, seccut, secT, secr, smoothing);
0147 }
0148
0149 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& t,
0150                                                               const reco::BeamSpot& s) const {
0151   return vertices(vector<reco::TransientTrack>(), t, s, false, true);
0152 }
0153
0154 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& primaries,
0155                                                               const vector<reco::TransientTrack>& tracks,
0156                                                               const reco::BeamSpot& s) const {
0157   return vertices(primaries, tracks, s, true, true);
0158 }
0159
0160 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& tracks) const {
0161   return vertices(vector<reco::TransientTrack>(), tracks, reco::BeamSpot(), false, false);
0162 }
0163
0164 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& primaries,
0165                                                               const vector<reco::TransientTrack>& tracks,
0166                                                               const reco::BeamSpot& s,
0167                                                               bool has_primaries,
0168                                                               bool usespot) const {
0169   vector<TransientVertex> ret;
0170   set<reco::TransientTrack> remainingtrks;
0171
0172   copy(tracks.begin(), tracks.end(), inserter(remainingtrks, remainingtrks.begin()));
0173
0174   unsigned int n_tracks = remainingtrks.size();
0175
0176   // cout << "[AdaptiveVertexReconstructor] DEBUG ::vertices!!" << endl;
0177   try {
0178     while (remainingtrks.size() > 1) {
0179       /*
0180       cout << "[AdaptiveVertexReconstructor] next round: "
0181            << remainingtrks.size() << endl;
0182            */
0183       const AdaptiveVertexFitter* fitter = theSecondaryFitter;
0184       if (ret.empty()) {
0185         fitter = thePrimaryFitter;
0186       };
0187       vector<reco::TransientTrack> fittrks;
0188       fittrks.reserve(remainingtrks.size());
0189
0190       copy(remainingtrks.begin(), remainingtrks.end(), back_inserter(fittrks));
0191
0192       TransientVertex tmpvtx;
0193       if ((ret.empty()) && has_primaries) {
0194         // add the primaries to the fitted tracks.
0195         copy(primaries.begin(), primaries.end(), back_inserter(fittrks));
0196       }
0197       if ((ret.empty()) && usespot) {
0198         tmpvtx = fitter->vertex(fittrks, s);
0199       } else {
0200         tmpvtx = fitter->vertex(fittrks);
0201       }
0202       TransientVertex newvtx = cleanUp(tmpvtx);
0203       ret.push_back(newvtx);
0204       erase(newvtx, remainingtrks, theMinWeight);
0205       if (n_tracks == remainingtrks.size()) {
0206         if (usespot) {
0207           // try once more without beamspot constraint!
0208           usespot = false;
0209           LogDebug("AdaptiveVertexReconstructor") << "no tracks in vertex. trying again without beamspot constraint!";
0210           continue;
0211         }
0212         LogDebug("AdaptiveVertexReconstructor")
0213             << "all tracks (" << n_tracks << ") would be recycled for next fit. Trying with low threshold!";
0214         erase(newvtx, remainingtrks, 1.e-5);
0215         if (n_tracks == remainingtrks.size()) {
0216           LogDebug("AdaptiveVertexReconstructor") << "low threshold didnt help! "
0217                                                   << "Discontinue procedure!";
0218           break;
0219         }
0220       };
0221
0222       // cout << "[AdaptiveVertexReconstructor] erased" << endl;
0223       n_tracks = remainingtrks.size();
0224     };
0225   } catch (VertexException& v) {
0226     // Will catch all (not enough significant tracks exceptions.
0227     // in this case, the iteration can safely terminate.
0228   };
0229
0230   return cleanUpVertices(ret);
0231 }
0232
0233 vector<TransientVertex> AdaptiveVertexReconstructor::cleanUpVertices(const vector<TransientVertex>& old) const {
0234   vector<TransientVertex> ret;
0235   for (vector<TransientVertex>::const_iterator i = old.begin(); i != old.end(); ++i) {
0236     if (!(i->hasTrackWeight())) {  // if we dont have track weights, we take the vtx
0237       ret.push_back(*i);
0238       continue;
0239     }
0240
0241     // maybe this should be replaced with asking for the ndf ...
0242     // e.g. if ( ndf > - 1. )
0243     int nsig = 0;  // number of significant tracks.
0244     TransientVertex::TransientTrackToFloatMap wm = i->weightMap();
0245     for (TransientVertex::TransientTrackToFloatMap::const_iterator w = wm.begin(); w != wm.end(); ++w) {
0246       if (w->second > theWeightThreshold)
0247         nsig++;
0248     }
0249     if (nsig > 1)
0250       ret.push_back(*i);
0251   }
0252
0253   return ret;
0254 }