AdaptiveVertexFinder/src/TracksClusteringFromDisplacedSeed.cc

0001 #include "RecoVertex/AdaptiveVertexFinder/interface/TracksClusteringFromDisplacedSeed.h"
0002 //#define VTXDEBUG 1
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004
0005 TracksClusteringFromDisplacedSeed::TracksClusteringFromDisplacedSeed(const edm::ParameterSet &params)
0006     :  //   maxNTracks(params.getParameter<unsigned int>("maxNTracks")),
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")),  //3
0013       distanceRatio(params.getParameter<double>("distanceRatio")),                    //was clusterScale/densityFactor
0014       clusterMinAngleCosine(params.getParameter<double>("clusterMinAngleCosine")),    //0.0
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         // apply only if time available and being used in vertexing
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       //SK:UNUSED//    float dotprodTrackSeed2D = trackDir2D.unit().dot(seedDir2D.unit());
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)) &&  //Angles between PV-PCAonSeed vectors and seed directions
0068            ((clusterMinAngleCosine > 0)
0069                 ? (dotprodTrack > clusterMinAngleCosine)
0070                 : (dotprodTrack < clusterMinAngleCosine)) &&  //Angles between PV-PCAonTrack vectors and track directions
0071            //dotprodTrackSeed2D > clusterMinAngleCosine && //Angle between track and seed
0072            //distance*clusterScale*tracks.size() < (distanceFromPV+pvDistance)*(distanceFromPV+pvDistance)/pvDistance && // cut scaling with track density
0073            distance * distanceRatio < distanceFromPV &&  // cut scaling with track density
0074            distance < clusterMaxDistance &&
0075            timeSig < maxTimeSignificance);  // absolute distance cut
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;  // cut scaling with track density
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  // VTXDEBUG
0124     std::pair<std::vector<reco::TransientTrack>, GlobalPoint> ntracks = nearTracks(*s, selectedTracks, pv);
0125     //          std::cout << ntracks.first.size() << " " << ntracks.first.size()  << std::endl;
0126     //                if(ntracks.first.size() == 0 || ntracks.first.size() > maxNTracks ) continue;
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 }