Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-06-24 02:11:19 

0001 #ifndef TrackVertexArbitration_H
0002 #define TrackVertexArbitration_H
0003 #include <memory>
0004 #include <set>
0005 #include <unordered_map>
0006 
0007 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0008 #include "TrackingTools/TransientTrack/interface/CandidatePtrTransientTrack.h"
0009 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0010 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0011 #include "TrackingTools/IPTools/interface/IPTools.h"
0012 
0013 #include "FWCore/Framework/interface/EDProducer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0024 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0025 
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "DataFormats/VertexReco/interface/Vertex.h"
0029 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0030 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0031 
0032 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0033 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h"
0034 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackCompatibilityEstimator.h"
0035 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
0036 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0037 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0038 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0039 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0040 
0041 #include "RecoVertex/AdaptiveVertexFinder/interface/SVTimeHelpers.h"
0042 
0043 #include "FWCore/Utilities/interface/isFinite.h"
0044 //#include "DataFormats/Math/interface/deltaR.h"
0045 
0046 //#define VTXDEBUG
0047 
0048 namespace svhelper {
0049   inline double cov33(const reco::Vertex &sv) { return sv.covariance(3, 3); }
0050   inline double cov33(const reco::VertexCompositePtrCandidate &sv) { return sv.vertexCovariance(3, 3); }
0051 }  // namespace svhelper
0052 
0053 template <class VTX>
0054 class TrackVertexArbitration {
0055 public:
0056   TrackVertexArbitration(const edm::ParameterSet &params);
0057 
0058   std::vector<VTX> trackVertexArbitrator(edm::Handle<reco::BeamSpot> &beamSpot,
0059                                          const reco::Vertex &pv,
0060                                          std::vector<reco::TransientTrack> &selectedTracks,
0061                                          std::vector<VTX> &secondaryVertices);
0062 
0063 private:
0064   bool trackFilterArbitrator(const reco::TransientTrack &track) const;
0065 
0066   edm::InputTag primaryVertexCollection;
0067   edm::InputTag secondaryVertexCollection;
0068   edm::InputTag trackCollection;
0069   edm::InputTag beamSpotCollection;
0070   double dRCut;
0071   double distCut;
0072   double sigCut;
0073   double dLenFraction;
0074   double fitterSigmacut;  // = 3.0;
0075   double fitterTini;      // = 256.;
0076   double fitterRatio;     //    = 0.25;
0077   int trackMinLayers;
0078   double trackMinPt;
0079   int trackMinPixels;
0080   double maxTimeSignificance;
0081   double sign;
0082 };
0083 
0084 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0085 template <class VTX>
0086 TrackVertexArbitration<VTX>::TrackVertexArbitration(const edm::ParameterSet &params)
0087     : primaryVertexCollection(params.getParameter<edm::InputTag>("primaryVertices")),
0088       secondaryVertexCollection(params.getParameter<edm::InputTag>("secondaryVertices")),
0089       trackCollection(params.getParameter<edm::InputTag>("tracks")),
0090       beamSpotCollection(params.getParameter<edm::InputTag>("beamSpot")),
0091       dRCut(params.getParameter<double>("dRCut")),
0092       distCut(params.getParameter<double>("distCut")),
0093       sigCut(params.getParameter<double>("sigCut")),
0094       dLenFraction(params.getParameter<double>("dLenFraction")),
0095       fitterSigmacut(params.getParameter<double>("fitterSigmacut")),
0096       fitterTini(params.getParameter<double>("fitterTini")),
0097       fitterRatio(params.getParameter<double>("fitterRatio")),
0098       trackMinLayers(params.getParameter<int32_t>("trackMinLayers")),
0099       trackMinPt(params.getParameter<double>("trackMinPt")),
0100       trackMinPixels(params.getParameter<int32_t>("trackMinPixels")),
0101       maxTimeSignificance(params.getParameter<double>("maxTimeSignificance")) {
0102   sign = 1.0;
0103   if (dRCut < 0) {
0104     sign = -1.0;
0105   }
0106   dRCut *= dRCut;
0107 }
0108 template <class VTX>
0109 bool TrackVertexArbitration<VTX>::trackFilterArbitrator(const reco::TransientTrack &track) const {
0110   if (!track.isValid())
0111     return false;
0112   if (track.track().hitPattern().trackerLayersWithMeasurement() < trackMinLayers)
0113     return false;
0114   if (track.track().pt() < trackMinPt)
0115     return false;
0116   if (track.track().hitPattern().numberOfValidPixelHits() < trackMinPixels)
0117     return false;
0118 
0119   return true;
0120 }
0121 
0122 inline float trackWeight(const reco::Vertex &sv, const reco::TransientTrack &track) {
0123   return sv.trackWeight(track.trackBaseRef());
0124 }
0125 inline float trackWeight(const reco::VertexCompositePtrCandidate &sv, const reco::TransientTrack &tt) {
0126   const reco::CandidatePtrTransientTrack *cptt =
0127       dynamic_cast<const reco::CandidatePtrTransientTrack *>(tt.basicTransientTrack());
0128   if (cptt == nullptr)
0129     edm::LogError("DynamicCastingFailed") << "Casting of TransientTrack to CandidatePtrTransientTrack failed!";
0130   else {
0131     const reco::CandidatePtr &c = cptt->candidate();
0132     if (std::find(sv.daughterPtrVector().begin(), sv.daughterPtrVector().end(), c) != sv.daughterPtrVector().end())
0133       return 1.0;
0134     else
0135       return 0;
0136   }
0137   return 0;
0138 }
0139 
0140 template <class VTX>
0141 std::vector<VTX> TrackVertexArbitration<VTX>::trackVertexArbitrator(edm::Handle<reco::BeamSpot> &beamSpot,
0142                                                                     const reco::Vertex &pv,
0143                                                                     std::vector<reco::TransientTrack> &selectedTracks,
0144                                                                     std::vector<VTX> &secondaryVertices) {
0145   using namespace reco;
0146 
0147   //std::cout << "PV: " << pv.position() << std::endl;
0148   VertexDistance3D dist;
0149   AdaptiveVertexFitter theAdaptiveFitter(GeometricAnnealing(fitterSigmacut, fitterTini, fitterRatio),
0150                                          DefaultLinearizationPointFinder(),
0151                                          KalmanVertexUpdator<5>(),
0152                                          KalmanVertexTrackCompatibilityEstimator<5>(),
0153                                          KalmanVertexSmoother());
0154 
0155   std::vector<VTX> recoVertices;
0156   VertexDistance3D vdist;
0157   GlobalPoint ppv(pv.position().x(), pv.position().y(), pv.position().z());
0158 
0159   std::unordered_map<unsigned int, Measurement1D> cachedIP;
0160   for (typename std::vector<VTX>::const_iterator sv = secondaryVertices.begin(); sv != secondaryVertices.end(); ++sv) {
0161     const double svTime(sv->t()), svTimeCov(svhelper::cov33(*sv));
0162     const bool svHasTime = svTimeCov > 0.;
0163 
0164     GlobalPoint ssv(sv->position().x(), sv->position().y(), sv->position().z());
0165     GlobalVector flightDir = ssv - ppv;
0166     //            std::cout << "Vertex : " << sv-secondaryVertices->begin() << " " << sv->position() << std::endl;
0167     Measurement1D dlen =
0168         vdist.distance(pv, VertexState(RecoVertex::convertPos(sv->position()), RecoVertex::convertError(sv->error())));
0169     std::vector<reco::TransientTrack> selTracks;
0170     for (unsigned int itrack = 0; itrack < selectedTracks.size(); itrack++) {
0171       TransientTrack &tt = (selectedTracks)[itrack];
0172       if (!trackFilterArbitrator(tt))
0173         continue;
0174       tt.setBeamSpot(*beamSpot);
0175       float w = trackWeight(*sv, tt);
0176       Measurement1D ipv;
0177       if (cachedIP.count(itrack)) {
0178         ipv = cachedIP[itrack];
0179       } else {
0180         std::pair<bool, Measurement1D> ipvp = IPTools::absoluteImpactParameter3D(tt, pv);
0181         cachedIP[itrack] = ipvp.second;
0182         ipv = ipvp.second;
0183       }
0184 
0185       AnalyticalImpactPointExtrapolator extrapolator(tt.field());
0186       TrajectoryStateOnSurface tsos =
0187           extrapolator.extrapolate(tt.impactPointState(), RecoVertex::convertPos(sv->position()));
0188       if (!tsos.isValid())
0189         continue;
0190       GlobalPoint refPoint = tsos.globalPosition();
0191       GlobalError refPointErr = tsos.cartesianError().position();
0192       Measurement1D isv =
0193           dist.distance(VertexState(RecoVertex::convertPos(sv->position()), RecoVertex::convertError(sv->error())),
0194                         VertexState(refPoint, refPointErr));
0195 
0196       float dR = Geom::deltaR2(((sign > 0) ? flightDir : -flightDir),
0197                                tt.track());  //.eta(), flightDir.phi(), tt.track().eta(), tt.track().phi());
0198 
0199       double timeSig = 0.;
0200       if (svHasTime && edm::isFinite(tt.timeExt())) {
0201         double tError = std::sqrt(std::pow(tt.dtErrorExt(), 2) + svTimeCov);
0202         timeSig = std::abs(tt.timeExt() - svTime) / tError;
0203       }
0204 
0205       if (w > 0 || (isv.significance() < sigCut && isv.value() < distCut && isv.value() < dlen.value() * dLenFraction &&
0206                     timeSig < maxTimeSignificance)) {
0207         if ((isv.value() < ipv.value()) && isv.value() < distCut && isv.value() < dlen.value() * dLenFraction &&
0208             dR < dRCut && timeSig < maxTimeSignificance) {
0209 #ifdef VTXDEBUG
0210           if (w > 0.5)
0211             std::cout << " = ";
0212           else
0213             std::cout << " + ";
0214 #endif
0215           selTracks.push_back(tt);
0216         } else {
0217 #ifdef VTXDEBUG
0218           if (w > 0.5 && isv.value() > ipv.value())
0219             std::cout << " - ";
0220           else
0221             std::cout << "   ";
0222 #endif
0223           //add also the tracks used in previous fitting that are still closer to Sv than Pv
0224           if (w > 0.5 && isv.value() <= ipv.value() && dR < dRCut && timeSig < maxTimeSignificance) {
0225             selTracks.push_back(tt);
0226 #ifdef VTXDEBUG
0227             std::cout << " = ";
0228 #endif
0229           }
0230           if (w > 0.5 && isv.value() <= ipv.value() && dR >= dRCut) {
0231 #ifdef VTXDEBUG
0232             std::cout << " - ";
0233 #endif
0234           }
0235         }
0236 #ifdef VTXDEBUG
0237 
0238         std::cout << "t : " << itrack << " ref " << ref.key() << " w: " << w << " svip: " << isv.significance() << " "
0239                   << isv.value() << " pvip: " << ipv.significance() << " " << ipv.value() << " res "
0240                   << tt.track().residualX(0) << "," << tt.track().residualY(0)
0241 //                  << " tpvip: " << itpv.second.significance() << " " << itpv.second.value()  << " dr: "   << dR << std::endl;
0242 #endif
0243 
0244       } else {
0245 #ifdef VTXDEBUG
0246 
0247         std::cout << " . t : " << itrack << " ref " << ref.key() << " w: " << w
0248                   << " svip: " << isv.second.significance() << " " << isv.second.value()
0249                   << " pvip: " << ipv.significance() << " " << ipv.value() << " dr: " << dR << std::endl;
0250 #endif
0251       }
0252     }
0253 
0254     if (selTracks.size() >= 2) {
0255       TransientVertex singleFitVertex;
0256       singleFitVertex = theAdaptiveFitter.vertex(selTracks, ssv);
0257 
0258       if (singleFitVertex.isValid()) {
0259         if (pv.covariance(3, 3) > 0.)
0260           svhelper::updateVertexTime(singleFitVertex);
0261         recoVertices.push_back(VTX(singleFitVertex));
0262 
0263 #ifdef VTXDEBUG
0264         const VTX &extVertex = recoVertices.back();
0265         GlobalVector vtxDir = GlobalVector(extVertex.p4().X(), extVertex.p4().Y(), extVertex.p4().Z());
0266 
0267         std::cout << " pointing : " << Geom::deltaR(extVertex.position() - pv.position(), vtxDir) << std::endl;
0268 #endif
0269       }
0270     }
0271   }
0272   return recoVertices;
0273 }
0274 
0275 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!