LinearizationPointFinders/src/CrossingPtBasedLinearizationPointFinder.cc

0001 #include "RecoVertex/LinearizationPointFinders/interface/CrossingPtBasedLinearizationPointFinder.h"
0002 #include "RecoVertex/LinearizationPointFinders/interface/LinPtException.h"
0003 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
0004 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0005 #include "RecoVertex/VertexTools/interface/ModeFinder3d.h"
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "RecoVertex/LinearizationPointFinders/interface/FallbackLinearizationPointFinder.h"
0008 #include <cmath>
0009 #include <algorithm>
0010
0011 // #define Use_GraphicsHarvester
0012 #ifdef Use_GraphicsHarvester
0013 #include "RecoVertex/VertexSimpleVis/interface/PrimitivesHarvester.h"
0014 #include "Utilities/UI/interface/SimpleConfigurable.h"
0015 #include "RecoVertex/VertexTools/interface/SubsetHsmModeFinder3d.h"
0016 #include "RecoVertex/VertexTools/interface/HsmModeFinder3d.h"
0017 #include "RecoVertex/VertexTools/interface/LmsModeFinder3d.h"
0018 #endif
0019
0020 namespace {
0021   inline GlobalPoint operator-(const GlobalPoint& a, const GlobalPoint& b) {
0022     return GlobalPoint(a.x() - b.x(), a.y() - b.y(), a.z() - b.z());
0023   }
0024
0025   inline GlobalPoint operator+(const GlobalPoint& a, const GlobalPoint& b) {
0026     return GlobalPoint(a.x() + b.x(), a.y() + b.y(), a.z() + b.z());
0027   }
0028
0029   inline GlobalPoint operator/(const GlobalPoint& a, const double b) {
0030     return GlobalPoint(a.x() / b, a.y() / b, a.z() / b);
0031   }
0032
0033 #ifndef __clang__
0034   inline GlobalPoint operator*(const GlobalPoint& a, const double b) {
0035     return GlobalPoint(a.x() * b, a.y() * b, a.z() * b);
0036   }
0037
0038   inline GlobalPoint operator*(const double b, const GlobalPoint& a) {
0039     return GlobalPoint(a.x() * b, a.y() * b, a.z() * b);
0040   }
0041 #endif
0042
0043   inline unsigned int sum(unsigned int nr) {
0044     /*
0045     int ret=0;
0046     for ( int i=1; i<= nr ; i++ )
0047     {
0048       ret+=i;
0049     }
0050     return ret;
0051     */
0052     return (nr * (nr + 1)) / 2;
0053   }
0054
0055 #ifdef Use_GraphicsHarvester
0056
0057   void graphicsDebug(const std::vector<PointAndDistance>& input) {
0058     bool sve = SimpleConfigurable<bool>(false, "CrossingPointStudy:Harvest").value();
0059     std::string fname = SimpleConfigurable<string>("crossingpoints.txt", "CrossingPointStudy:FileName").value();
0060     if (sve) {
0061       for (std::vector<PointAndDistance>::const_iterator i = input.begin(); i != input.end(); ++i) {
0062         std::map<std::string, MultiType> attrs;
0063         attrs["name"] = "GlobalPoint";
0064         attrs["color"] = "yellow";
0065         attrs["mag"] = .7;
0066         attrs["comment"] = "Input for CrossingPtBasedLinearizationPointFinder";
0067         PrimitivesHarvester(fname).save(i->first, attrs);
0068       }
0069       std::map<std::string, MultiType> attrs;
0070       attrs["name"] = "The Mode(*theAlgo)";
0071       attrs["color"] = "green";
0072       attrs["comment"] = "Output from CrossingPtBasedLinearizationPointFinder";
0073       PrimitivesHarvester(fname).save(ret, attrs);
0074       GlobalPoint subsethsm = SubsetHsmModeFinder3d()(input);
0075       attrs["name"] = "The Mode(SubsetHSM)";
0076       attrs["color"] = "red";
0077       PrimitivesHarvester(fname).save(subsethsm, attrs);
0078       GlobalPoint hsm = HsmModeFinder3d()(input);
0079       attrs["name"] = "The Mode(HSM)";
0080       attrs["color"] = "blue";
0081       PrimitivesHarvester(fname).save(hsm, attrs);
0082       GlobalPoint lms = LmsModeFinder3d()(input);
0083       attrs["name"] = "The Mode(LMS)";
0084       attrs["color"] = "gold";
0085       PrimitivesHarvester(fname).save(lms, attrs);
0086     }
0087   }
0088 #endif
0089 }  // namespace
0090
0091 CrossingPtBasedLinearizationPointFinder::CrossingPtBasedLinearizationPointFinder(const ModeFinder3d& algo,
0092                                                                                  const signed int n_pairs)
0093     : useMatrix(false), theNPairs(n_pairs), theMatrix(nullptr), theAlgo(algo.clone()) {}
0094
0095 CrossingPtBasedLinearizationPointFinder::CrossingPtBasedLinearizationPointFinder(const RecTracksDistanceMatrix* m,
0096                                                                                  const ModeFinder3d& algo,
0097                                                                                  const signed int n_pairs)
0098     : useMatrix(m->hasCrossingPoints()), theNPairs(n_pairs), theMatrix(m), theAlgo(algo.clone()) {}
0099
0100 CrossingPtBasedLinearizationPointFinder::CrossingPtBasedLinearizationPointFinder(
0101     const CrossingPtBasedLinearizationPointFinder& o)
0102     : useMatrix(o.useMatrix),
0103       theNPairs(o.theNPairs),
0104       theMatrix(o.theMatrix /* nope, we dont clone!! */),
0105       theAlgo(o.theAlgo->clone()) {}
0106
0107 CrossingPtBasedLinearizationPointFinder::~CrossingPtBasedLinearizationPointFinder() { delete theAlgo; }
0108
0109 std::vector<reco::TransientTrack> CrossingPtBasedLinearizationPointFinder::getBestTracks(
0110     const std::vector<reco::TransientTrack>& tracks) const {
0111   unsigned int n_tracks = (2 * (unsigned int)(theNPairs)) < tracks.size() ? 2 * theNPairs : tracks.size();
0112
0113   std::vector<reco::TransientTrack> newtracks = tracks;
0114   partial_sort(newtracks.begin(), newtracks.begin() + n_tracks, newtracks.end(), CompareTwoTracks());
0115   newtracks.erase(newtracks.begin() + n_tracks, newtracks.end());
0116
0117   /*
0118      * old version, when we have a default constructor
0119      *
0120
0121     std::vector <reco::TransientTrack> newtracks ( n_tracks );
0122     partial_sort_copy ( tracks.begin(), tracks.end(), newtracks.begin(),
0123                         newtracks.begin() + n_tracks  , CompareTwoTracks() );
0124                         */
0125
0126   return newtracks;
0127 }
0128
0129 typedef std::pair<GlobalPoint, float> PointAndDistance;
0130
0131 GlobalPoint CrossingPtBasedLinearizationPointFinder::useAllTracks(
0132     const std::vector<reco::TransientTrack>& tracks) const {
0133   std::vector<PointAndDistance> vgp;
0134   // vgp.reserve ( ( tracks.size() * ( tracks.size()-1 ) ) / 2 - 1 );
0135   TwoTrackMinimumDistance ttmd;
0136   bool status;
0137   std::vector<reco::TransientTrack>::const_iterator end = tracks.end();
0138   std::vector<reco::TransientTrack>::const_iterator endm1 = (end - 1);
0139   for (std::vector<reco::TransientTrack>::const_iterator x = tracks.begin(); x != endm1; ++x) {
0140     for (std::vector<reco::TransientTrack>::const_iterator y = x + 1; y != end; ++y) {
0141       status = ttmd.calculate((*x).impactPointState(), (*y).impactPointState());
0142       if (status) {
0143         std::pair<GlobalPoint, GlobalPoint> pts = ttmd.points();
0144         std::pair<GlobalPoint, float> v((pts.second + pts.first) / 2., (pts.second - pts.first).mag());
0145         vgp.push_back(v);
0146       }  // If sth goes wrong, we just dont add. Who cares?
0147     }
0148   }
0149   if (vgp.empty()) {
0150     return FallbackLinearizationPointFinder().getLinearizationPoint(tracks);
0151   }
0152   return find(vgp);
0153 }
0154
0155 GlobalPoint CrossingPtBasedLinearizationPointFinder::useFullMatrix(
0156     const std::vector<reco::TransientTrack>& tracks) const {
0157   std::vector<PointAndDistance> vgp;
0158   vgp.reserve((int)(tracks.size() * (tracks.size() - 1) / 2. - 1));
0159   std::vector<reco::TransientTrack>::const_iterator end = tracks.end();
0160   std::vector<reco::TransientTrack>::const_iterator endm1 = (end - 1);
0161   for (std::vector<reco::TransientTrack>::const_iterator x = tracks.begin(); x != endm1; ++x) {
0162     for (std::vector<reco::TransientTrack>::const_iterator y = x + 1; y != end; ++y) {
0163       PointAndDistance v(theMatrix->crossingPoint(*x, *y), theMatrix->distance(*x, *y));
0164       vgp.push_back(v);
0165     }
0166   }
0167   if (vgp.empty()) {
0168     return FallbackLinearizationPointFinder().getLinearizationPoint(tracks);
0169   }
0170   return find(vgp);
0171 }
0172
0173 GlobalPoint CrossingPtBasedLinearizationPointFinder::getLinearizationPoint(
0174     const std::vector<FreeTrajectoryState>& tracks) const {
0175   return LinearizationPointFinder::getLinearizationPoint(tracks);
0176 }
0177
0178 GlobalPoint CrossingPtBasedLinearizationPointFinder::find(const std::vector<PointAndDistance>& input) const {
0179   /*
0180     std::cout << "[CrossingPtBasedLinearizationPointFinder] input size="
0181          << input.size() << std::endl;*/
0182   GlobalPoint ret = (*theAlgo)(input);
0183 #ifdef Use_GraphicsHarvester
0184
0185   graphicsDebug(input);
0186 #endif
0187
0188   return ret;
0189 }
0190
0191 GlobalPoint CrossingPtBasedLinearizationPointFinder::getLinearizationPoint(
0192     const std::vector<reco::TransientTrack>& tracks) const {
0193   if (tracks.size() < 2)
0194     throw LinPtException("CrossingPtBasedLinPtFinder: too few tracks given.");
0195   std::vector<PointAndDistance> vgp;
0196   if (theNPairs == -1) {
0197     if (useMatrix) {
0198       return useFullMatrix(tracks);
0199     } else {
0200       return useAllTracks(tracks);
0201     }
0202   }
0203
0204   if (sum(tracks.size() - 1) < ((unsigned int)(theNPairs))) {
0205     /*
0206             std::cout << "[CrossingPtBasedLinearizationPointFinder] we exploit all track pairs"
0207                  << std::endl;*/
0208     // we have fewer track pair options than is foreseen
0209     // in the quota.
0210     // so we exploit all tracks pairs.
0211     return useAllTracks(tracks);
0212   }
0213
0214   std::vector<reco::TransientTrack> goodtracks = getBestTracks(tracks);
0215
0216   // we sort according to momenta.
0217   if (goodtracks.size() < 2)
0218     throw LinPtException("CrossingPtBasedLinPtFinder: less than two tracks given");
0219   // vgp.reserve ( theNPairs - 1 );
0220   unsigned int t_first = 0;
0221   unsigned int t_interval = goodtracks.size() / 2;
0222   unsigned int lim = goodtracks.size() - 1;
0223
0224   /*
0225         std::cout << "[CrossingPtBasedLinearizationPointFinder] we start: npairs=" << theNPairs
0226              << std::endl;
0227         std::cout << "[CrossingPtBasedLinearizationPointFinder] t_interval=" << t_interval << std::endl;
0228         std::cout << "[CrossingPtBasedLinearizationPointFinder goodtracks.size=" << goodtracks.size()
0229              << std::endl;*/
0230
0231   // the 'direction' false: intervals will expand
0232   // true: intervals will shrink
0233   bool dir = false;
0234
0235   while (vgp.size() < ((unsigned int)(theNPairs))) {
0236     reco::TransientTrack rt1 = goodtracks[t_first];
0237     reco::TransientTrack rt2 = goodtracks[t_first + t_interval];
0238     // std::cout << "Considering now: " << t_first << ", " << t_first+t_interval << std::endl;
0239     if (useMatrix) {
0240       PointAndDistance v(theMatrix->crossingPoint(rt1, rt2), theMatrix->distance(rt1, rt2));
0241       vgp.push_back(v);
0242     } else {  // No DistanceMatrix available
0243       TwoTrackMinimumDistance ttmd;
0244       bool status = ttmd.calculate(rt1.impactPointState(), rt2.impactPointState());
0245       if (status) {
0246         std::pair<GlobalPoint, GlobalPoint> pts = ttmd.points();
0247         PointAndDistance v((pts.second + pts.first) / 2., (pts.second - pts.first).mag());
0248         vgp.push_back(v);
0249       }
0250     }
0251     if ((t_first + t_interval) < lim) {
0252       t_first++;
0253     } else if (dir) {
0254       t_first = 0;
0255       t_interval--;
0256       if (t_interval == 0) {
0257         /* std::cout << "[CrossingPtBasedLinearizationPointFinder] t_interval=0. break."
0258                          << std::endl;*/
0259         break;
0260       }
0261     } else {
0262       t_first = 0;
0263       t_interval++;
0264       if (t_interval == goodtracks.size()) {
0265         /* std::cout << "[CrossingPtBasedLinearizationPointFinder] t_interval="
0266                          << goodtracks.size() << "(max). start to decrease intervals"
0267                          << std::endl; */
0268         dir = true;
0269         t_interval = goodtracks.size() / 2 - 1;
0270       }
0271     }
0272   }
0273   if (vgp.empty()) {
0274     // no crossing points? Fallback to a crossingpoint-less lin pt finder!
0275     return FallbackLinearizationPointFinder().getLinearizationPoint(tracks);
0276   }
0277   return find(vgp);
0278
0279   return GlobalPoint(0., 0., 0.);  // if nothing else, then return 0,0,0
0280 }