Page MenuHomeHEPForge
Files F11221790

SequentialVertexFitter.cc
No OneTemporary
Actions

SequentialVertexFitter.cc
View Options

#include "RecoVertex/VertexTools/interface/SequentialVertexFitter.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <algorithm>
using namespace std;
using namespace reco;
namespace {
// FIXME
// hard-coded tracker bounds
// workaround while waiting for Geometry service
static const float TrackerBoundsRadius = 112;
static const float TrackerBoundsHalfLength = 273.5;
bool insideTrackerBounds(const GlobalPoint& point) {
return ((point.transverse() < TrackerBoundsRadius)
&& (abs(point.z()) < TrackerBoundsHalfLength));
}
}
template <unsigned int N>
SequentialVertexFitter<N>::SequentialVertexFitter(
const LinearizationPointFinder & linP,
const VertexUpdator<N> & updator, const VertexSmoother<N> & smoother,
const AbstractLTSFactory<N> & ltsf ) :
theLinP(linP.clone()), theUpdator(updator.clone()),
theSmoother(smoother.clone()), theLTrackFactory ( ltsf.clone() )
{
setDefaultParameters();
}
template <unsigned int N>
SequentialVertexFitter<N>::SequentialVertexFitter(
const edm::ParameterSet& pSet, const LinearizationPointFinder & linP,
const VertexUpdator<N> & updator, const VertexSmoother<N> & smoother,
const AbstractLTSFactory<N> & ltsf) :
thePSet(pSet), theLinP(linP.clone()), theUpdator(updator.clone()),
theSmoother(smoother.clone()), theLTrackFactory ( ltsf.clone() )
{
readParameters();
}
template <unsigned int N>
SequentialVertexFitter<N>::SequentialVertexFitter(
const SequentialVertexFitter & original)
{
thePSet = original.parameterSet();
theLinP = original.linearizationPointFinder()->clone();
theUpdator = original.vertexUpdator()->clone();
theSmoother = original.vertexSmoother()->clone();
theMaxShift = original.maxShift();
theMaxStep = original.maxStep();
theLTrackFactory = original.linearizedTrackStateFactory()->clone();
}
template <unsigned int N>
SequentialVertexFitter<N>::~SequentialVertexFitter()
{
delete theLinP;
delete theUpdator;
delete theSmoother;
delete theLTrackFactory;
}
template <unsigned int N>
void SequentialVertexFitter<N>::readParameters()
{
theMaxShift = thePSet.getParameter<double>("maxDistance"); //0.01
theMaxStep = thePSet.getParameter<int>("maxNbrOfIterations"); //10
}
template <unsigned int N>
void SequentialVertexFitter<N>::setDefaultParameters()
{
thePSet.addParameter<double>("maxDistance", 0.01);
thePSet.addParameter<int>("maxNbrOfIterations", 10); //10
readParameters();
}
template <unsigned int N>
CachingVertex<N>
SequentialVertexFitter<N>::vertex(const vector<reco::TransientTrack> & tracks) const
{
// Linearization Point
GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
// cout << "[SequentialVertexFitter] called now w/o linpt" << endl;
if (!insideTrackerBounds(linP)) linP = GlobalPoint(0,0,0);
// Initial vertex state, with a very large error matrix
AlgebraicSymMatrix we(3,1);
GlobalError error(we*10000);
VertexState state(linP, error);
vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
return fit(vtContainer, state, false);
}
template <unsigned int N>
CachingVertex<N> SequentialVertexFitter<N>::vertex(
const vector<RefCountedVertexTrack> & tracks,
const reco::BeamSpot & spot ) const
{
VertexState state(spot);
return fit(tracks, state, true );
}
template <unsigned int N>
CachingVertex<N>
SequentialVertexFitter<N>::vertex(const vector<RefCountedVertexTrack> & tracks) const
{
// Initial vertex state, with a very small weight matrix
GlobalPoint linP = tracks[0]->linearizedTrack()->linearizationPoint();
AlgebraicSymMatrix we(3,1);
GlobalError error(we*10000);
VertexState state(linP, error);
return fit(tracks, state, false);
}
// Fit vertex out of a set of RecTracks.
// Uses the specified linearization point.
//
template <unsigned int N>
CachingVertex<N>
SequentialVertexFitter<N>::vertex(const vector<reco::TransientTrack> & tracks,
const GlobalPoint& linPoint) const
{
// Initial vertex state, with a very large error matrix
AlgebraicSymMatrix we(3,1);
GlobalError error(we*10000);
VertexState state(linPoint, error);
vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
/* cout << "[SequentialVertexFitter] hey, we use vertex with "
<< state.position() << endl; */
return fit(vtContainer, state, false);
}
/** Fit vertex out of a set of TransientTracks.
* The specified BeamSpot will be used as priot, but NOT for the linearization.
* The specified LinearizationPointFinder will be used to find the linearization point.
*/
template <unsigned int N>
CachingVertex<N>
SequentialVertexFitter<N>::vertex(const vector<reco::TransientTrack> & tracks,
const BeamSpot& beamSpot) const
{
VertexState beamSpotState(beamSpot);
vector<RefCountedVertexTrack> vtContainer;
if (tracks.size() > 1) {
// Linearization Point search if there are more than 1 track
GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
if (!insideTrackerBounds(linP)) linP = GlobalPoint(0,0,0);
AlgebraicSymMatrix we(3,1);
GlobalError error(we*10000);
VertexState lpState(linP, error);
vtContainer = linearizeTracks(tracks, lpState);
} else {
// otherwise take the beamspot position.
vtContainer = linearizeTracks(tracks, beamSpotState);
}
return fit(vtContainer, beamSpotState, true);
}
// Fit vertex out of a set of RecTracks.
// Uses the position as both the linearization point AND as prior
// estimate of the vertex position. The error is used for the
// weight of the prior estimate.
//
template <unsigned int N>
CachingVertex<N> SequentialVertexFitter<N>::vertex(
const vector<reco::TransientTrack> & tracks,
const GlobalPoint& priorPos,
const GlobalError& priorError) const
{
VertexState state(priorPos, priorError);
vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
return fit(vtContainer, state, true);
}
// Fit vertex out of a set of VertexTracks
// Uses the position and error for the prior estimate of the vertex.
// This position is not used to relinearize the tracks.
//
template <unsigned int N>
CachingVertex<N> SequentialVertexFitter<N>::vertex(
const vector<RefCountedVertexTrack> & tracks,
const GlobalPoint& priorPos,
const GlobalError& priorError) const
{
VertexState state(priorPos, priorError);
return fit(tracks, state, true);
}
// Construct a container of VertexTrack from a set of RecTracks.
//
template <unsigned int N>
vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack>
SequentialVertexFitter<N>::linearizeTracks(
const vector<reco::TransientTrack> & tracks,
const VertexState state) const
{
GlobalPoint linP = state.position();
vector<RefCountedVertexTrack> finalTracks;
finalTracks.reserve(tracks.size());
for(vector<reco::TransientTrack>::const_iterator i = tracks.begin();
i != tracks.end(); i++) {
RefCountedLinearizedTrackState lTrData
= theLTrackFactory->linearizedTrackState(linP, *i);
RefCountedVertexTrack vTrData = theVTrackFactory.vertexTrack(lTrData,state);
finalTracks.push_back(vTrData);
}
return finalTracks;
}
// Construct new a container of VertexTrack with a new linearization point
// and vertex state, from an existing set of VertexTrack, from which only the
// recTracks will be used.
//
template <unsigned int N>
vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack>
SequentialVertexFitter<N>::reLinearizeTracks(
const vector<RefCountedVertexTrack> & tracks,
const VertexState state) const
{
GlobalPoint linP = state.position();
vector<RefCountedVertexTrack> finalTracks;
finalTracks.reserve(tracks.size());
for(typename vector<RefCountedVertexTrack>::const_iterator i = tracks.begin();
i != tracks.end(); i++) {
RefCountedLinearizedTrackState lTrData =
(**i).linearizedTrack()->stateWithNewLinearizationPoint(linP);
// RefCountedLinearizedTrackState lTrData =
// theLTrackFactory->linearizedTrackState(linP,
// (**i).linearizedTrack()->track());
RefCountedVertexTrack vTrData =
theVTrackFactory.vertexTrack(lTrData,state, (**i).weight() );
finalTracks.push_back(vTrData);
}
return finalTracks;
}
// The method where the vertex fit is actually done!
//
template <unsigned int N>
CachingVertex<N>
SequentialVertexFitter<N>::fit(const vector<RefCountedVertexTrack> & tracks,
const VertexState priorVertex, bool withPrior ) const
{
// cout << "[SequentialVertexFitter] priorVertex in fit =" << priorVertex.position() << endl;
vector<RefCountedVertexTrack> initialTracks;
GlobalPoint priorVertexPosition = priorVertex.position();
GlobalError priorVertexError = priorVertex.error();
CachingVertex<N> returnVertex(priorVertexPosition,priorVertexError,initialTracks,0);
if (withPrior) {
returnVertex = CachingVertex<N>(priorVertexPosition,priorVertexError,
priorVertexPosition,priorVertexError,initialTracks,0);
}
CachingVertex<N> initialVertex = returnVertex;
vector<RefCountedVertexTrack> globalVTracks = tracks;
// main loop through all the VTracks
bool validVertex = true;
int step = 0;
GlobalPoint newPosition = priorVertexPosition;
GlobalPoint previousPosition;
do {
CachingVertex<N> fVertex = initialVertex;
// make new linearized and vertex tracks for the next iteration
if(step != 0) globalVTracks = reLinearizeTracks(tracks,
returnVertex.vertexState());
// update sequentially the vertex estimate
for (typename vector<RefCountedVertexTrack>::const_iterator i
= globalVTracks.begin(); i != globalVTracks.end(); i++) {
fVertex = theUpdator->add(fVertex,*i);
if (!fVertex.isValid()) break;
}
validVertex = fVertex.isValid();
// check tracker bounds and NaN in position
if (validVertex && hasNan(fVertex.position())) {
LogDebug("RecoVertex/SequentialVertexFitter")
<< "Fitted position is NaN.\n";
validVertex = false;
}
if (validVertex && !insideTrackerBounds(fVertex.position())) {
LogDebug("RecoVertex/SequentialVertexFitter")
<< "Fitted position is out of tracker bounds.\n";
validVertex = false;
}
if (!validVertex) {
// reset initial vertex position to (0,0,0) and force new iteration
// if number of steps not exceeded
GlobalError error(AlgebraicSymMatrix(3,1)*10000);
fVertex = CachingVertex<N>(GlobalPoint(0,0,0), error,
initialTracks, 0);
}
previousPosition = newPosition;
newPosition = fVertex.position();
returnVertex = fVertex;
globalVTracks.clear();
step++;
} while ( (step != theMaxStep) &&
(((previousPosition - newPosition).transverse() > theMaxShift) ||
(!validVertex) ) );
if (!validVertex) {
LogDebug("RecoVertex/SequentialVertexFitter")
<< "Fitted position is invalid (out of tracker bounds or has NaN). Returned vertex is invalid\n";
return CachingVertex<N>(); // return invalid vertex
}
if (step >= theMaxStep) {
LogDebug("RecoVertex/SequentialVertexFitter")
<< "The maximum number of steps has been exceeded. Returned vertex is invalid\n";
return CachingVertex<N>(); // return invalid vertex
}
// smoothing
returnVertex = theSmoother->smooth(returnVertex);
return returnVertex;
}
template class SequentialVertexFitter<5>;
template class SequentialVertexFitter<6>;

File Metadata

Mime Type
text/x-c++
Expires
Wed, May 14, 10:54 AM (20 h, 41 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111291
Default Alt Text
SequentialVertexFitter.cc (11 KB)

Event Timeline