diff --git a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetImprovedJetFitterVxFinder.h b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetImprovedJetFitterVxFinder.h
index 5ef9dccce8b08bf99c26e335adaa57edbb5ed2fa..2bfae6ffa44f16d40ac677b419b8e55a102cd8b0 100755
--- a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetImprovedJetFitterVxFinder.h
+++ b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetImprovedJetFitterVxFinder.h
@@ -2,6 +2,7 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
+
 /***************************************************************************
                           InDetImprovedJetFitterVxFinder.cxx  -  Description
                              -------------------
@@ -128,6 +129,18 @@ namespace InDet {
     double m_vertexProbCut;
     int m_maxClusteringIterations;
     double m_vertexClusteringProbabilityCut;
+    double m_vertexClusteringProbabilityCutWithMass;
+    double m_vertexClusteringProbabilityCutWithMass0010;
+    double m_vertexClusteringProbabilityCutWithMass1015;
+    double m_vertexClusteringProbabilityCutWithMass1520;
+    double m_vertexClusteringProbabilityCutWithMass2025;
+    double m_vertexClusteringProbabilityCutWithMass2530;
+    double m_vertexClusteringProbabilityCutWithMass3040;
+    double m_vertexClusteringProbabilityCutWithMass4050;
+    double m_vertexClusteringProbabilityCutWithMass5060;
+    double m_vertexClusteringProbabilityCutWithMass6070;
+    double m_vertexClusteringTwoVtxMassForProbCut;
+
     bool m_useFastClustering;
 
     double m_cutCompatibilityPrimaryVertexForPositiveLifetimeTracks;
@@ -173,6 +186,8 @@ namespace InDet {
     double m_cutIPD0SingleTrackForBSecondSelection;
     double m_cutIPZ0SingleTrackForBSecondSelection;
     double m_cutPtSingleTrackForBSecondSelection;
+    double m_cutIPD0SigBoxSingleTrackForBSecondSelection; //box cut for PU rejection in 1-t vtx
+    double m_cutIPZ0SigBoxSingleTrackForBSecondSelection; //box cut for PU rejection in 1-t vtx
     double m_cutCompatibilityPrimaryVertexSinglePositiveLifetimeTrackForBSecondSelection;
     double m_cutCompatibilityPrimaryVertexSingleNegativeLifetimeTrackForBSecondSelection;
 //don't use the primary vertex finder information to estimate the compatibility! 
diff --git a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetJetFitterUtils.h b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetJetFitterUtils.h
index 6c69f0253e370b3eb3edb00090b34e04a5441e14..b126a9c81580ac7946fa71138d6b9c9c43098488 100644
--- a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetJetFitterUtils.h
+++ b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/InDetSecVxFinderTool/InDetJetFitterUtils.h
@@ -26,6 +26,7 @@
 #include "GaudiKernel/ToolHandle.h"
 #include "CLHEP/Matrix/SymMatrix.h"
 #include "CLHEP/Matrix/Matrix.h"
+#include "CLHEP/Vector/LorentzVector.h"
 //#include "TrkParticleBase/LinkToTrackParticleBase.h"
 //#include "TrkParticleBase/TrackParticleBaseCollection.h"
 //#include "TrkParticleBase/TrackParticleBase.h"
@@ -43,6 +44,7 @@ namespace Trk {
   class LinkToTrackParticleBase;
   class ITrackLink;
   class TrackParticleBase;
+  class VxVertexOnJetAxis;
 }
 
 namespace InDet {
@@ -97,6 +99,9 @@ namespace InDet {
     
     std::pair<double,double> getD0andZ0IP(const Trk::TrackParameters & trackPerigee,
                                           const Trk::Vertex & vertexToExtrapolateTo) const;
+
+    std::pair<double,double> getD0andZ0IPSig(const Trk::TrackParameters & trackPerigee,
+					     const Trk::RecVertex & vertex) const;
     
     const Trk::LinkToTrackParticleBase* findNeutralTrackParticleBase(const std::vector<const Trk::LinkToTrackParticleBase*> &,
                                                                      const Trk::VxCandidate &) const;
@@ -119,7 +124,8 @@ namespace InDet {
     bool checkIfVxCandidateIsInVector(const xAOD::Vertex * vertexToCheck,
                                       const std::vector<const xAOD::Vertex*> & vectorOfCandidates) const;
 
-   
+    CLHEP::HepLorentzVector fourMomentumAtVertex(const Trk::VxVertexOnJetAxis &) const;
+    
     
 
   private:
diff --git a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetImprovedJetFitterVxFinder.cxx b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetImprovedJetFitterVxFinder.cxx
index 92df540fceb807c2c01e3a8f4a646059b62faf98..4134c66a50826228172dc93e692b70997cf1be5b 100755
--- a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetImprovedJetFitterVxFinder.cxx
+++ b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetImprovedJetFitterVxFinder.cxx
@@ -101,7 +101,18 @@ namespace InDet
     m_maxNumDeleteIterations(30),
     m_vertexProbCut(0.001),
     m_maxClusteringIterations(30),
-    m_vertexClusteringProbabilityCut(0.002),
+    m_vertexClusteringProbabilityCut(0.005),//GP 11/6/16: changed to default in job option file
+    m_vertexClusteringProbabilityCutWithMass(0.05),//GP 11/7/16: keep default to not active
+    m_vertexClusteringProbabilityCutWithMass0010(0.002),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass1015(0.002),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass1520(0.050),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass2025(0.100),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass2530(0.200),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass3040(0.500),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass4050(0.700),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass5060(0.900),//MB Mass dependent cut
+    m_vertexClusteringProbabilityCutWithMass6070(0.900),//MB Mass dependent cut
+    m_vertexClusteringTwoVtxMassForProbCut(2000),
     m_useFastClustering(false),
     m_cutCompatibilityPrimaryVertexForPositiveLifetimeTracks(1e-1),
     m_cutCompatibilityPrimaryVertexForNegativeLifetimeTracks(5e-2),
@@ -128,19 +139,21 @@ namespace InDet
     m_cutTwoTrkVtxVertexProbForBFirstSelectionFirstCriterium(0.05),
     m_cutTwoTrkVtxLifetimeSignificanceForBFirstSelectionSecondCriterium(1.5),
     m_cutTwoTrkVtxVertexProbForBFirstSelectionSecondCriterium(0.034),
-    m_firstBeam_min(29.35-1.5*0.68),
-    m_firstBeam_max(29.35+1.5*0.68),
-    m_secondBeam_min(34.18-1.5*0.65),
-    m_secondBeam_max(34.18+1.5*0.65),
-    m_firstLayer_min(57.5),
-    m_firstLayer_max(46.5),
-    m_secondLayer_min(70.5-1.5*2.1),
-    m_secondLayer_max(70.5+1.5*2.1),
+    m_firstBeam_min(23), //29.35-1.5*0.68),
+    m_firstBeam_max(25), //35+1.5*0.68),
+    m_secondBeam_min(23), //34.18-1.5*0.65),
+    m_secondBeam_max(25), //34.18+1.5*0.65),
+    m_firstLayer_min(34.0-2.5), //57.5)
+    m_firstLayer_max(34.0+2.5), //46.5),
+    m_secondLayer_min(51.5-3), //70.5-1.5*2.1),
+    m_secondLayer_max(51.5+3), //70.5+1.5*2.1),
     m_cutCompatibilityToPrimarySingleTrackForMatInteractions(1e-4),
     m_cutCompatibilityToPrimaryBothTracksForMatInteractions(1e-6),
-    m_cutIPD0SingleTrackForBSecondSelection(1.5),
-    m_cutIPZ0SingleTrackForBSecondSelection(3.),
+    m_cutIPD0SingleTrackForBSecondSelection(1.5), 
+    m_cutIPZ0SingleTrackForBSecondSelection(3.),  
     m_cutPtSingleTrackForBSecondSelection(750.),
+    m_cutIPD0SigBoxSingleTrackForBSecondSelection(2.), //0 effectively disables the cut
+    m_cutIPZ0SigBoxSingleTrackForBSecondSelection(5.),
     m_cutCompatibilityPrimaryVertexSinglePositiveLifetimeTrackForBSecondSelection(5e-2),
     m_cutCompatibilityPrimaryVertexSingleNegativeLifetimeTrackForBSecondSelection(1e-2),
     m_doNotUsePrimaryVertexCombatibilityInfo(false),
@@ -179,7 +192,19 @@ namespace InDet
     declareProperty("VertexProbCut",m_vertexProbCut);
     declareProperty("MaxClusteringIterations",m_maxClusteringIterations);
     declareProperty("VertexClusteringProbabilityCut",m_vertexClusteringProbabilityCut);
+    declareProperty("VertexClusteringProbabilityCutWithMass",m_vertexClusteringProbabilityCutWithMass);
+    declareProperty("vertexClusteringTwoVtxMassForProbCut",m_vertexClusteringTwoVtxMassForProbCut);
     declareProperty("UseFastClustering",m_useFastClustering);
+    
+    declareProperty("VertexClusteringProbabilityCutWithMass0010",m_vertexClusteringProbabilityCutWithMass0010);
+    declareProperty("VertexClusteringProbabilityCutWithMass1015",m_vertexClusteringProbabilityCutWithMass1015);
+    declareProperty("VertexClusteringProbabilityCutWithMass1520",m_vertexClusteringProbabilityCutWithMass1520);
+    declareProperty("VertexClusteringProbabilityCutWithMass2025",m_vertexClusteringProbabilityCutWithMass2025);
+    declareProperty("VertexClusteringProbabilityCutWithMass2530",m_vertexClusteringProbabilityCutWithMass2530);
+    declareProperty("VertexClusteringProbabilityCutWithMass3040",m_vertexClusteringProbabilityCutWithMass3040);
+    declareProperty("VertexClusteringProbabilityCutWithMass4050",m_vertexClusteringProbabilityCutWithMass4050);
+    declareProperty("VertexClusteringProbabilityCutWithMass5060",m_vertexClusteringProbabilityCutWithMass5060);
+    declareProperty("VertexClusteringProbabilityCutWithMass6070",m_vertexClusteringProbabilityCutWithMass6070);
 
     //Cuts which steer the finding in the seeding phase before JetFitter (yes, so many!!!)
     declareProperty("cutCompPrimaryVertexForPosLifetimeTracks",m_cutCompatibilityPrimaryVertexForPositiveLifetimeTracks);
@@ -203,6 +228,8 @@ namespace InDet
     declareProperty("cutIPD0BothTracksForBFirstSelection",m_cutIPD0BothTracksForBFirstSelection);
     declareProperty("cutIPZ0BothTracksForBFirstSelection",m_cutIPZ0BothTracksForBFirstSelection);
     declareProperty("cutPtBothTracksForBFirstSelection",m_cutPtBothTracksForBFirstSelection);
+    declareProperty("cutIPD0SigBoxSingleTrackForBSecondSelection",m_cutIPD0SigBoxSingleTrackForBSecondSelection);
+    declareProperty("cutIPZ0SigBoxSingleTrackForBSecondSelection",m_cutIPZ0SigBoxSingleTrackForBSecondSelection);
     declareProperty("cutTwoTrkVtxLifeSignForBFirstSelectCriteriumA",m_cutTwoTrkVtxLifetimeSignificanceForBFirstSelectionFirstCriterium);
     declareProperty("cutTwoTrkVtxVtxProbForBFirstSelectCriteriumA",m_cutTwoTrkVtxVertexProbForBFirstSelectionFirstCriterium);
     declareProperty("cutTwoTrkVtxLifeSignForBFirstSelectCriteriumB",m_cutTwoTrkVtxLifetimeSignificanceForBFirstSelectionSecondCriterium);
@@ -1843,10 +1870,13 @@ namespace InDet
        
        std::pair<double,double> track_IPd0z0=m_jetFitterUtils->getD0andZ0IP(*perigee,
                                                                             primaryVertexRecVertex);
+       std::pair<double,double> track_IPd0z0Sig=m_jetFitterUtils->getD0andZ0IPSig(*perigee,
+										  primaryVertexRecVertex);
 
        const double & IPd0=track_IPd0z0.first;
        const double & IPz0=track_IPd0z0.second;
-
+       const double & IPd0Sig=track_IPd0z0Sig.first;
+       const double & IPz0Sig=track_IPd0z0Sig.second;
        const double pT=perigee->momentum().perp();
        
        double cutCompatibilityPrimaryVertexSinglePositiveLifetimeTrackForBSecondSelection=m_cutCompatibilityPrimaryVertexSinglePositiveLifetimeTrackForBSecondSelection;
@@ -1858,12 +1888,13 @@ namespace InDet
          cutCompatibilityPrimaryVertexSingleNegativeLifetimeTrackForBSecondSelection=m_cutCompatibilityPrimaryVertexSinglePositiveLifetimeTrackForBSecondSelection;
        }
        
+       bool passBoxCut=( fabs(IPd0Sig)<m_cutIPD0SigBoxSingleTrackForBSecondSelection && fabs(IPz0Sig)>m_cutIPZ0SigBoxSingleTrackForBSecondSelection);
 
        if (fabs(IPd0)>m_cutIPD0SingleTrackForBSecondSelection||
            fabs(IPz0)>m_cutIPZ0SingleTrackForBSecondSelection||
            ( compatibilityOfActualTrack>=0 && TMath::Prob(fabs(compatibilityOfActualTrack),2)>cutCompatibilityPrimaryVertexSinglePositiveLifetimeTrackForBSecondSelection) ||
            ( compatibilityOfActualTrack<0 && TMath::Prob(fabs(compatibilityOfActualTrack),2)>cutCompatibilityPrimaryVertexSingleNegativeLifetimeTrackForBSecondSelection) ||
-           pT< m_cutPtSingleTrackForBSecondSelection) 
+           pT< m_cutPtSingleTrackForBSecondSelection || passBoxCut) 
        {
          if (msgLvl(MSG::DEBUG)) msg() << " Candidate didn't pass one of the selection cuts " << endmsg;
          continue;
@@ -2249,19 +2280,84 @@ namespace InDet
 	Trk::PairOfVxVertexOnJetAxis pairOfVxVertexOnJetAxis=clusteringTablePtr->getMostCompatibleVertices(probVertex);
 	//a PairOfVxVertexOnJetAxis is a std::pair<VxVertexOnJetAxis*,VxVertexOnJetAxis*>
 	
-	if (probVertex>0.&&probVertex>m_vertexClusteringProbabilityCut) {
+	float probVertexExcludingPrimary(0.);
+	Trk::PairOfVxVertexOnJetAxis pairOfVxVertexOnJetAxisExcludingPrimary=clusteringTablePtr->getMostCompatibleVerticesExcludingPrimary(probVertexExcludingPrimary);
+
+	bool firstProbIsWithPrimary= ( fabs(probVertex-probVertexExcludingPrimary)>1e-6 );
+
+	if (probVertex>0.&&probVertex>m_vertexClusteringProbabilityCut&&firstProbIsWithPrimary) {
 	  if (msgLvl(MSG::VERBOSE)) msg() <<" merging vtx number " << (*pairOfVxVertexOnJetAxis.first).getNumVertex() << 
-	    " and " << (*pairOfVxVertexOnJetAxis.second).getNumVertex() << endmsg;
+                                        " and " << (*pairOfVxVertexOnJetAxis.second).getNumVertex() << " (should be PV)." << endmsg;
           //	  const Trk::VxVertexOnJetAxis & mergedVertex=
           m_helper->mergeVerticesInJetCandidate(*pairOfVxVertexOnJetAxis.first,
                                                 *pairOfVxVertexOnJetAxis.second,
                                                 *myJetCandidate);
 	  //now you need to update the numbering scheme
 	  m_initializationHelper->updateTrackNumbering(myJetCandidate);//maybe this should be moved to a lower level...
+          continue;
+	}
 
-	} else {
-	  noMoreVerticesToCluster=true;
+        if (probVertexExcludingPrimary>0.)
+	{
+
+	  //GP suggested by Marco Battaglia, use vertex mass in order to decide wether to split or not, so derive vertex masses first
+	  const Trk::VxVertexOnJetAxis* firstVertex=pairOfVxVertexOnJetAxisExcludingPrimary.first;
+	  const Trk::VxVertexOnJetAxis* secondVertex=pairOfVxVertexOnJetAxisExcludingPrimary.second;
+          
+	  CLHEP::HepLorentzVector massVector1=m_jetFitterUtils->fourMomentumAtVertex(*firstVertex);//MeV
+	  CLHEP::HepLorentzVector massVector2=m_jetFitterUtils->fourMomentumAtVertex(*secondVertex);//MeV
+
+	  CLHEP::HepLorentzVector sumMassVector=massVector1+massVector2;
+
+	  double massTwoVertex=sumMassVector.mag();//MeV
+	  
+	  bool doMerge(false);
+
+	  double vertexClusteringProbabilityCutWithMass;
+
+	  if(massTwoVertex< 1000.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass0010;
+	  }else if(massTwoVertex< 1500.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass1015;
+	  }else if(massTwoVertex< 2000.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass1520;
+	  }else if(massTwoVertex< 2500.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass2025;
+	  }else if(massTwoVertex< 3000.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass2530;
+	  }else if(massTwoVertex< 4000.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass3040;
+	  }else if(massTwoVertex< 5000.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass4050;
+	  }else if(massTwoVertex< 6000.){
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass5060;
+	  }else{
+	    vertexClusteringProbabilityCutWithMass = m_vertexClusteringProbabilityCutWithMass6070;
+	  }
+
+	  if (probVertexExcludingPrimary>vertexClusteringProbabilityCutWithMass)
+          {
+            doMerge=true;
+          }
+
+	  if (doMerge)
+	  {
+
+	    if (msgLvl(MSG::VERBOSE)) msg() <<" merging vtx number " << (*pairOfVxVertexOnJetAxis.first).getNumVertex() <<
+					" and " << (*pairOfVxVertexOnJetAxis.second).getNumVertex() << " mass merged vertex: " << massTwoVertex << endmsg;
+	    
+	    m_helper->mergeVerticesInJetCandidate(*pairOfVxVertexOnJetAxisExcludingPrimary.first,
+						  *pairOfVxVertexOnJetAxisExcludingPrimary.second,
+						  *myJetCandidate);
+	    
+	    m_initializationHelper->updateTrackNumbering(myJetCandidate);//maybe this should be moved to a lower level...                                                                   
+	    continue;//go to next cycle, after a succesful merging
+	  }
 	}
+	  
+	noMoreVerticesToCluster=true;
+
+
       }
       numClusteringLoops+=1;
     } while (numClusteringLoops<m_maxClusteringIterations&&!(noMoreVerticesToCluster));
diff --git a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetJetFitterUtils.cxx b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetJetFitterUtils.cxx
index 9a37b025c24ad494644e7d70a178db529fe4f68c..a85736bd38884c26f5a56b8feadc633fa2e5c4f6 100644
--- a/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetJetFitterUtils.cxx
+++ b/InnerDetector/InDetRecTools/InDetSecVxFinderTool/src/InDetJetFitterUtils.cxx
@@ -23,6 +23,7 @@
 #include "VxVertex/VxCandidate.h"
 #include "VxVertex/ExtendedVxCandidate.h"
 #include "VxVertex/LinearizedTrack.h"
+#include "VxJetVertex/VxVertexOnJetAxis.h"
 #include "TrkSurfaces/PerigeeSurface.h"
 #include "TrkEventPrimitives/JacobianPxyzToPhiThetaQoverPspherical.h"
 #include "TrkEventPrimitives/ParamDefs.h"
@@ -518,7 +519,47 @@ namespace InDet
     return std::pair<double,double>(IPd0,IPz0);
   }
     
+
+  std::pair<double,double> InDetJetFitterUtils::getD0andZ0IPSig(const Trk::TrackParameters & trackPerigee,
+								const Trk::RecVertex & vertex) const
+  {
+    
+    if (m_linearizedTrackFactoryIsAvailable==false)
+    {
+      msg(MSG::ERROR) << "Cannot perform requested extrapolation. No extrapolator defined...Returning 0 compatibility..." << endreq;
+      return std::pair<double,double>(0,0);
+    }
+    
+    Trk::LinearizedTrack* myLinearizedTrack=m_LinearizedTrackFactory->linearizedTrack(&trackPerigee,vertex.position());
+    Amg::Vector3D vertexPosition;
+    vertexPosition[0]=vertex.position()[0];
+    vertexPosition[1]=vertex.position()[1];
+    vertexPosition[2]=vertex.position()[2];
+    
+    const AmgSymMatrix(5) & ExpectedCovariance=myLinearizedTrack->expectedCovarianceAtPCA();
+    AmgSymMatrix(2) weightReduced=ExpectedCovariance.block<2,2>(0,0);
+    //AmgSymMatrix(2) errorVertexReduced=vertex.covariancePosition().similarity(myLinearizedTrack->positionJacobian()).block<2,2>(0,0);
+    //weightReduced+=errorVertexReduced;
+
+    double IPd0Sig=(myLinearizedTrack->expectedParametersAtPCA())[0]/sqrt( weightReduced(0,0) );
+    double IPz0Sig=(myLinearizedTrack->expectedParametersAtPCA())[1]/sqrt( weightReduced(1,1) );
+   
+    /*
+    std::cout << " " << std::endl;
+    std::cout << " --> ExpectedCovariance : " << sqrt( weightReduced(0,0)) << " , " << sqrt( weightReduced(1,1)) << std::endl;
+    std::cout << " --> Covarianceposition : " << sqrt(vertex.covariancePosition()(0,0)) << " , " << sqrt(vertex.covariancePosition()(1,1)) << std::endl;
+    std::cout << " --> d0/z0              : " << myLinearizedTrack->expectedParametersAtPCA()[0] << " , " << myLinearizedTrack->expectedParametersAtPCA()[1] << std::endl;
+    std::cout << " --> d0Sig/z0Sig        : " << IPd0Sig << " , " << IPz0Sig << std::endl;
+    */    
+
+    delete myLinearizedTrack;
+    myLinearizedTrack=0;
+        
+    return std::pair<double,double>(IPd0Sig,IPz0Sig);
+  }
     
+
+
   const Trk::LinkToTrackParticleBase* InDetJetFitterUtils::findNeutralTrackParticleBase(const std::vector<const Trk::LinkToTrackParticleBase*> & /*neutralTracks*/,
                                                                                         const xAOD::Vertex & /*myVxCandidate*/) const 
   {
@@ -675,4 +716,34 @@ namespace InDet
     
   }
   
-}//end namespace InDet
+  CLHEP::HepLorentzVector  InDetJetFitterUtils::fourMomentumAtVertex(const Trk::VxVertexOnJetAxis & myVxVertexOnJetAxis) const
+  {
+    
+ 
+    const double s_pion=139.57018;
+    //hard coded pion mass
+ 
+    CLHEP::HepLorentzVector massVector(0,0,0,0);
+ 
+    const std::vector<Trk::VxTrackAtVertex*> & tracksOfVertex=myVxVertexOnJetAxis.getTracksAtVertex();
+    std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackBegin=tracksOfVertex.begin();
+    std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackEnd=tracksOfVertex.end();
+    for (std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackIter=clustersOfTrackBegin;
+         clustersOfTrackIter!=clustersOfTrackEnd;
+         ++clustersOfTrackIter)
+    {
+      if (dynamic_cast<const Trk::Perigee*>((*clustersOfTrackIter)->perigeeAtVertex())!=0)
+      {
+        
+        const Trk::TrackParameters* aMeasPer=(*clustersOfTrackIter)->perigeeAtVertex();
+        Amg::Vector3D mytrack(aMeasPer->momentum());
+        massVector+=CLHEP::HepLorentzVector(mytrack.x(),mytrack.y(),mytrack.z(),TMath::Sqrt(s_pion*s_pion+mytrack.mag()*mytrack.mag()));
+      }
+    }
+   
+    return massVector;
+  }
+  
+ 
+}
+//end namespace InDet