GaussianSumVertexFit/src/GsfVertexMerger.cc

0001 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexMerger.h"
0002 #include "TrackingTools/GsfTools/interface/KullbackLeiblerDistance.h"
0003 // #include "TrackingTools/GsfTools/interface/MahalanobisDistance.h"
0004 #include "TrackingTools/GsfTools/interface/CloseComponentsMerger.h"
0005 // #include "CommonReco/GSFTools/interface/KeepingNonZeroWeightsMerger.h"
0006 // #include "TrackingTools/GsfTools/interface/LargestWeightsStateMerger.h"
0007 // #include "TrackingTools/GsfTools/interface/RCMultiGaussianState.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0010
0011 #include "RecoVertex/GaussianSumVertexFit/interface/VertexGaussianStateConversions.h"
0012
0013 GsfVertexMerger::GsfVertexMerger(const edm::ParameterSet& pSet) {
0014   maxComponents = pSet.getParameter<int>("maxNbrComponents");
0015   std::string mergerName = pSet.getParameter<std::string>("merger");
0016   std::string distanceName = pSet.getParameter<std::string>("distance");
0017
0018   if (mergerName == "CloseComponentsMerger") {
0019     DistanceBetweenComponents<3>* distance;
0020
0021     if (distanceName == "KullbackLeiblerDistance")
0022       distance = new KullbackLeiblerDistance<3>();
0023     //     else if ( distanceName=="MahalanobisDistance" )
0024     //       distance = new MahalanobisDistance();
0025     else
0026       throw VertexException("GsfVertexMerger: Distance type " + distanceName +
0027                             " unknown. Check distance parameter in GsfMergerParameters PSet");
0028
0029     merger = new CloseComponentsMerger<3>(maxComponents, distance);
0030     delete distance;
0031   }
0032   //   else if ( mergerName=="LargestWeightsStateMerger" )
0033   //     merger = new LargestWeightsStateMerger(maxComponents);
0034   else
0035     throw VertexException("GsfVertexMerger: Merger type " + mergerName +
0036                           " unknown. Check merger parameter in GsfMergerParameters PSet");
0037
0038   //   else if ( mergerName=="KeepingNonZeroWeightsMerger" )
0039   //     merger = new KeepingNonZeroWeightsMerger();
0040
0041   //   std::string mergerName = "test";
0042   //   edm::ESHandle<MultiGaussianStateMerger> mergerProducer;
0043   //   iRecord.get(mergerName,mergerProducer);
0044   //   merger = (MultiGaussianStateMerger *) mergerProducer.product();
0045 }
0046
0047 CachingVertex<5> GsfVertexMerger::merge(const CachingVertex<5>& oldVertex) const {
0048   if (oldVertex.vertexState().components().size() <= maxComponents)
0049     return oldVertex;
0050
0051   VertexState newVertexState = merge(oldVertex.vertexState());
0052
0053   if (oldVertex.hasPrior()) {
0054     return CachingVertex<5>(oldVertex.priorPosition(),
0055                             oldVertex.priorError(),
0056                             newVertexState.weightTimesPosition(),
0057                             newVertexState.weight(),
0058                             oldVertex.tracks(),
0059                             oldVertex.totalChiSquared());
0060   } else {
0061     return CachingVertex<5>(newVertexState, oldVertex.tracks(), oldVertex.totalChiSquared());
0062   }
0063 }
0064
0065 VertexState GsfVertexMerger::merge(const VertexState& oldVertex) const {
0066   using namespace GaussianStateConversions;
0067
0068   if (oldVertex.components().size() <= maxComponents)
0069     return oldVertex;
0070
0071   MultiGaussianState<3> multiGaussianState(multiGaussianStateFromVertex(oldVertex));
0072   MultiGaussianState<3> finalState(merger->merge(multiGaussianState));
0073   return vertexFromMultiGaussianState(finalState);
0074 }