CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SequentialVertexFitter.cc
Go to the documentation of this file.
4 
5 #include <algorithm>
6 using namespace std;
7 using namespace reco;
8 
9 namespace {
10  // FIXME
11  // hard-coded tracker bounds
12  // workaround while waiting for Geometry service
13  static const float TrackerBoundsRadius = 112;
14  static const float TrackerBoundsHalfLength = 273.5;
15  bool insideTrackerBounds(const GlobalPoint& point) {
16  return ((point.transverse() < TrackerBoundsRadius)
17  && (abs(point.z()) < TrackerBoundsHalfLength));
18  }
19 }
20 
21 
22 template <unsigned int N>
24  const LinearizationPointFinder & linP,
25  const VertexUpdator<N> & updator, const VertexSmoother<N> & smoother,
26  const AbstractLTSFactory<N> & ltsf ) :
27  theLinP(linP.clone()), theUpdator(updator.clone()),
28  theSmoother(smoother.clone()), theLTrackFactory ( ltsf.clone() )
29 {
31 }
32 
33 template <unsigned int N>
35  const edm::ParameterSet& pSet, const LinearizationPointFinder & linP,
36  const VertexUpdator<N> & updator, const VertexSmoother<N> & smoother,
37  const AbstractLTSFactory<N> & ltsf) :
38  thePSet(pSet), theLinP(linP.clone()), theUpdator(updator.clone()),
39  theSmoother(smoother.clone()), theLTrackFactory ( ltsf.clone() )
40 {
42 }
43 
44 
45 template <unsigned int N>
48 {
49  thePSet = original.parameterSet();
50  theLinP = original.linearizationPointFinder()->clone();
51  theUpdator = original.vertexUpdator()->clone();
52  theSmoother = original.vertexSmoother()->clone();
53  theMaxShift = original.maxShift();
54  theMaxStep = original.maxStep();
55  theLTrackFactory = original.linearizedTrackStateFactory()->clone();
56 }
57 
58 
59 template <unsigned int N>
61 {
62  delete theLinP;
63  delete theUpdator;
64  delete theSmoother;
65  delete theLTrackFactory;
66 }
67 
68 
69 template <unsigned int N>
71 {
72  theMaxShift = thePSet.getParameter<double>("maxDistance"); //0.01
73  theMaxStep = thePSet.getParameter<int>("maxNbrOfIterations"); //10
74 }
75 
76 template <unsigned int N>
78 {
79  thePSet.addParameter<double>("maxDistance", 0.01);
80  thePSet.addParameter<int>("maxNbrOfIterations", 10); //10
81  readParameters();
82 }
83 
84 template <unsigned int N>
86 SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack> & tracks) const
87 {
88  // Linearization Point
89  GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
90  if (!insideTrackerBounds(linP)) linP = GlobalPoint(0,0,0);
91 
92  // Initial vertex state, with a very large error matrix
93  ROOT::Math::SMatrixIdentity id;
94  AlgebraicSymMatrix33 we(id);
95  GlobalError error(we*10000);
96  VertexState state(linP, error);
97  std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
98  return fit(vtContainer, state, false);
99 }
100 
101 template <unsigned int N>
103  const std::vector<RefCountedVertexTrack> & tracks,
104  const reco::BeamSpot & spot ) const
105 {
106  VertexState state(spot);
107  return fit(tracks, state, true );
108 }
109 
110 template <unsigned int N>
112 SequentialVertexFitter<N>::vertex(const std::vector<RefCountedVertexTrack> & tracks) const
113 {
114  // Initial vertex state, with a very small weight matrix
115  GlobalPoint linP = tracks[0]->linearizedTrack()->linearizationPoint();
116  ROOT::Math::SMatrixIdentity id;
117  AlgebraicSymMatrix33 we(id);
118  GlobalError error(we*10000);
119  VertexState state(linP, error);
120  return fit(tracks, state, false);
121 }
122 
123 
124 // Fit vertex out of a set of RecTracks.
125 // Uses the specified linearization point.
126 //
127 template <unsigned int N>
129 SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack> & tracks,
130  const GlobalPoint& linPoint) const
131 {
132  // Initial vertex state, with a very large error matrix
133  ROOT::Math::SMatrixIdentity id;
134  AlgebraicSymMatrix33 we(id);
135  GlobalError error(we*10000);
136  VertexState state(linPoint, error);
137  std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
138  return fit(vtContainer, state, false);
139 }
140 
141 
146 template <unsigned int N>
148 SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack> & tracks,
149  const BeamSpot& beamSpot) const
150 {
151  VertexState beamSpotState(beamSpot);
152  std::vector<RefCountedVertexTrack> vtContainer;
153 
154  if (tracks.size() > 1) {
155  // Linearization Point search if there are more than 1 track
156  GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
157  if (!insideTrackerBounds(linP)) linP = GlobalPoint(0,0,0);
158  ROOT::Math::SMatrixIdentity id;
159  AlgebraicSymMatrix33 we(id);
160  GlobalError error(we*10000);
161  VertexState lpState(linP, error);
162  vtContainer = linearizeTracks(tracks, lpState);
163  } else {
164  // otherwise take the beamspot position.
165  vtContainer = linearizeTracks(tracks, beamSpotState);
166  }
167 
168  return fit(vtContainer, beamSpotState, true);
169 }
170 
171 
172 // Fit vertex out of a set of RecTracks.
173 // Uses the position as both the linearization point AND as prior
174 // estimate of the vertex position. The error is used for the
175 // weight of the prior estimate.
176 //
177 template <unsigned int N>
179  const std::vector<reco::TransientTrack> & tracks,
180  const GlobalPoint& priorPos,
181  const GlobalError& priorError) const
182 {
183  VertexState state(priorPos, priorError);
184  std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
185  return fit(vtContainer, state, true);
186 }
187 
188 // Fit vertex out of a set of VertexTracks
189 // Uses the position and error for the prior estimate of the vertex.
190 // This position is not used to relinearize the tracks.
191 //
192 template <unsigned int N>
194  const std::vector<RefCountedVertexTrack> & tracks,
195  const GlobalPoint& priorPos,
196  const GlobalError& priorError) const
197 {
198  VertexState state(priorPos, priorError);
199  return fit(tracks, state, true);
200 }
201 
202 
203 // Construct a container of VertexTrack from a set of RecTracks.
204 //
205 template <unsigned int N>
206 vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack>
208  const std::vector<reco::TransientTrack> & tracks,
209  const VertexState state) const
210 {
211  GlobalPoint linP = state.position();
212  std::vector<RefCountedVertexTrack> finalTracks;
213  finalTracks.reserve(tracks.size());
214  for(vector<reco::TransientTrack>::const_iterator i = tracks.begin();
215  i != tracks.end(); i++) {
217  = theLTrackFactory->linearizedTrackState(linP, *i);
218  RefCountedVertexTrack vTrData = theVTrackFactory.vertexTrack(lTrData,state);
219  finalTracks.push_back(vTrData);
220  }
221  return finalTracks;
222 }
223 
224 
225 // Construct new a container of VertexTrack with a new linearization point
226 // and vertex state, from an existing set of VertexTrack, from which only the
227 // recTracks will be used.
228 //
229 template <unsigned int N>
230 vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack>
232  const std::vector<RefCountedVertexTrack> & tracks,
233  const VertexState state) const
234 {
235 
236  GlobalPoint linP = state.position();
237  std::vector<RefCountedVertexTrack> finalTracks;
238  finalTracks.reserve(tracks.size());
239  for(typename std::vector<RefCountedVertexTrack>::const_iterator i = tracks.begin();
240  i != tracks.end(); i++) {
242  (**i).linearizedTrack()->stateWithNewLinearizationPoint(linP);
243  // RefCountedLinearizedTrackState lTrData =
244  // theLTrackFactory->linearizedTrackState(linP,
245  // (**i).linearizedTrack()->track());
246  RefCountedVertexTrack vTrData =
247  theVTrackFactory.vertexTrack(lTrData,state, (**i).weight() );
248  finalTracks.push_back(vTrData);
249  }
250  return finalTracks;
251 }
252 
253 
254 // The method where the vertex fit is actually done!
255 //
256 template <unsigned int N>
258 SequentialVertexFitter<N>::fit(const std::vector<RefCountedVertexTrack> & tracks,
259  const VertexState priorVertex, bool withPrior ) const
260 {
261  std::vector<RefCountedVertexTrack> initialTracks;
262  GlobalPoint priorVertexPosition = priorVertex.position();
263  GlobalError priorVertexError = priorVertex.error();
264 
265  CachingVertex<N> returnVertex(priorVertexPosition,priorVertexError,initialTracks,0);
266  if (withPrior) {
267  returnVertex = CachingVertex<N>(priorVertexPosition,priorVertexError,
268  priorVertexPosition,priorVertexError,initialTracks,0);
269  }
270  CachingVertex<N> initialVertex = returnVertex;
271  std::vector<RefCountedVertexTrack> globalVTracks = tracks;
272 
273  // main loop through all the VTracks
274  bool validVertex = true;
275  int step = 0;
276  GlobalPoint newPosition = priorVertexPosition;
277  GlobalPoint previousPosition;
278  do {
279  CachingVertex<N> fVertex = initialVertex;
280  // make new linearized and vertex tracks for the next iteration
281  if(step != 0) globalVTracks = reLinearizeTracks(tracks,
282  returnVertex.vertexState());
283 
284  // update sequentially the vertex estimate
285  for (typename std::vector<RefCountedVertexTrack>::const_iterator i
286  = globalVTracks.begin(); i != globalVTracks.end(); i++) {
287  fVertex = theUpdator->add(fVertex,*i);
288  if (!fVertex.isValid()) break;
289  }
290 
291  validVertex = fVertex.isValid();
292  // check tracker bounds and NaN in position
293  if (validVertex && hasNan(fVertex.position())) {
294  LogDebug("RecoVertex/SequentialVertexFitter")
295  << "Fitted position is NaN.\n";
296  validVertex = false;
297  }
298 
299  if (validVertex && !insideTrackerBounds(fVertex.position())) {
300  LogDebug("RecoVertex/SequentialVertexFitter")
301  << "Fitted position is out of tracker bounds.\n";
302  validVertex = false;
303  }
304 
305  if (!validVertex) {
306  // reset initial vertex position to (0,0,0) and force new iteration
307  // if number of steps not exceeded
308  ROOT::Math::SMatrixIdentity id;
309  AlgebraicSymMatrix33 we(id);
310  GlobalError error(we*10000);
311  fVertex = CachingVertex<N>(GlobalPoint(0,0,0), error,
312  initialTracks, 0);
313  }
314 
315  previousPosition = newPosition;
316  newPosition = fVertex.position();
317 
318  returnVertex = fVertex;
319  globalVTracks.clear();
320  step++;
321  } while ( (step != theMaxStep) &&
322  (((previousPosition - newPosition).transverse() > theMaxShift) ||
323  (!validVertex) ) );
324 
325  if (!validVertex) {
326  LogDebug("RecoVertex/SequentialVertexFitter")
327  << "Fitted position is invalid (out of tracker bounds or has NaN). Returned vertex is invalid\n";
328  return CachingVertex<N>(); // return invalid vertex
329  }
330 
331  if (step >= theMaxStep) {
332  LogDebug("RecoVertex/SequentialVertexFitter")
333  << "The maximum number of steps has been exceeded. Returned vertex is invalid\n";
334  return CachingVertex<N>(); // return invalid vertex
335  }
336 
337  // smoothing
338  returnVertex = theSmoother->smooth(returnVertex);
339 
340  return returnVertex;
341 }
342 
343 template class SequentialVertexFitter<5>;
344 template class SequentialVertexFitter<6>;
#define LogDebug(id)
int i
Definition: DBlmapReader.cc:9
std::vector< RefCountedVertexTrack > reLinearizeTracks(const std::vector< RefCountedVertexTrack > &tracks, const VertexState state) const
const VertexUpdator< N > * vertexUpdator() const
CachingVertex< N > fit(const std::vector< RefCountedVertexTrack > &tracks, const VertexState priorVertex, bool withPrior) const
VertexState const & vertexState() const
Definition: CachingVertex.h:85
list original
Definition: definitions.py:57
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
const AbstractLTSFactory< N > * linearizedTrackStateFactory() const
virtual CachingVertex< N > vertex(const std::vector< reco::TransientTrack > &tracks) const
T transverse() const
Definition: PV3DBase.h:73
const VertexSmoother< N > * vertexSmoother() const
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
T z() const
Definition: PV3DBase.h:64
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const edm::ParameterSet parameterSet() const
tuple tracks
Definition: testEve_cfg.py:39
TEveGeoShape * clone(const TEveElement *element, TEveElement *parent)
Definition: eve_macros.cc:135
GlobalPoint position() const
virtual LinearizationPointFinder * clone() const =0
T transverse() const
Another name for perp()
bool isValid() const
Definition: CachingVertex.h:96
std::vector< RefCountedVertexTrack > linearizeTracks(const std::vector< reco::TransientTrack > &tracks, const VertexState state) const
const float maxShift() const
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
const LinearizationPointFinder * linearizationPointFinder() const