Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:23:11 

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