VertexPrimitives/src/CachingVertex.cc

0001 #include "RecoVertex/VertexPrimitives/interface/CachingVertex.h"
0002 #include "RecoVertex/VertexPrimitives/interface/LinearizedTrackState.h"
0003 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0004 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0005 #include "RecoVertex/VertexPrimitives/interface/TTtoTTmap.h"
0006 #include <map>
0007
0008 //to be removed
0009 template <unsigned int N>
0010 CachingVertex<N>::CachingVertex(const GlobalPoint& pos,
0011                                 const GlobalError& posErr,
0012                                 const std::vector<RefCountedVertexTrack>& tks,
0013                                 float totalChiSq)
0014     : theVertexState(pos, posErr),
0015       theChiSquared(totalChiSq),
0016       theNDF(0),
0017       theNDFAvailable(false),
0018       theTracks(tks),
0019       theCovMapAvailable(false),
0020       withPrior(false),
0021       theValid(true),
0022       vertexIs4D(false)
0023
0024 {}
0025
0026 //to be removed
0027 template <unsigned int N>
0028 CachingVertex<N>::CachingVertex(const GlobalPoint& pos,
0029                                 const double time,
0030                                 const GlobalError& posTimeErr,
0031                                 const std::vector<RefCountedVertexTrack>& tks,
0032                                 float totalChiSq)
0033     : theVertexState(pos, time, posTimeErr),
0034       theChiSquared(totalChiSq),
0035       theNDF(0),
0036       theNDFAvailable(false),
0037       theTracks(tks),
0038       theCovMapAvailable(false),
0039       withPrior(false),
0040       theValid(true),
0041       vertexIs4D(true)
0042
0043 {}
0044
0045 //to be removed
0046 template <unsigned int N>
0047 CachingVertex<N>::CachingVertex(const GlobalPoint& pos,
0048                                 const GlobalWeight& posWeight,
0049                                 const std::vector<RefCountedVertexTrack>& tks,
0050                                 float totalChiSq)
0051     : theVertexState(pos, posWeight),
0052       theChiSquared(totalChiSq),
0053       theNDF(0),
0054       theNDFAvailable(false),
0055       theTracks(tks),
0056       theCovMapAvailable(false),
0057       withPrior(false),
0058       theValid(true),
0059       vertexIs4D(false) {}
0060
0061 //to be removed
0062 template <unsigned int N>
0063 CachingVertex<N>::CachingVertex(const GlobalPoint& pos,
0064                                 const double time,
0065                                 const GlobalWeight& posTimeWeight,
0066                                 const std::vector<RefCountedVertexTrack>& tks,
0067                                 float totalChiSq)
0068     : theVertexState(pos, time, posTimeWeight),
0069       theChiSquared(totalChiSq),
0070       theNDF(0),
0071       theNDFAvailable(false),
0072       theTracks(tks),
0073       theCovMapAvailable(false),
0074       withPrior(false),
0075       theValid(true),
0076       vertexIs4D(true) {}
0077
0078 //to be removed
0079 template <unsigned int N>
0080 CachingVertex<N>::CachingVertex(const AlgebraicVector3& weightTimesPosition,
0081                                 const GlobalWeight& posWeight,
0082                                 const std::vector<RefCountedVertexTrack>& tks,
0083                                 float totalChiSq)
0084     : theVertexState(weightTimesPosition, posWeight),
0085       theChiSquared(totalChiSq),
0086       theNDF(0),
0087       theNDFAvailable(false),
0088       theTracks(tks),
0089       theCovMapAvailable(false),
0090       withPrior(false),
0091       theValid(true),
0092       vertexIs4D(false) {}
0093
0094 //to be removed
0095 template <unsigned int N>
0096 CachingVertex<N>::CachingVertex(const AlgebraicVector4& weightTimesPosition,
0097                                 const GlobalWeight& posWeight,
0098                                 const std::vector<RefCountedVertexTrack>& tks,
0099                                 float totalChiSq)
0100     : theVertexState(weightTimesPosition, posWeight),
0101       theChiSquared(totalChiSq),
0102       theNDF(0),
0103       theNDFAvailable(false),
0104       theTracks(tks),
0105       theCovMapAvailable(false),
0106       withPrior(false),
0107       theValid(true),
0108       vertexIs4D(false) {}
0109
0110 template <unsigned int N>
0111 CachingVertex<N>::CachingVertex(const VertexState& aVertexState,
0112                                 const std::vector<RefCountedVertexTrack>& tks,
0113                                 float totalChiSq)
0114     : theVertexState(aVertexState),
0115       theChiSquared(totalChiSq),
0116       theNDF(0),
0117       theNDFAvailable(false),
0118       theTracks(tks),
0119       theCovMapAvailable(false),
0120       withPrior(false),
0121       theValid(true),
0122       vertexIs4D(aVertexState.is4D()) {}
0123
0124 template <unsigned int N>
0125 CachingVertex<N>::CachingVertex(const VertexState& aVertexState,
0126                                 const std::vector<RefCountedVertexTrack>& tks,
0127                                 float totalChiSq,
0128                                 const TrackToTrackMap& covMap)
0129     : theVertexState(aVertexState),
0130       theChiSquared(totalChiSq),
0131       theNDF(0),
0132       theNDFAvailable(false),
0133       theTracks(tks),
0134       theCovMap(covMap),
0135       theCovMapAvailable(true),
0136       withPrior(false),
0137       theValid(true),
0138       vertexIs4D(aVertexState.is4D()) {
0139   if (theCovMap.empty())
0140     theCovMapAvailable = false;
0141 }
0142
0143 template <unsigned int N>
0144 CachingVertex<N>::CachingVertex(const VertexState& priorVertexState,
0145                                 const VertexState& aVertexState,
0146                                 const std::vector<RefCountedVertexTrack>& tks,
0147                                 float totalChiSq)
0148     : theVertexState(aVertexState),
0149       theChiSquared(totalChiSq),
0150       theNDF(0),
0151       theNDFAvailable(false),
0152       theTracks(tks),
0153       theCovMapAvailable(false),
0154       thePriorVertexState(priorVertexState),
0155       withPrior(true),
0156       theValid(true),
0157       vertexIs4D(priorVertexState.is4D() && aVertexState.is4D()) {}
0158
0159 //to be removed
0160 template <unsigned int N>
0161 CachingVertex<N>::CachingVertex(const GlobalPoint& priorPos,
0162                                 const GlobalError& priorErr,
0163                                 const GlobalPoint& pos,
0164                                 const GlobalError& posErr,
0165                                 const std::vector<RefCountedVertexTrack>& tks,
0166                                 float totalChiSq)
0167     : theVertexState(pos, posErr),
0168       theChiSquared(totalChiSq),
0169       theNDF(0),
0170       theNDFAvailable(false),
0171       theTracks(tks),
0172       theCovMapAvailable(false),
0173       thePriorVertexState(priorPos, priorErr),
0174       withPrior(true),
0175       theValid(true),
0176       vertexIs4D(false) {}
0177
0178 //to be removed
0179 template <unsigned int N>
0180 CachingVertex<N>::CachingVertex(const GlobalPoint& priorPos,
0181                                 const double priorTime,
0182                                 const GlobalError& priorErr,
0183                                 const GlobalPoint& pos,
0184                                 const double time,
0185                                 const GlobalError& posErr,
0186                                 const std::vector<RefCountedVertexTrack>& tks,
0187                                 float totalChiSq)
0188     : theVertexState(pos, time, posErr),
0189       theChiSquared(totalChiSq),
0190       theNDF(0),
0191       theNDFAvailable(false),
0192       theTracks(tks),
0193       theCovMapAvailable(false),
0194       thePriorVertexState(priorPos, priorTime, priorErr),
0195       withPrior(true),
0196       theValid(true),
0197       vertexIs4D(true) {}
0198
0199 //to be removed
0200 template <unsigned int N>
0201 CachingVertex<N>::CachingVertex(const GlobalPoint& priorPos,
0202                                 const GlobalError& priorErr,
0203                                 const GlobalPoint& pos,
0204                                 const GlobalWeight& posWeight,
0205                                 const std::vector<RefCountedVertexTrack>& tks,
0206                                 float totalChiSq)
0207     : theVertexState(pos, posWeight),
0208       theChiSquared(totalChiSq),
0209       theNDF(0),
0210       theNDFAvailable(false),
0211       theTracks(tks),
0212       theCovMapAvailable(false),
0213       thePriorVertexState(priorPos, priorErr),
0214       withPrior(true),
0215       theValid(true),
0216       vertexIs4D(false) {}
0217
0218 //to be removed
0219 template <unsigned int N>
0220 CachingVertex<N>::CachingVertex(const GlobalPoint& priorPos,
0221                                 const double priorTime,
0222                                 const GlobalError& priorErr,
0223                                 const GlobalPoint& pos,
0224                                 const double time,
0225                                 const GlobalWeight& posWeight,
0226                                 const std::vector<RefCountedVertexTrack>& tks,
0227                                 float totalChiSq)
0228     : theVertexState(pos, priorTime, posWeight),
0229       theChiSquared(totalChiSq),
0230       theNDF(0),
0231       theNDFAvailable(false),
0232       theTracks(tks),
0233       theCovMapAvailable(false),
0234       thePriorVertexState(priorPos, priorTime, priorErr),
0235       withPrior(true),
0236       theValid(true),
0237       vertexIs4D(true) {}
0238
0239 //to be removed
0240 template <unsigned int N>
0241 CachingVertex<N>::CachingVertex(const GlobalPoint& priorPos,
0242                                 const GlobalError& priorErr,
0243                                 const AlgebraicVector3& weightTimesPosition,
0244                                 const GlobalWeight& posWeight,
0245                                 const std::vector<RefCountedVertexTrack>& tks,
0246                                 float totalChiSq)
0247     : theVertexState(weightTimesPosition, posWeight),
0248       theChiSquared(totalChiSq),
0249       theNDF(0),
0250       theNDFAvailable(false),
0251       theTracks(tks),
0252       theCovMapAvailable(false),
0253       thePriorVertexState(priorPos, priorErr),
0254       withPrior(true),
0255       theValid(true) {}
0256
0257 //to be removed
0258 template <unsigned int N>
0259 CachingVertex<N>::CachingVertex(const GlobalPoint& priorPos,
0260                                 const GlobalError& priorErr,
0261                                 const AlgebraicVector4& weightTimesPosition,
0262                                 const GlobalWeight& posWeight,
0263                                 const std::vector<RefCountedVertexTrack>& tks,
0264                                 float totalChiSq)
0265     : theVertexState(weightTimesPosition, posWeight),
0266       theChiSquared(totalChiSq),
0267       theNDF(0),
0268       theNDFAvailable(false),
0269       theTracks(tks),
0270       theCovMapAvailable(false),
0271       thePriorVertexState(priorPos, priorErr),
0272       withPrior(true),
0273       theValid(true) {}
0274
0275 template <unsigned int N>
0276 CachingVertex<N>::CachingVertex(const VertexState& priorVertexState,
0277                                 const VertexState& aVertexState,
0278                                 const std::vector<RefCountedVertexTrack>& tks,
0279                                 float totalChiSq,
0280                                 const TrackToTrackMap& covMap)
0281     : theVertexState(aVertexState),
0282       theChiSquared(totalChiSq),
0283       theNDF(0),
0284       theNDFAvailable(false),
0285       theTracks(tks),
0286       theCovMap(covMap),
0287       theCovMapAvailable(true),
0288       thePriorVertexState(priorVertexState),
0289       withPrior(true),
0290       theValid(true),
0291       vertexIs4D(priorVertexState.is4D() && aVertexState.is4D()) {
0292   if (theCovMap.empty())
0293     theCovMapAvailable = false;
0294 }
0295
0296 template <unsigned int N>
0297 CachingVertex<N>::CachingVertex()
0298     : theChiSquared(-1),
0299       theNDF(0),
0300       theNDFAvailable(false),
0301       theTracks(),
0302       theCovMapAvailable(false),
0303       withPrior(false),
0304       theValid(false),
0305       vertexIs4D(false) {}
0306
0307 template <unsigned int N>
0308 GlobalPoint CachingVertex<N>::position() const {
0309   return theVertexState.position();
0310 }
0311
0312 template <unsigned int N>
0313 double CachingVertex<N>::time() const {
0314   return theVertexState.time();
0315 }
0316
0317 template <unsigned int N>
0318 GlobalError CachingVertex<N>::error() const {
0319   return theVertexState.error();
0320 }
0321
0322 template <unsigned int N>
0323 GlobalError CachingVertex<N>::error4D() const {
0324   return theVertexState.error4D();
0325 }
0326
0327 template <unsigned int N>
0328 GlobalWeight CachingVertex<N>::weight() const {
0329   return theVertexState.weight();
0330 }
0331
0332 template <unsigned int N>
0333 GlobalWeight CachingVertex<N>::weight4D() const {
0334   return theVertexState.weight4D();
0335 }
0336
0337 template <unsigned int N>
0338 AlgebraicVector3 CachingVertex<N>::weightTimesPosition() const {
0339   return theVertexState.weightTimesPosition();
0340 }
0341
0342 template <unsigned int N>
0343 AlgebraicVector4 CachingVertex<N>::weightTimesPosition4D() const {
0344   return theVertexState.weightTimesPosition4D();
0345 }
0346
0347 template <unsigned int N>
0348 float CachingVertex<N>::degreesOfFreedom() const {
0349   if (!theNDFAvailable)
0350     computeNDF();
0351   return theNDF;
0352 }
0353
0354 template <unsigned int N>
0355 void CachingVertex<N>::computeNDF() const {
0356   theNDF = 0;
0357   for (typename std::vector<RefCountedVertexTrack>::const_iterator itk = theTracks.begin(); itk != theTracks.end();
0358        ++itk) {
0359     theNDF += (**itk).weight();  // adds up weights
0360   }
0361   theNDF *= 2.;                                // times 2df for each track
0362   const double adjust = (vertexIs4D ? 4 : 3);  // ndf adjust is 3 or 4
0363   if (!withPrior)
0364     theNDF -= adjust;
0365   theNDFAvailable = true;
0366 }
0367
0368 template <unsigned int N>
0369 typename CachingVertex<N>::AlgebraicMatrixMM CachingVertex<N>::tkToTkCovariance(const RefCountedVertexTrack t1,
0370                                                                                 const RefCountedVertexTrack t2) const {
0371   if (!tkToTkCovarianceIsAvailable()) {
0372     throw VertexException("CachingVertex::TkTkCovariance requested before been calculated");
0373   } else {
0374     RefCountedVertexTrack tr1;
0375     RefCountedVertexTrack tr2;
0376     bool transp = false;
0377     if (t1 < t2) {
0378       tr1 = t1;
0379       tr2 = t2;
0380     } else {
0381       tr1 = t2;
0382       tr2 = t1;
0383       transp = true;
0384     }
0385     typename TrackToTrackMap::const_iterator it = theCovMap.find(tr1);
0386     if (it != theCovMap.end()) {
0387       const TrackMap& tm = it->second;
0388       typename TrackMap::const_iterator nit = tm.find(tr2);
0389       if (nit != tm.end()) {
0390         if (transp)
0391           return (ROOT::Math::Transpose(nit->second));
0392         else
0393           return (nit->second);
0394       } else {
0395         throw VertexException("CachingVertex::requested TkTkCovariance does not exist");
0396       }
0397     } else {
0398       throw VertexException("CachingVertex::requested TkTkCovariance does not exist");
0399     }
0400   }
0401 }
0402
0403 template <unsigned int N>
0404 CachingVertex<N>::operator TransientVertex() const {
0405   //If the vertex is invalid, return an invalid TV !
0406   if (!isValid())
0407     return TransientVertex();
0408
0409   typedef std::map<reco::TransientTrack, float> TransientTrackToFloatMap;
0410
0411   // Construct Track vector
0412   std::vector<reco::TransientTrack> ttVect;
0413   ttVect.reserve(theTracks.size());
0414   std::vector<reco::TransientTrack> refTTVect;
0415   TransientTrackToFloatMap theWeightMap;
0416   TTtoTTmap ttCovMap;
0417   // float theMinWeight = 0.5;
0418
0419   for (typename std::vector<RefCountedVertexTrack>::const_iterator i = theTracks.begin(); i != theTracks.end(); ++i) {
0420     // discard tracks with too low weight
0421     // if ((**i).weight() < theMinWeight) continue;
0422
0423     reco::TransientTrack t1((**i).linearizedTrack()->track());
0424     ttVect.push_back(t1);
0425     //Fill in the weight map
0426     theWeightMap[t1] = (**i).weight();
0427
0428     //Fill in the tk-to-tk covariance map
0429     if (theCovMapAvailable) {
0430       for (typename std::vector<RefCountedVertexTrack>::const_iterator j = (i + 1); j != theTracks.end(); ++j) {
0431         reco::TransientTrack t2((**j).linearizedTrack()->track());
0432         ttCovMap[t1][t2] = tkToTkCovariance(*i, *j);
0433       }
0434     }
0435     if ((**i).refittedStateAvailable()) {
0436       refTTVect.push_back((**i).refittedState()->transientTrack());
0437     }
0438   }
0439   TransientVertex tv;
0440   if (withPrior) {
0441     tv = TransientVertex(priorVertexState(), vertexState(), ttVect, totalChiSquared(), degreesOfFreedom());
0442   } else {
0443     tv = TransientVertex(vertexState(), ttVect, totalChiSquared(), degreesOfFreedom());
0444   }
0445   tv.weightMap(theWeightMap);
0446   if (theCovMapAvailable)
0447     tv.tkToTkCovariance(ttCovMap);
0448   if (!refTTVect.empty())
0449     tv.refittedTracks(refTTVect);
0450   return tv;
0451 }
0452
0453 template class CachingVertex<5>;
0454 template class CachingVertex<6>;