File indexing completed on 2023-03-17 11:23:11
0001 #include "RecoVertex/AdaptiveVertexFinder/interface/TracksClusteringFromDisplacedSeed.h"
0002
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004
0005 TracksClusteringFromDisplacedSeed::TracksClusteringFromDisplacedSeed(const edm::ParameterSet ¶ms)
0006 :
0007 max3DIPSignificance(params.getParameter<double>("seedMax3DIPSignificance")),
0008 max3DIPValue(params.getParameter<double>("seedMax3DIPValue")),
0009 min3DIPSignificance(params.getParameter<double>("seedMin3DIPSignificance")),
0010 min3DIPValue(params.getParameter<double>("seedMin3DIPValue")),
0011 clusterMaxDistance(params.getParameter<double>("clusterMaxDistance")),
0012 clusterMaxSignificance(params.getParameter<double>("clusterMaxSignificance")),
0013 distanceRatio(params.getParameter<double>("distanceRatio")),
0014 clusterMinAngleCosine(params.getParameter<double>("clusterMinAngleCosine")),
0015 maxTimeSignificance(params.getParameter<double>("maxTimeSignificance"))
0016
0017 {}
0018
0019 std::pair<std::vector<reco::TransientTrack>, GlobalPoint> TracksClusteringFromDisplacedSeed::nearTracks(
0020 const reco::TransientTrack &seed,
0021 const std::vector<reco::TransientTrack> &tracks,
0022 const reco::Vertex &primaryVertex) const {
0023 VertexDistance3D distanceComputer;
0024 GlobalPoint pv(primaryVertex.position().x(), primaryVertex.position().y(), primaryVertex.position().z());
0025 std::vector<reco::TransientTrack> result;
0026 TwoTrackMinimumDistance dist;
0027 GlobalPoint seedingPoint;
0028 float sumWeights = 0;
0029 std::pair<bool, Measurement1D> ipSeed = IPTools::absoluteImpactParameter3D(seed, primaryVertex);
0030 float pvDistance = ipSeed.second.value();
0031 for (std::vector<reco::TransientTrack>::const_iterator tt = tracks.begin(); tt != tracks.end(); ++tt) {
0032 if (*tt == seed)
0033 continue;
0034
0035 if (dist.calculate(tt->impactPointState(), seed.impactPointState())) {
0036 GlobalPoint ttPoint = dist.points().first;
0037 GlobalError ttPointErr = tt->impactPointState().cartesianError().position();
0038 GlobalPoint seedPosition = dist.points().second;
0039 GlobalError seedPositionErr = seed.impactPointState().cartesianError().position();
0040 Measurement1D m =
0041 distanceComputer.distance(VertexState(seedPosition, seedPositionErr), VertexState(ttPoint, ttPointErr));
0042 GlobalPoint cp(dist.crossingPoint());
0043
0044 double timeSig = 0.;
0045 if (primaryVertex.covariance(3, 3) > 0. && edm::isFinite(seed.timeExt()) && edm::isFinite(tt->timeExt())) {
0046
0047 const double tError = std::sqrt(std::pow(seed.dtErrorExt(), 2) + std::pow(tt->dtErrorExt(), 2));
0048 timeSig = std::abs(seed.timeExt() - tt->timeExt()) / tError;
0049 }
0050
0051 float distanceFromPV = (dist.points().second - pv).mag();
0052 float distance = dist.distance();
0053 GlobalVector trackDir2D(
0054 tt->impactPointState().globalDirection().x(), tt->impactPointState().globalDirection().y(), 0.);
0055 GlobalVector seedDir2D(
0056 seed.impactPointState().globalDirection().x(), seed.impactPointState().globalDirection().y(), 0.);
0057
0058
0059 float dotprodTrack = (dist.points().first - pv).unit().dot(tt->impactPointState().globalDirection().unit());
0060 float dotprodSeed = (dist.points().second - pv).unit().dot(seed.impactPointState().globalDirection().unit());
0061
0062 float w = distanceFromPV * distanceFromPV / (pvDistance * distance);
0063 bool selected =
0064 (m.significance() < clusterMaxSignificance &&
0065 ((clusterMinAngleCosine > 0)
0066 ? (dotprodSeed > clusterMinAngleCosine)
0067 : (dotprodSeed < clusterMinAngleCosine)) &&
0068 ((clusterMinAngleCosine > 0)
0069 ? (dotprodTrack > clusterMinAngleCosine)
0070 : (dotprodTrack < clusterMinAngleCosine)) &&
0071
0072
0073 distance * distanceRatio < distanceFromPV &&
0074 distance < clusterMaxDistance &&
0075 timeSig < maxTimeSignificance);
0076
0077 #ifdef VTXDEBUG
0078 std::cout << tt->trackBaseRef().key() << " : " << (selected ? "+" : " ") << " " << m.significance() << " < "
0079 << clusterMaxSignificance << " && " << dotprodSeed << " > " << clusterMinAngleCosine << " && "
0080 << dotprodTrack << " > " << clusterMinAngleCosine << " && " << dotprodTrackSeed2D << " > "
0081 << clusterMinAngleCosine << " && " << distance * distanceRatio << " < " << distanceFromPV
0082 << " crossingtoPV: " << distanceFromPV << " dis*scal " << distance * distanceRatio << " < "
0083 << distanceFromPV << " dist: " << distance << " < " << clusterMaxDistance < < < <
0084 "timeSig: " << timeSig << std::endl;
0085 #endif
0086 if (selected) {
0087 result.push_back(*tt);
0088 seedingPoint =
0089 GlobalPoint(cp.x() * w + seedingPoint.x(), cp.y() * w + seedingPoint.y(), cp.z() * w + seedingPoint.z());
0090 sumWeights += w;
0091 }
0092 }
0093 }
0094
0095 seedingPoint =
0096 GlobalPoint(seedingPoint.x() / sumWeights, seedingPoint.y() / sumWeights, seedingPoint.z() / sumWeights);
0097 return std::pair<std::vector<reco::TransientTrack>, GlobalPoint>(result, seedingPoint);
0098 }
0099
0100 std::vector<TracksClusteringFromDisplacedSeed::Cluster> TracksClusteringFromDisplacedSeed::clusters(
0101 const reco::Vertex &pv, const std::vector<reco::TransientTrack> &selectedTracks) {
0102 using namespace reco;
0103 std::vector<TransientTrack> seeds;
0104 for (std::vector<TransientTrack>::const_iterator it = selectedTracks.begin(); it != selectedTracks.end(); it++) {
0105 std::pair<bool, Measurement1D> ip = IPTools::absoluteImpactParameter3D(*it, pv);
0106 if (ip.first && ip.second.value() >= min3DIPValue && ip.second.significance() >= min3DIPSignificance &&
0107 ip.second.value() <= max3DIPValue && ip.second.significance() <= max3DIPSignificance) {
0108 #ifdef VTXDEBUG
0109 std::cout << "new seed " << it - selectedTracks.begin() << " ref " << it->trackBaseRef().key() << " "
0110 << ip.second.value() << " " << ip.second.significance() << " "
0111 << it->track().hitPattern().trackerLayersWithMeasurement() << " " << it->track().pt() << " "
0112 << it->track().eta() << std::endl;
0113 #endif
0114 seeds.push_back(*it);
0115 }
0116 }
0117
0118 std::vector<Cluster> clusters;
0119 int i = 0;
0120 for (std::vector<TransientTrack>::const_iterator s = seeds.begin(); s != seeds.end(); ++s, ++i) {
0121 #ifdef VTXDEBUG
0122 std::cout << "Seed N. " << i << std::endl;
0123 #endif
0124 std::pair<std::vector<reco::TransientTrack>, GlobalPoint> ntracks = nearTracks(*s, selectedTracks, pv);
0125
0126
0127 ntracks.first.push_back(*s);
0128 Cluster aCl;
0129 aCl.seedingTrack = *s;
0130 aCl.seedPoint = ntracks.second;
0131 aCl.tracks = ntracks.first;
0132 clusters.push_back(aCl);
0133 }
0134
0135 return clusters;
0136 }