Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:28:03

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   int ctr = 0;
0175   unsigned int n_tracks = remainingtrks.size();
0176 
0177   // cout << "[AdaptiveVertexReconstructor] DEBUG ::vertices!!" << endl;
0178   try {
0179     while (remainingtrks.size() > 1) {
0180       /*
0181       cout << "[AdaptiveVertexReconstructor] next round: "
0182            << remainingtrks.size() << endl;
0183            */
0184       ctr++;
0185       const AdaptiveVertexFitter* fitter = theSecondaryFitter;
0186       if (ret.empty()) {
0187         fitter = thePrimaryFitter;
0188       };
0189       vector<reco::TransientTrack> fittrks;
0190       fittrks.reserve(remainingtrks.size());
0191 
0192       copy(remainingtrks.begin(), remainingtrks.end(), back_inserter(fittrks));
0193 
0194       TransientVertex tmpvtx;
0195       if ((ret.empty()) && has_primaries) {
0196         // add the primaries to the fitted tracks.
0197         copy(primaries.begin(), primaries.end(), back_inserter(fittrks));
0198       }
0199       if ((ret.empty()) && usespot) {
0200         tmpvtx = fitter->vertex(fittrks, s);
0201       } else {
0202         tmpvtx = fitter->vertex(fittrks);
0203       }
0204       TransientVertex newvtx = cleanUp(tmpvtx);
0205       ret.push_back(newvtx);
0206       erase(newvtx, remainingtrks, theMinWeight);
0207       if (n_tracks == remainingtrks.size()) {
0208         if (usespot) {
0209           // try once more without beamspot constraint!
0210           usespot = false;
0211           LogDebug("AdaptiveVertexReconstructor") << "no tracks in vertex. trying again without beamspot constraint!";
0212           continue;
0213         }
0214         LogDebug("AdaptiveVertexReconstructor")
0215             << "all tracks (" << n_tracks << ") would be recycled for next fit. Trying with low threshold!";
0216         erase(newvtx, remainingtrks, 1.e-5);
0217         if (n_tracks == remainingtrks.size()) {
0218           LogDebug("AdaptiveVertexReconstructor") << "low threshold didnt help! "
0219                                                   << "Discontinue procedure!";
0220           break;
0221         }
0222       };
0223 
0224       // cout << "[AdaptiveVertexReconstructor] erased" << endl;
0225       n_tracks = remainingtrks.size();
0226     };
0227   } catch (VertexException& v) {
0228     // Will catch all (not enough significant tracks exceptions.
0229     // in this case, the iteration can safely terminate.
0230   };
0231 
0232   return cleanUpVertices(ret);
0233 }
0234 
0235 vector<TransientVertex> AdaptiveVertexReconstructor::cleanUpVertices(const vector<TransientVertex>& old) const {
0236   vector<TransientVertex> ret;
0237   for (vector<TransientVertex>::const_iterator i = old.begin(); i != old.end(); ++i) {
0238     if (!(i->hasTrackWeight())) {  // if we dont have track weights, we take the vtx
0239       ret.push_back(*i);
0240       continue;
0241     }
0242 
0243     // maybe this should be replaced with asking for the ndf ...
0244     // e.g. if ( ndf > - 1. )
0245     int nsig = 0;  // number of significant tracks.
0246     TransientVertex::TransientTrackToFloatMap wm = i->weightMap();
0247     for (TransientVertex::TransientTrackToFloatMap::const_iterator w = wm.begin(); w != wm.end(); ++w) {
0248       if (w->second > theWeightThreshold)
0249         nsig++;
0250     }
0251     if (nsig > 1)
0252       ret.push_back(*i);
0253   }
0254 
0255   return ret;
0256 }