GsfVertexUpdator.cc

Go to the documentation of this file.

#include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexUpdator.h"
#include "RecoVertex/GaussianSumVertexFit/interface/BasicMultiVertexState.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <cfloat>

GsfVertexUpdator::GsfVertexUpdator(bool limit, const GsfVertexMerger * merger) :
  limitComponents (limit)
{
  if (limitComponents) theMerger = merger->clone();
}


CachingVertex<5> GsfVertexUpdator::add(const CachingVertex<5> & oldVertex,
                    const RefCountedVertexTrack track) const
{

  VSC prevVtxComponents = oldVertex.vertexState().components();

//   cout << "GsfVertexUpdator::Add new Track with "
//       << track->linearizedTrack()->components().size() << " components to vertex of "
//       << prevVtxComponents.size() << " components.\n";

  if (prevVtxComponents.empty()) {
  throw VertexException
    ("GsfVertexUpdator::(Previous) Vertex to update has no components");
  }

  LTC ltComponents = track->linearizedTrack()->components();
  if (ltComponents.empty()) {
  throw VertexException
    ("GsfVertexUpdator::Track to add to vertex has no components");
  }

  if ((ltComponents.size()==1) && (prevVtxComponents.size()==1)) 
    return kalmanVertexUpdator.add(oldVertex, track);

  float trackWeight = track->weight();

  std::vector<VertexComponent> newVertexComponents;
  newVertexComponents.reserve(prevVtxComponents.size()*ltComponents.size());

//     for (LTC::iterator trackCompIter = ltComponents.begin();
//      trackCompIter != ltComponents.end(); trackCompIter++ ) {
//   cout <<(**trackCompIter).state().globalPosition()<<endl;
//     }

  for (VSC::iterator vertexCompIter = prevVtxComponents.begin();
    vertexCompIter != prevVtxComponents.end(); vertexCompIter++ ) {
    for (LTC::iterator trackCompIter = ltComponents.begin();
    trackCompIter != ltComponents.end(); trackCompIter++ ) {
      newVertexComponents.push_back(
        createNewComponent(*vertexCompIter, *trackCompIter, trackWeight, +1));
     // return invalid vertex in case one of the updated vertex-components is invalid
      if (!newVertexComponents.back().first.isValid()) return CachingVertex<5>();
    }
  }
//   cout << "updator components: "<<newVertexComponents.size()<<endl;

  // Update tracks vector

  std::vector<RefCountedVertexTrack> newVertexTracks = oldVertex.tracks();
  newVertexTracks.push_back(track);
//   cout << "a \n ";

  // Assemble VertexStates and compute Chi**2

  VertexChi2Pair vertexChi2Pair = assembleVertexComponents(newVertexComponents);
//   cout << "b \n ";
  VertexState newVertexState = vertexChi2Pair.first;
//   cout << "c \n ";
  double chi2 = oldVertex.totalChiSquared() + vertexChi2Pair.second;

  // Merge:
  if (limitComponents) newVertexState = theMerger->merge(newVertexState);

  if  (oldVertex.hasPrior()) {
    return CachingVertex<5>(oldVertex.priorPosition(), oldVertex.priorError(),
            newVertexState.weightTimesPosition(),
        newVertexState.weight(), newVertexTracks, chi2);
  } else {
    return CachingVertex<5>(newVertexState, newVertexTracks, chi2);
  }


}




CachingVertex<5> GsfVertexUpdator::remove(const CachingVertex<5> & oldVertex, 
    const RefCountedVertexTrack track) const
{
  throw VertexException("GsfVertexUpdator::Remove Methode not yet done");
//  return CachingVertex<5>();
}


GsfVertexUpdator::VertexComponent 
GsfVertexUpdator::createNewComponent(const VertexState & oldVertex,
     const RefCountedLinearizedTrackState linTrack, float weight, int sign) const
{

  if(abs(sign) != 1)
    throw VertexException ("GsfVertexUpdator::sign not equal to 1.");

  if(sign == -1)
    throw VertexException("GsfVertexUpdator::sign of -1 not yet implemented.");


  // Weight of the component in the mixture (non-normalized)
  double weightInMixture = theWeightCalculator.calculate(oldVertex, linTrack, 1.E9);
  if (weightInMixture < 0.) return VertexComponent(VertexState(), WeightChi2Pair(0.,0.));

  // position estimate of the component
  VertexState newVertex = kalmanVertexUpdator.positionUpdate(oldVertex, 
                linTrack, weight, sign);
  if (!newVertex.isValid()) return VertexComponent(newVertex, WeightChi2Pair(0.,0.));

  //Chi**2 contribution of the component
  std::pair <bool, double> chi2P = kalmanVertexUpdator.chi2Increment(oldVertex, newVertex, 
                linTrack, weight);
  if (!chi2P.first) return VertexComponent(VertexState(), WeightChi2Pair(0.,0.));
//         cout << "Update: "<<oldVertex.position()<<" "<<newVertex.position()<<" "<<chi2
//       <<" "<<linTrack->weightInMixture()<<" "<<weightInMixture<<endl;

  return VertexComponent(newVertex, WeightChi2Pair(weightInMixture, chi2P.second));
}

GsfVertexUpdator::VertexChi2Pair GsfVertexUpdator::assembleVertexComponents(
    const std::vector<GsfVertexUpdator::VertexComponent> & newVertexComponents) const
{
  VSC vertexComponents;
  vertexComponents.reserve(newVertexComponents.size());
  
  //renormalize weights
// cout << "assemble "<<newVertexComponents.size()<<endl;
  double totalWeight = 0.;
  double totalChi2 = 0.;

  for (std::vector<VertexComponent>::const_iterator iter = newVertexComponents.begin();
    iter != newVertexComponents.end(); iter ++) {
    totalWeight += iter->second.first;
  }
// cout << "totalWeight "<<totalWeight<<endl;
  if (totalWeight<DBL_MIN) {
    edm::LogWarning("GsfVertexUpdator") << "Updated Vertex has total weight of 0. "
    <<"The track is probably very far away.";
    return VertexChi2Pair( VertexState(), 0.);
  }

  for (std::vector<VertexComponent>::const_iterator iter = newVertexComponents.begin();
    iter != newVertexComponents.end(); iter ++) {
    double weight = iter->second.first/totalWeight;
    if (iter->second.first>DBL_MIN) {
      vertexComponents.push_back(VertexState(iter->first.weightTimesPosition(),
       iter->first.weight(), weight));
      totalChi2 += iter->second.second * weight;
    }
  }
// cout << "totalChi2 "<<totalChi2<<endl;
// cout << "vertexComponents "<<vertexComponents.size()<<endl;

  if (vertexComponents.empty()){
    edm::LogWarning("GsfVertexUpdator") << "No Vertex State left after reweighting.";
    return VertexChi2Pair( VertexState(), 0.);
  }

  return VertexChi2Pair( VertexState( new BasicMultiVertexState( vertexComponents)),
            totalChi2);
}