VertexPrimitives/src/BasicSingleVertexState.cc

0001 #include "RecoVertex/VertexPrimitives/interface/BasicSingleVertexState.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003 #include <limits>
0004
0005 namespace {
0006   constexpr double dNaN = std::numeric_limits<double>::quiet_NaN();
0007 }
0008
0009 BasicSingleVertexState::BasicSingleVertexState()
0010     : thePos(GlobalPoint(0, 0, 0)),
0011       theTime(dNaN),
0012       theErr(AlgebraicSymMatrix44()),
0013       theWeight(AlgebraicSymMatrix44()),
0014       theWeightTimesPos(AlgebraicVector4()),
0015       theWeightInMix(0.),
0016       thePosAvailable(false),
0017       theTimeAvailable(false),
0018       theErrAvailable(false),
0019       theWeightAvailable(false),
0020       theWeightTimesPosAvailable(false),
0021       valid(false),
0022       vertexIs4D(false) {}
0023
0024 BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
0025                                                const GlobalError& posErr,
0026                                                const double& weightInMix)
0027     : thePos(pos),
0028       theTime(dNaN),
0029       theErr(posErr),
0030       theWeight(AlgebraicSymMatrix44()),
0031       theWeightTimesPos(AlgebraicVector4()),
0032       theWeightInMix(weightInMix),
0033       thePosAvailable(true),
0034       theTimeAvailable(false),
0035       theErrAvailable(true),
0036       theWeightAvailable(false),
0037       theWeightTimesPosAvailable(false),
0038       valid(true),
0039       vertexIs4D(false) {}
0040
0041 BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
0042                                                const GlobalWeight& posWeight,
0043                                                const double& weightInMix)
0044     : thePos(pos),
0045       theTime(dNaN),
0046       theErr(AlgebraicSymMatrix44()),
0047       theWeight(posWeight),
0048       theWeightTimesPos(AlgebraicVector4()),
0049       theWeightInMix(weightInMix),
0050       thePosAvailable(true),
0051       theTimeAvailable(false),
0052       theErrAvailable(false),
0053       theWeightAvailable(true),
0054       theWeightTimesPosAvailable(false),
0055       valid(true),
0056       vertexIs4D(false) {}
0057
0058 BasicSingleVertexState::BasicSingleVertexState(const AlgebraicVector3& weightTimesPosition,
0059                                                const GlobalWeight& posWeight,
0060                                                const double& weightInMix)
0061     : thePos(GlobalPoint(0, 0, 0)),
0062       theTime(dNaN),
0063       theErr(AlgebraicSymMatrix44()),
0064       theWeight(posWeight),
0065       theWeightTimesPos(weightTimesPosition[0], weightTimesPosition[1], weightTimesPosition[2], 0),
0066       theWeightInMix(weightInMix),
0067       thePosAvailable(false),
0068       theTimeAvailable(false),
0069       theErrAvailable(false),
0070       theWeightAvailable(true),
0071       theWeightTimesPosAvailable(true),
0072       valid(true),
0073       vertexIs4D(false) {}
0074
0075 // no-offdiags for time
0076 BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
0077                                                const GlobalError& posErr,
0078                                                const double time,
0079                                                const double timeErr,
0080                                                const double& weightInMix)
0081     : thePos(pos),
0082       theTime(time),
0083       theErr(posErr),
0084       theWeight(AlgebraicSymMatrix44()),
0085       theWeightTimesPos(AlgebraicVector4()),
0086       theWeightInMix(weightInMix),
0087       thePosAvailable(true),
0088       theTimeAvailable(true),
0089       theErrAvailable(true),
0090       theWeightAvailable(false),
0091       theWeightTimesPosAvailable(false),
0092       valid(true),
0093       vertexIs4D(true) {
0094   // You dumb bastard. It's not a schooner, its a sailboat.
0095   GlobalError timeErrMat(0., 0., 0., 0., 0., 0., 0., 0., 0., timeErr * timeErr);
0096   theErr = theErr + timeErrMat;
0097 }
0098
0099 BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
0100                                                const GlobalWeight& posWeight,
0101                                                const double time,
0102                                                const double timeWeight,
0103                                                const double& weightInMix)
0104     : thePos(pos),
0105       theTime(time),
0106       theErr(AlgebraicSymMatrix44()),
0107       theWeight(posWeight),
0108       theWeightTimesPos(AlgebraicVector4()),
0109       theWeightInMix(weightInMix),
0110       thePosAvailable(true),
0111       theTimeAvailable(true),
0112       theErrAvailable(false),
0113       theWeightAvailable(true),
0114       theWeightTimesPosAvailable(false),
0115       valid(true),
0116       vertexIs4D(true) {
0117   GlobalWeight timeWeightMat(0., 0., 0., 0., 0., 0., 0., 0., 0., timeWeight);
0118   theWeight = theWeight + timeWeightMat;
0119 }
0120
0121 BasicSingleVertexState::BasicSingleVertexState(const AlgebraicVector3& weightTimesPosition,
0122                                                const GlobalWeight& posWeight,
0123                                                const double weightTimesTime,
0124                                                const double timeWeight,
0125                                                const double& weightInMix)
0126     : thePos(GlobalPoint(0, 0, 0)),
0127       theTime(dNaN),
0128       theErr(AlgebraicSymMatrix44()),
0129       theWeight(posWeight),
0130       theWeightTimesPos(weightTimesPosition[0], weightTimesPosition[1], weightTimesPosition[2], weightTimesTime),
0131       theWeightInMix(weightInMix),
0132       thePosAvailable(false),
0133       theTimeAvailable(false),
0134       theErrAvailable(false),
0135       theWeightAvailable(true),
0136       theWeightTimesPosAvailable(true),
0137       valid(true),
0138       vertexIs4D(true) {
0139   GlobalWeight timeWeightMat(0., 0., 0., 0., 0., 0., 0., 0., 0., timeWeight);
0140   theWeight = theWeight + timeWeightMat;
0141 }
0142
0143 // off-diags for time
0144 BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
0145                                                const double time,
0146                                                const GlobalError& posTimeErr,  // fully filled 4x4 matrix
0147                                                const double& weightInMix)
0148     : thePos(pos),
0149       theTime(time),
0150       theErr(posTimeErr),
0151       theWeight(AlgebraicSymMatrix44()),
0152       theWeightTimesPos(AlgebraicVector4()),
0153       theWeightInMix(weightInMix),
0154       thePosAvailable(true),
0155       theTimeAvailable(true),
0156       theErrAvailable(true),
0157       theWeightAvailable(false),
0158       theWeightTimesPosAvailable(false),
0159       valid(true),
0160       vertexIs4D(true) {}
0161
0162 BasicSingleVertexState::BasicSingleVertexState(const GlobalPoint& pos,
0163                                                const double time,
0164                                                const GlobalWeight& posTimeWeight,
0165                                                const double& weightInMix)
0166     : thePos(pos),
0167       theTime(time),
0168       theErr(AlgebraicSymMatrix44()),
0169       theWeight(posTimeWeight),
0170       theWeightTimesPos(AlgebraicVector4()),
0171       theWeightInMix(weightInMix),
0172       thePosAvailable(true),
0173       theTimeAvailable(true),
0174       theErrAvailable(false),
0175       theWeightAvailable(true),
0176       theWeightTimesPosAvailable(false),
0177       valid(true),
0178       vertexIs4D(true) {}
0179
0180 BasicSingleVertexState::BasicSingleVertexState(const AlgebraicVector4& weightTimesPosition,
0181                                                const GlobalWeight& posWeight,
0182                                                const double& weightInMix)
0183     : thePos(GlobalPoint(0, 0, 0)),
0184       theTime(dNaN),
0185       theErr(AlgebraicSymMatrix44()),
0186       theWeight(posWeight),
0187       theWeightTimesPos(weightTimesPosition),
0188       theWeightInMix(weightInMix),
0189       thePosAvailable(false),
0190       theTimeAvailable(false),
0191       theErrAvailable(false),
0192       theWeightAvailable(true),
0193       theWeightTimesPosAvailable(true),
0194       valid(true),
0195       vertexIs4D(true) {}
0196
0197 GlobalPoint BasicSingleVertexState::position() const {
0198   if (!valid)
0199     throw VertexException("BasicSingleVertexState::position::invalid");
0200   if (!thePosAvailable)
0201     computePosition();
0202   return thePos;
0203 }
0204
0205 GlobalError BasicSingleVertexState::error() const {
0206   if (!valid)
0207     throw VertexException("BasicSingleVertexState::error::invalid");
0208   if (!theErrAvailable)
0209     computeError();
0210   return GlobalError(theErr.matrix());
0211 }
0212
0213 GlobalError BasicSingleVertexState::error4D() const {
0214   if (!valid)
0215     throw VertexException("BasicSingleVertexState::error4D::invalid");
0216   if (!theErrAvailable)
0217     computeError();
0218   return theErr;
0219 }
0220
0221 double BasicSingleVertexState::time() const {
0222   if (!valid)
0223     throw VertexException("BasicSingleVertexState::time::invalid");
0224   if (!theTimeAvailable)
0225     computePosition();  // time computed with position (4-vector)
0226   return theTime;
0227 }
0228
0229 double BasicSingleVertexState::timeError() const {
0230   if (!valid)
0231     throw VertexException("BasicSingleVertexState::timeError::invalid");
0232   if (!theTimeAvailable)
0233     computeError();
0234   return std::sqrt(theErr.matrix4D()(3, 3));
0235 }
0236
0237 GlobalWeight BasicSingleVertexState::weight() const {
0238   if (!valid)
0239     throw VertexException("BasicSingleVertexState::weight::invalid");
0240   if (!theWeightAvailable)
0241     computeWeight();
0242   return GlobalWeight(theWeight.matrix());
0243 }
0244
0245 GlobalWeight BasicSingleVertexState::weight4D() const {
0246   if (!valid)
0247     throw VertexException("BasicSingleVertexState::weight4D::invalid");
0248   if (!theWeightAvailable)
0249     computeWeight();
0250   return theWeight;
0251 }
0252
0253 AlgebraicVector3 BasicSingleVertexState::weightTimesPosition() const {
0254   if (!valid)
0255     throw VertexException("BasicSingleVertexState::weightTimesPosition::invalid");
0256   if (!theWeightTimesPosAvailable)
0257     computeWeightTimesPos();
0258   return AlgebraicVector3(theWeightTimesPos[0], theWeightTimesPos[1], theWeightTimesPos[2]);
0259 }
0260
0261 AlgebraicVector4 BasicSingleVertexState::weightTimesPosition4D() const {
0262   if (!valid)
0263     throw VertexException("BasicSingleVertexState::weightTimesPosition4D::invalid");
0264   if (!theWeightTimesPosAvailable)
0265     computeWeightTimesPos();
0266   return theWeightTimesPos;
0267 }
0268
0269 double BasicSingleVertexState::weightInMixture() const {
0270   if (!valid)
0271     throw VertexException("BasicSingleVertexState::weightInMixture::invalid");
0272   return theWeightInMix;
0273 }
0274 // RefCountedVertexSeed BasicSingleVertexState::seedWithoutTracks() const
0275 // {
0276 //   RefCountedVertexSeed v = new VertexSeed(position(), error());
0277 //   return v;
0278 // }
0279
0280 void BasicSingleVertexState::computePosition() const {
0281   if (!valid)
0282     throw VertexException("BasicSingleVertexState::computePosition::invalid");
0283   AlgebraicVector4 pos = error4D().matrix4D() * weightTimesPosition4D();
0284   thePos = GlobalPoint(pos[0], pos[1], pos[2]);
0285   theTime = pos[3];
0286   thePosAvailable = true;
0287   theTimeAvailable = true;
0288 }
0289
0290 void BasicSingleVertexState::computeError() const {
0291   if (!valid)
0292     throw VertexException("BasicSingleVertexState::computeError::invalid");
0293   int ifail;
0294   if (vertexIs4D) {
0295     theErr = weight4D().matrix4D().Inverse(ifail);
0296     if (ifail != 0)
0297       throw VertexException("BasicSingleVertexState::could not invert weight matrix");
0298   } else {
0299     theErr = weight4D().matrix().Inverse(ifail);
0300     if (ifail != 0)
0301       throw VertexException("BasicSingleVertexState::could not invert weight matrix");
0302   }
0303   theErrAvailable = true;
0304 }
0305
0306 void BasicSingleVertexState::computeWeight() const {
0307   if (!valid)
0308     throw VertexException("BasicSingleVertexState::computeWeight::invalid");
0309   int ifail;
0310   if (vertexIs4D) {
0311     theWeight = error4D().matrix4D().Inverse(ifail);
0312     if (ifail != 0)
0313       throw VertexException("BasicSingleVertexState::could not invert error matrix");
0314   } else {
0315     theWeight = error4D().matrix().Inverse(ifail);
0316     if (ifail != 0)
0317       throw VertexException("BasicSingleVertexState::could not invert error matrix");
0318   }
0319   theWeightAvailable = true;
0320 }
0321
0322 void BasicSingleVertexState::computeWeightTimesPos() const {
0323   if (!valid)
0324     throw VertexException("BasicSingleVertexState::computeWeightTimesPos::invalid");
0325   AlgebraicVector4 pos;
0326   pos(0) = position().x();
0327   pos(1) = position().y();
0328   pos(2) = position().z();
0329   if (vertexIs4D) {
0330     pos(3) = theTime;
0331   } else {
0332     pos(3) = 0.;
0333   }
0334   theWeightTimesPos = weight4D().matrix4D() * pos;
0335   theWeightTimesPosAvailable = true;
0336 }