GaussianSumVertexFit/src/MultiVertexStateCombiner.cc

0001 #include "RecoVertex/GaussianSumVertexFit/interface/MultiVertexStateCombiner.h"
0002 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0003
0004 VertexState MultiVertexStateCombiner::combine(const VSC& theMixture) const {
0005   if (theMixture.empty()) {
0006     throw VertexException("MultiVertexStateCombiner:: VertexState container to collapse is empty");
0007   }
0008
0009   if (theMixture.size() == 1) {
0010     //     #ifndef CMS_NO_COMPLEX_RETURNS
0011     return VertexState(theMixture.front().position(), theMixture.front().error(), 1.0);
0012     //     #else
0013     //       VertexState theFinalVM(theMixture.front().position(), theMixture.front().error(), 1.0);
0014     //       return theFinalVM;
0015     //     #endif
0016   }
0017
0018   AlgebraicVector3 meanPosition;
0019   double weightSum = 0.;
0020   AlgebraicSymMatrix33 measCovar1, measCovar2;
0021   for (VSC::const_iterator mixtureIter1 = theMixture.begin(); mixtureIter1 != theMixture.end(); mixtureIter1++) {
0022     double vtxWeight = mixtureIter1->weightInMixture();
0023
0024     GlobalPoint vertexPosition = mixtureIter1->position();
0025     AlgebraicVector3 vertexCoord1;
0026     vertexCoord1[0] = vertexPosition.x();
0027     vertexCoord1[1] = vertexPosition.y();
0028     vertexCoord1[2] = vertexPosition.z();
0029
0030     //    AlgebraicVector position = mixtureIter1->position().vector(); //???
0031     weightSum += vtxWeight;
0032     meanPosition += vtxWeight * vertexCoord1;
0033
0034     measCovar1 += vtxWeight * mixtureIter1->error().matrix();
0035     for (VSC::const_iterator mixtureIter2 = mixtureIter1 + 1; mixtureIter2 != theMixture.end(); mixtureIter2++) {
0036       GlobalPoint vertexPosition2 = mixtureIter2->position();
0037       AlgebraicVector3 vertexCoord2;
0038       vertexCoord2[0] = vertexPosition2.x();
0039       vertexCoord2[1] = vertexPosition2.y();
0040       vertexCoord2[2] = vertexPosition2.z();
0041       AlgebraicVector3 posDiff = vertexCoord1 - vertexCoord2;
0042       AlgebraicMatrix13 tmp;
0043       tmp.Place_in_row(posDiff, 0, 0);
0044       AlgebraicSymMatrix11 s = AlgebraicMatrixID();
0045       measCovar2 += vtxWeight * mixtureIter2->weightInMixture() * ROOT::Math::SimilarityT(tmp, s);
0046     }
0047   }
0048   meanPosition /= weightSum;
0049   AlgebraicSymMatrix33 measCovar = measCovar1 / weightSum + measCovar2 / weightSum / weightSum;
0050
0051   GlobalPoint newPos(meanPosition[0], meanPosition[1], meanPosition[2]);
0052
0053   // #ifndef CMS_NO_COMPLEX_RETURNS
0054   return VertexState(newPos, GlobalError(measCovar), weightSum);
0055   // #else
0056   //   VertexState theFinalVS(newPos, GlobalError(measCovar), weightSum);
0057   //   return theFinalVS;
0058   // #endif
0059 }