ATLAS Offline Software
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Private Attributes | List of all members
Trk::FullVertexFitter Class Reference

This class implements a full vertex fitting algorithm as proposed by P. More...

#include <FullVertexFitter.h>

Inheritance diagram for Trk::FullVertexFitter:
Collaboration diagram for Trk::FullVertexFitter:

Public Types

enum  FitError {
  FITOK , MATINV , NEGTRCHI2 , MAXCHI2 ,
  MAXTRCHI2 , NOTRKS , NOFIT
}

Public Member Functions

virtual StatusCode initialize () override
 FullVertexFitter (const std::string &t, const std::string &n, const IInterface *p)
virtual ~FullVertexFitter ()
 standard destructor
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const Amg::Vector3D &startingPoint) const override
 Interface for ParametersBase with starting point.
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const TrackParameters * > &perigeeList, const std::vector< const Trk::NeutralParameters * > &, const Amg::Vector3D &startingPoint) const override
 Interface for TrackParameters and NeutralParameters with starting point.
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const xAOD::Vertex &constraint) const override
 Interface for ParametersBase with vertex constraint.
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const TrackParameters * > &perigeeList, const std::vector< const Trk::NeutralParameters * > &, const xAOD::Vertex &constraint) const override
 Interface for TrackParameters and NeutralParameters with RecVertex starting point.
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList) const override
 Fit interface with no starting point or constraint.
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const TrackParameters * > &perigeeList, const std::vector< const Trk::NeutralParameters * > &) const override
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const Amg::Vector3D &startingPoint) const override
 *Interface for xAOD::TrackParticle with starting point
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const std::vector< const xAOD::NeutralParticle * > &, const Amg::Vector3D &startingPoint) const override
 *Interface for xAOD::TrackParticle and NeutralParticle with starting point
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const xAOD::Vertex &constraint) const override
virtual std::unique_ptr< xAOD::Vertex > fit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const std::vector< const xAOD::NeutralParticle * > &, const xAOD::Vertex &constraint) const override

Private Attributes

unsigned int m_maxIterations
double m_maxDchi2PerNdf
ToolHandle< Trk::IVertexLinearizedTrackFactory > m_linFactory
 Data members to store the results.

Detailed Description

This class implements a full vertex fitting algorithm as proposed by P.

Billoir. this algorithm tries to estimate the vertex position with refitting the track parameters.

Changes:

David Shope david.nosp@m..ric.nosp@m.hard..nosp@m.shop.nosp@m.e@cer.nosp@m.n.ch (2016-04-19)

EDM Migration to xAOD - from Trk::VxCandidate to xAOD::Vertex, from Trk::RecVertex to xAOD::Vertex, from Trk::Vertex to Amg::Vector3D

Definition at line 37 of file FullVertexFitter.h.

Member Enumeration Documentation

◆ FitError

enum Trk::FullVertexFitter::FitError
Enumerator
FITOK 
MATINV 
NEGTRCHI2 
MAXCHI2 
MAXTRCHI2 
NOTRKS 
NOFIT 

Definition at line 46 of file FullVertexFitter.h.

47 {
48 FITOK,
49 MATINV,
50 NEGTRCHI2,
51 MAXCHI2,
52 MAXTRCHI2,
53 NOTRKS,
54 NOFIT
55 };

Constructor & Destructor Documentation

◆ FullVertexFitter()

Trk::FullVertexFitter::FullVertexFitter ( const std::string & t,
const std::string & n,
const IInterface * p )

Definition at line 108 of file FullVertexFitter.cxx.

108 : base_class ( t,n,p ),
109 m_maxIterations ( 5 ),
110 m_maxDchi2PerNdf ( 0.000001 )
111
112 {
113 declareProperty ( "MaxIterations", m_maxIterations );
114 declareProperty ( "MaxChi2PerNdf", m_maxDchi2PerNdf );
115 declareProperty ( "LinearizedTrackFactory", m_linFactory );
116 declareInterface<IVertexFitter> ( this );
117 }
Trk::FullVertexFitter::m_linFactory
ToolHandle< Trk::IVertexLinearizedTrackFactory > m_linFactory
Data members to store the results.
Definition FullVertexFitter.h:180

◆ ~FullVertexFitter()

Trk::FullVertexFitter::~FullVertexFitter ( )
virtualdefault

standard destructor

Member Function Documentation

◆ fit() [1/10]

virtual std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const TrackParameters * > & perigeeList,
const std::vector< const Trk::NeutralParameters * > &  ) const
inlineoverridevirtual

Definition at line 115 of file FullVertexFitter.h.

120 {
121 msg(MSG::WARNING) << "FullVertexFitter::fit(fit(const std::vector<const "
122 "TrackParameters*>&,const std::vector<const "
123 "Trk::NeutralParameters*>&) ignoring neutrals"
124 << endmsg;
125 return fit(ctx, perigeeList);
126 };
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:63
Trk::FullVertexFitter::fit
virtual std::unique_ptr< xAOD::Vertex > fit(const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const Amg::Vector3D &startingPoint) const override
Interface for ParametersBase with starting point.
Definition FullVertexFitter.cxx:122
msg
MsgStream & msg
Definition testRead.cxx:32

◆ fit() [2/10]

virtual std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const TrackParameters * > & perigeeList,
const std::vector< const Trk::NeutralParameters * > & ,
const Amg::Vector3D & startingPoint ) const
inlineoverridevirtual

Interface for TrackParameters and NeutralParameters with starting point.

Definition at line 69 of file FullVertexFitter.h.

74 {
75 msg(MSG::WARNING)
76 << "FullVertexFitter::fit(fit(const std::vector<const "
77 "TrackParameters*>&,const std::vector<const "
78 "Trk::NeutralParameters*>&,const Amg::Vector3D&) ignoring neutrals"
79 << endmsg;
80 return fit(ctx, perigeeList, startingPoint);
81 };

◆ fit() [3/10]

virtual std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const TrackParameters * > & perigeeList,
const std::vector< const Trk::NeutralParameters * > & ,
const xAOD::Vertex & constraint ) const
inlineoverridevirtual

Interface for TrackParameters and NeutralParameters with RecVertex starting point.

Definition at line 92 of file FullVertexFitter.h.

97 {
98 msg(MSG::WARNING)
99 << "FullVertexFitter::fit(fit(const std::vector<const "
100 "TrackParameters*>&,const std::vector<const "
101 "Trk::NeutralParameters*>&,const xAOD::Vertex&) ignoring neutrals"
102 << endmsg;
103 return fit(ctx, perigeeList, constraint);
104 };

◆ fit() [4/10]

std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const Trk::TrackParameters * > & perigeeList ) const
overridevirtual

Fit interface with no starting point or constraint.

(0,0,0) will be assumed.

Definition at line 479 of file FullVertexFitter.cxx.

480 {
481 Amg::Vector3D tmpVtx(0.,0.,0.);
482 return fit ( ctx, perigeeList, tmpVtx );
483 }
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47

◆ fit() [5/10]

std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const Trk::TrackParameters * > & perigeeList,
const Amg::Vector3D & startingPoint ) const
overridevirtual

Interface for ParametersBase with starting point.

Definition at line 122 of file FullVertexFitter.cxx.

125 {
126 xAOD::Vertex constraint;
127 constraint.makePrivateStore();
128 constraint.setPosition( firstStartingPoint );
129 constraint.setCovariancePosition( AmgSymMatrix(3)(3,3) );
130 constraint.setFitQuality( 0.,0.);
131 return fit ( ctx, originalPerigees, constraint );
132 }
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition EventPrimitives.h:50
xAOD::Vertex_v1::setCovariancePosition
void setCovariancePosition(const AmgSymMatrix(3)&covariancePosition)
Sets the vertex covariance matrix.
xAOD::Vertex_v1::setPosition
void setPosition(const Amg::Vector3D &position)
Sets the 3-position.
Definition Vertex_v1.cxx:121
xAOD::Vertex_v1::setFitQuality
void setFitQuality(float chiSquared, float numberDoF)
Set the 'Fit Quality' information.
Definition Vertex_v1.cxx:172
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16

◆ fit() [6/10]

std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const Trk::TrackParameters * > & originalPerigees,
const xAOD::Vertex & firstStartingPoint ) const
overridevirtual

Interface for ParametersBase with vertex constraint.

the position of the constraint is ALWAYS the starting point

Definition at line 136 of file FullVertexFitter.cxx.

139 {
140 if ( originalPerigees.empty() )
141 {
142 ATH_MSG_VERBOSE("No tracks to fit in this event.");
143 return nullptr;
144 }
145
146 /* Initialisation of variables */
147 double chi2 = 2000000000000.;
148 unsigned int nRP = originalPerigees.size(); // Number of tracks to fit
149 int ndf = nRP * ( 5-3 ) - 3; // Number of degrees of freedom
150 if ( ndf < 0 ) {ndf = 1;}
151
152 /* Determine if we are doing a constraint fit.*/
153 bool constraint = false;
154 if ( firstStartingPoint.covariancePosition().trace() != 0. )
155 {
156 constraint = true;
157 ndf += 3;
158 ATH_MSG_DEBUG("Fitting with constraint.");
159 ATH_MSG_VERBOSE(firstStartingPoint.covariancePosition().inverse());
160 }
161
162 double chi2New=0.;
163
164 Amg::Vector3D linPoint ( firstStartingPoint.position() ); // linearization point for track parameters (updated for every iteration)
165
166 Amg::Vector3D p0;
167
168 std::vector<Amg::Vector3D> mom_at_Origin;
169
170 auto fittedVertex = std::make_unique<xAOD::Vertex>();
171 fittedVertex->makePrivateStore(); // xAOD::VertexContainer will take ownership of AuxStore when ActualVertex is added to it
172
173 std::vector<VxTrackAtVertex> tracksAtVertex;
174 std::vector<BilloirTrack> billoirTracks;
175
176 /* Iterate fits until the fit criteria are met, or the number of max
177 iterations is reached. */
178 for ( unsigned int niter=0; niter < m_maxIterations; ++niter )
179 {
180 ATH_MSG_VERBOSE("Start of iteration " << niter << ", starting point ("
181 << linPoint [0] << ", " << linPoint [1] << ", " << linPoint [2]
182 << ") and " << originalPerigees.size() << " tracks.");
183
184 billoirTracks.clear();
185 chi2New = 0.;
186
187 AmgMatrix(2,3) D;
188 BilloirVertex billoirVertex;
189
190 /* Linearize the track parameters wrt. starting point of the fit */
191 Amg::Vector3D globalPosition = linPoint;
192 Trk::PerigeeSurface perigeeSurface ( globalPosition );
193
194 /* Extrapolate the perigees to the startpoint of the fit */
195 unsigned int count(0);
196 for (const auto *originalPerigee : originalPerigees)
197 {
198
199 if ( niter == 0 )
200 {
201 // need to cast to access parameters() (this is to guarantee that the code knows if it is neutrals or charged parameters)
202 const Trk::TrackParameters* chargedParameters ( nullptr );
203 chargedParameters = dynamic_cast<const Trk::TrackParameters*> ( originalPerigee );
204 if ( chargedParameters==nullptr )
205 {
206 ATH_MSG_ERROR("Track parameters are not charged tracks ... full fit aborted (this will be handled correctly soon)");
207 return nullptr;
208 }
209
210
211 p0.x() = chargedParameters->parameters() [Trk::phi];
212 p0.y() = chargedParameters->parameters() [Trk::theta];
213 p0.z() = chargedParameters->parameters() [Trk::qOverP] ;
214 mom_at_Origin.push_back ( p0 );
215
216 }
217
218 // FOR NOW THE ONLY WAY TO CHECK IF EXTRAPOLATION WORKED IS TO CHECK THE RETURN POINTER
219 // if it is ZERO then take the original perigee
220 // there is a tiny tiny chance that the extrapolation fails because of a non converging runge kutta procedure
221 // in all other cases it did not extrapolate because the reference surface of the original perigee
222 // is already given to the extrapolation point (or very close nearby)
223
224 std::unique_ptr<LinearizedTrack> linTrack ( m_linFactory->linearizedTrack ( originalPerigee, linPoint ) );
225 if ( linTrack==nullptr )
226 {
227 ATH_MSG_DEBUG("Could not linearize track! Skipping this track!");
228 }
229 else
230 {
231 BilloirTrack locBilloirTrack;
232
233 locBilloirTrack.originalPerigee = originalPerigee;
234 double d0 = linTrack->expectedParametersAtPCA()[Trk::d0];
235 double z0 = linTrack->expectedParametersAtPCA()[Trk::z0];
236 double phi = linTrack->expectedParametersAtPCA()[Trk::phi];
237 double theta = linTrack->expectedParametersAtPCA()[Trk::theta];
238 double qOverP = linTrack->expectedParametersAtPCA()[Trk::qOverP];
239
240 // calculate f(V_0,p_0) f_d0 = f_z0 = 0
241 double f_phi = mom_at_Origin[count] ( 0 );
242 double f_theta = mom_at_Origin[count] ( 1 );
243 double f_qOverP = mom_at_Origin[count] ( 2 );
244
245 //calculate Dper[i]
246 locBilloirTrack.Dper[0] = d0;
247 locBilloirTrack.Dper[1] = z0;
248 locBilloirTrack.Dper[2] = phi - f_phi;
249 locBilloirTrack.Dper[3] = theta - f_theta;
250 locBilloirTrack.Dper[4] = qOverP -f_qOverP;
251
252 // the D matrix
253 AmgMatrix(5,3) D_mat;
254 D_mat.setZero();
255 D_mat = linTrack->positionJacobian();
256
257 // calculate the E matrix
258 AmgMatrix(5,3) E_mat;
259 E_mat.setZero();
260 E_mat = linTrack->momentumJacobian();
261
262 // and cache martix operations related to these matrices
263 AmgMatrix(3,5) Dt_W_mat;
264 Dt_W_mat.setZero();
265 AmgMatrix(3,5) Et_W_mat;
266 Et_W_mat.setZero();
267 Dt_W_mat = D_mat.transpose() * ( linTrack->expectedWeightAtPCA() );
268 Et_W_mat = E_mat.transpose() * ( linTrack->expectedWeightAtPCA() );
269
270 // these are the Billoir computed results
271 locBilloirTrack.Di_mat = D_mat;
272 locBilloirTrack.Ei_mat = E_mat;
273 locBilloirTrack.Gi_mat = Et_W_mat*E_mat;
274 locBilloirTrack.Bi_mat = Dt_W_mat * E_mat; // Di.T * Wi * Ei
275 locBilloirTrack.Ui_vec = Et_W_mat * locBilloirTrack.Dper; // Ei.T * Wi * dqi
276 locBilloirTrack.Ci_inv = Et_W_mat * E_mat ; // (Ei.T * Wi * Ei)^-1
277 // we need the inverse matrix here
278 locBilloirTrack.Ci_inv = locBilloirTrack.Ci_inv.inverse().eval();
279 // sum up over all tracks
280 billoirVertex.T_vec += Dt_W_mat * locBilloirTrack.Dper; // sum{Di.T * Wi * dqi}
281 billoirVertex.A_mat = billoirVertex.A_mat + Dt_W_mat * D_mat ; // sum{Di.T * Wi * Di}
282
283 // remember those results for all tracks
284 locBilloirTrack.BCi_mat = locBilloirTrack.Bi_mat * locBilloirTrack.Ci_inv; // BCi = Bi * Ci^-1
285
286 // and some summed results
287 billoirVertex.BCU_vec += locBilloirTrack.BCi_mat * locBilloirTrack.Ui_vec; // sum{Bi * Ci^-1 * Ui}
288 billoirVertex.BCB_mat = billoirVertex.BCB_mat + locBilloirTrack.BCi_mat * locBilloirTrack.Bi_mat.transpose() ;// sum{Bi * Ci^-1 * Bi.T}
289 locBilloirTrack.linTrack = std::move(linTrack);
290 billoirTracks.push_back ( std::move(locBilloirTrack) );
291 count++;
292
293 }
294 }
295
296 // calculate delta (billoirFrameOrigin-position), might be changed by the beam-const
297 Amg::Vector3D V_del = billoirVertex.T_vec - billoirVertex.BCU_vec; // V_del = T-sum{Bi*Ci^-1*Ui}
298 AmgMatrix(3,3) V_wgt_mat = billoirVertex.A_mat - billoirVertex.BCB_mat; // V_wgt = A-sum{Bi*Ci^-1*Bi.T}
299
300 if ( constraint )
301 {
302 // add V_del += wgtconst * (billoirFrameOrigin - Vconst) and V_wgt +=wgtconst
303 Amg::Vector3D constraintPosInBilloirFrame;
304 constraintPosInBilloirFrame.setZero();
305 // this will be 0 for first iteration but != 0 from second on
306 constraintPosInBilloirFrame[0] = firstStartingPoint.position() [0]-linPoint [0];
307 constraintPosInBilloirFrame[1] = firstStartingPoint.position() [1]-linPoint [1];
308 constraintPosInBilloirFrame[2] = firstStartingPoint.position() [2]-linPoint [2];
309
310 V_del += firstStartingPoint.covariancePosition().inverse() *constraintPosInBilloirFrame;
311 V_wgt_mat += firstStartingPoint.covariancePosition().inverse();
312 }
313 // cov(delta_V) = V_wgt^-1
314 AmgMatrix(3,3) cov_delta_V_mat = V_wgt_mat.inverse().eval() ;
315
316 // delta_V = cov_(delta_V) * V_del;
317 Amg::Vector3D delta_V = cov_delta_V_mat * V_del;
318 /* Momentum Calculation follows here. */
319 /* ---------------------------------------------------------------------------------------- */
320 std::vector<AmgMatrix(5,5)> cov_delta_P_mat ( nRP );
321 std::vector<double> chi2PerTrack;
322 chi2PerTrack.clear();
323
324 // ===> loop over tracks
325 std::vector<BilloirTrack>::iterator BTIter;
326 int iRP =0;
327 for ( BTIter = billoirTracks.begin(); BTIter != billoirTracks.end() ; ++BTIter )
328 {
329 // delta_P calculation
330 BilloirTrack locP ( ( *BTIter ) );
331 Amg::Vector3D delta_P = ( locP.Ci_inv ) * ( ( locP.Ui_vec )- ( locP.Bi_mat.transpose() ) *delta_V );
332 mom_at_Origin[iRP] ( 0 ) += delta_P[0];
333 if ( fabs ( mom_at_Origin[iRP] ( 0 ) ) > 10*M_PI )
334 { // protect while loop
335 if (msgLvl(MSG::WARNING))
336 {
337 msg(MSG::WARNING) << " Track direction angles have numerical problems, stop perigee parameter update." << endmsg;
338 msg(MSG::WARNING) << " Phi value: "<<mom_at_Origin[iRP] ( 0 ) <<endmsg;
339 }
340 return nullptr;
341 }
342 mom_at_Origin[iRP] ( 1 ) += delta_P[1];
343 if ( fabs ( mom_at_Origin[iRP] ( 0 ) ) > 5*M_PI )
344 { // protect while loop
345 if (msgLvl(MSG::WARNING))
346 {
347 msg(MSG::WARNING) << " Track direction angles have numerical problems, stop perigee parameter update." << endmsg;
348 msg(MSG::WARNING) << " Theta value: "<<mom_at_Origin[iRP] ( 1 ) <<endmsg;
349 }
350 return nullptr;
351 }
352 mom_at_Origin[iRP] ( 2 ) += delta_P[2];
353 // mom_at_Origin[iRP](0) -= M_PI > M_PI ? M_PI : 0;
354 // mom_at_Origin[iRP](0) += M_PI < -M_PI ? M_PI : 0;
355 //correct phi, theta coordinates
356 while ( fabs ( mom_at_Origin[iRP] ( 0 ) ) > M_PI ) mom_at_Origin[iRP] ( 0 ) += ( mom_at_Origin[iRP] ( 0 ) > 0 ) ? -2*M_PI : 2*M_PI;
357 while ( mom_at_Origin[iRP] ( 1 ) > 2*M_PI ) mom_at_Origin[iRP] ( 1 ) -= 2*M_PI;
358 while ( mom_at_Origin[iRP] ( 1 ) < -M_PI ) mom_at_Origin[iRP] ( 1 ) += M_PI;
359 if ( mom_at_Origin[iRP] ( 1 ) > M_PI )
360 {
361 mom_at_Origin[iRP] ( 1 ) = 2*M_PI - mom_at_Origin[iRP] ( 1 );
362 if ( mom_at_Origin[iRP] ( 0 ) >= 0 ) mom_at_Origin[iRP] ( 0 ) += ( mom_at_Origin[iRP] ( 0 ) >0 ) ? -M_PI : M_PI;
363 }
364 if ( mom_at_Origin[iRP] ( 1 ) < 0.0 )
365 {
366 mom_at_Origin[iRP] ( 1 ) = - mom_at_Origin[iRP] ( 1 );
367 if ( mom_at_Origin[iRP] ( 0 ) >= 0 ) mom_at_Origin[iRP] ( 0 ) += ( mom_at_Origin[iRP] ( 0 ) >0 ) ? -M_PI : M_PI;
368 }
369
370 //calculate 5x5 cov_delta_P matrix
371 // d(d0,z0,phi,theta,qOverP)/d(x,y,z,phi,theta,qOverP)-transformation matrix
372 AmgMatrix(5,6) trans_mat;
373 trans_mat.setZero();
374 trans_mat ( 0,0 ) = locP.Di_mat ( 0,0 ); trans_mat ( 0,1 ) = locP.Di_mat ( 0,1 );
375 trans_mat ( 1,0 ) = locP.Di_mat ( 1,0 ); trans_mat ( 1,1 ) = locP.Di_mat ( 1,1 ); trans_mat ( 1,2 ) = 1.;
376 trans_mat ( 2,3 ) = 1.; trans_mat ( 3,4 ) = 1.; trans_mat ( 4,5 ) = 1.;
377
378 //some intermediate calculations to get 5x5 matrix
379 //cov(V,V)
380 AmgMatrix(3,3) V_V_mat; V_V_mat.setZero(); V_V_mat = cov_delta_V_mat ;
381 // std::cout<<"V_V_mat = "<<V_V_mat<<std::endl;
382 //cov(V,P)
383 AmgMatrix(3,3) V_P_mat; V_P_mat.setZero(); V_P_mat = -cov_delta_V_mat*locP.Gi_mat*locP.Ci_inv ;
384 // std::cout<<"V_P_mat = "<<V_P_mat<<std::endl;
385 //cov(P,P)
386 AmgMatrix(3,3) P_P_mat; P_P_mat.setZero(); P_P_mat = locP.Ci_inv + locP.BCi_mat.transpose() *cov_delta_V_mat*locP.BCi_mat ;
387 // std::cout<<"P_P_mat = "<<P_P_mat<<std::endl;
388
389 AmgMatrix(6,6) cov_mat;
390 cov_mat.setZero();
391 cov_mat.block<3,3>(0,3) = V_P_mat;
392 cov_mat.block<3,3>(3,0) = V_P_mat.transpose() ;
393 cov_mat.block<3,3>(0,0) = V_V_mat ;
394 cov_mat.block<3,3>(3,3) = P_P_mat ;
395
396 // cov_delta_P calculation
397 // std::cout<<"trans_mat = "<<trans_mat<<std::endl;
398 cov_delta_P_mat[iRP] = trans_mat*cov_mat*trans_mat.transpose() ;
399 // std::cout<<"cov_delta_P_mat[iRP] = "<<cov_delta_P_mat[iRP]<<std::endl;
400
401 ++iRP;
402
403 /* Calculate chi2 per track. */
404 ( *BTIter ).chi2= ( ( ( *BTIter ).Dper- ( *BTIter ).Di_mat*delta_V - ( *BTIter ).Ei_mat*delta_P ).transpose() * ( *BTIter ).linTrack->expectedWeightAtPCA() * ( ( *BTIter ).Dper - ( *BTIter ).Di_mat*delta_V - ( *BTIter ).Ei_mat*delta_P ) ) [0];
405
406 if ( ( *BTIter ).chi2 < 0 )
407 {
408 ATH_MSG_WARNING( "FullVertexFitter::calculate: error in chi2_per_track" );
409 return nullptr;
410 }
411 chi2New += ( *BTIter ).chi2;
412 }
413
414 if ( constraint )
415 {
416 Amg::Vector3D deltaTrk;
417 deltaTrk.setZero();
418 // last term will also be 0 again but only in the first iteration
419 // = calc. vtx in billoir frame - ( constraint pos. in billoir frame )
420 deltaTrk[0] = delta_V[0] - ( firstStartingPoint.position() [0] - linPoint [0] );
421 deltaTrk[1] = delta_V[1] - ( firstStartingPoint.position() [1] - linPoint [1] );
422 deltaTrk[2] = delta_V[2] - ( firstStartingPoint.position() [2] - linPoint [2] );
423 chi2New += ( deltaTrk.transpose() * firstStartingPoint.covariancePosition().inverse() * deltaTrk ) [0];
424 }
425
426 /* assign new linearization point (= new vertex position in global frame) */
427 Amg::Vector3D tmpPos ( linPoint );
428 tmpPos[0] += delta_V[0]; tmpPos[1] += delta_V[1]; tmpPos[2] += delta_V[2];
429 linPoint = tmpPos;
430
431 if ( chi2New < chi2 )
432 {
433 /* Store the vertex */
434 chi2 = chi2New;
435 //const AmgMatrix(3,3) * newCovarianceMatrix = &cov_delta_V_mat ;
436 //const AmgMatrix(3,3) newErrorMatrix = newCovarianceMatrix->inverse().eval();
437 //fittedVertex = RecVertex ( linPoint.position(), newErrorMatrix, ndf, chi2 );
438
439 //the cov_delta_V_mat is already the inverted form. -katy 2/2/16
440 fittedVertex->setPosition( linPoint );
441 fittedVertex->setCovariancePosition( cov_delta_V_mat );
442 fittedVertex->setFitQuality( chi2, ndf );
443
444 // new go through vector and delete entries
445 /* // TODO: not needed anymore, tracksAtVertex doesn't store pointers - just the objects themselves <David Shope> (EDM Migration) 03/21/16
446 for ( std::vector<Trk::VxTrackAtVertex*>::const_iterator itr = tracksAtVertex.begin();
447 itr != tracksAtVertex.end(); ++itr )
448 {
449 delete ( *itr );
450 }
451 */
452
453 tracksAtVertex.clear();
454 /* Store the tracks at vertex */
455 Amg::Vector3D Vertex(linPoint[0], linPoint[1], linPoint[2]);
456 const Trk::PerigeeSurface Surface ( Vertex );
457 TrackParameters * refittedPerigee ( nullptr );
458 unsigned int iter= 0;
459 std::vector<BilloirTrack>::iterator BTIter;
460 for ( BTIter = billoirTracks.begin(); BTIter != billoirTracks.end() ; ++BTIter )
461 {
462 const AmgMatrix(5,5) * newTrackCovarianceMatrix = &cov_delta_P_mat[iter] ;
463 //Covariance matrix does not need to be inverted:
464 // AmgMatrix(5,5) newTrackErrorMatrix = (AmgMatrix(5,5)) newTrackCovarianceMatrix->inverse().eval();
465 AmgMatrix(5,5) newTrackErrorMatrix = (AmgMatrix(5,5)) newTrackCovarianceMatrix->eval();
466 refittedPerigee = new Trk::Perigee ( 0.,0.,mom_at_Origin[iter][0],mom_at_Origin[iter][1],mom_at_Origin[iter][2],
467 Surface, std::move(newTrackErrorMatrix) );
468 Trk::VxTrackAtVertex* tmpVxTrkAtVtx = new Trk::VxTrackAtVertex ( ( *BTIter ).chi2, refittedPerigee, ( *BTIter ).originalPerigee );
469 tracksAtVertex.push_back ( *tmpVxTrkAtVtx );
470 iter ++;
471 }
472 }
473 } // end of iteration
474 fittedVertex->vxTrackAtVertex() = tracksAtVertex;
475 //ATH_MSG_VERBOSE("Final Vertex Fitted: " << fittedVxCandidate->recVertex()); // TODO: can no longer print vertex after converting to xAOD
476 return fittedVertex;
477 }
M_PI
#define M_PI
Definition ActiveFraction.h:11
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
AmgMatrix
#define AmgMatrix(rows, cols)
Definition EventPrimitives.h:49
Vertex
boost::graph_traits< boost::adjacency_list< boost::vecS, boost::vecS, boost::bidirectionalS > >::vertex_descriptor Vertex
Definition FPGATrackSimGNNRoadMakerTool.h:35
if
if(pathvar)
Definition SealSharedLib.cxx:167
clear
void clear()
Empty the pool.
Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
xAOD::Vertex_v1::position
const Amg::Vector3D & position() const
Returns the 3-pos.
Definition Vertex_v1.cxx:106
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:525
count
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:148
TRTCalib_cfilter.p0
p0
Definition TRTCalib_cfilter.py:129
Trk::Perigee
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:33
Trk::theta
@ theta
Definition ParamDefs.h:66
Trk::qOverP
@ qOverP
perigee
Definition ParamDefs.h:67
Trk::phi
@ phi
Definition ParamDefs.h:75
Trk::d0
@ d0
Definition ParamDefs.h:63
Trk::z0
@ z0
Definition ParamDefs.h:64
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27

◆ fit() [7/10]

std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const xAOD::TrackParticle * > & vectorTrk,
const Amg::Vector3D & startingPoint ) const
overridevirtual

*Interface for xAOD::TrackParticle with starting point

Definition at line 488 of file FullVertexFitter.cxx.

489 {
490 xAOD::Vertex constraint;
491 constraint.makePrivateStore();
492 constraint.setPosition( startingPoint );
493 constraint.setCovariancePosition( AmgSymMatrix(3)(3,3) );
494 constraint.setFitQuality( 0.,0.);
495 return fit(ctx, vectorTrk, constraint);
496 }//end of the xAOD starting point fit method

◆ fit() [8/10]

virtual std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const xAOD::TrackParticle * > & vectorTrk,
const std::vector< const xAOD::NeutralParticle * > & ,
const Amg::Vector3D & startingPoint ) const
inlineoverridevirtual

*Interface for xAOD::TrackParticle and NeutralParticle with starting point

Definition at line 136 of file FullVertexFitter.h.

141 {
142 msg(MSG::WARNING)
143 << "FullVertexFitter::fit(fit(const std::vector<const "
144 "TrackParticle*>&,const std::vector<const "
145 "Trk::NeutralParticle*>&,const Amg::Vector3D&) ignoring neutrals"
146 << endmsg;
147 return fit(ctx, vectorTrk, startingPoint);
148 };

◆ fit() [9/10]

virtual std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const xAOD::TrackParticle * > & vectorTrk,
const std::vector< const xAOD::NeutralParticle * > & ,
const xAOD::Vertex & constraint ) const
inlineoverridevirtual

Definition at line 160 of file FullVertexFitter.h.

165 {
166 msg(MSG::WARNING)
167 << "FullVertexFitter::fit(fit(const std::vector<const "
168 "TrackParticle*>&,const std::vector<const "
169 "Trk::NeutralParticle*>&,const xAOD::Vertex&) ignoring neutrals"
170 << endmsg;
171 return fit(ctx, vectorTrk, constraint);
172 };

◆ fit() [10/10]

std::unique_ptr< xAOD::Vertex > Trk::FullVertexFitter::fit ( const EventContext & ctx,
const std::vector< const xAOD::TrackParticle * > & vectorTrk,
const xAOD::Vertex & constraint ) const
overridevirtual

Definition at line 499 of file FullVertexFitter.cxx.

500 {
501 if(vectorTrk.empty())
502 {
503 msg(MSG::INFO)<<"Empty vector of tracks passed"<<endmsg;
504 return nullptr;
505 }
506
507 //making a list of perigee out of the vector of tracks
508 std::vector<const Trk::TrackParameters*> measuredPerigees;
509
510 for(const auto *i : vectorTrk)
511 {
512 const Trk::TrackParameters * tmpMeasPer = &(i->perigeeParameters());
513
514 if(tmpMeasPer!=nullptr) measuredPerigees.push_back(tmpMeasPer);
515 else msg(MSG::INFO)<<"Failed to dynamic_cast this track parameters to perigee"<<endmsg; //TODO: Failed to implicit cast the perigee parameters to track parameters?
516 }
517
518
519 std::unique_ptr<xAOD::Vertex> fittedVertex = fit( ctx, measuredPerigees, constraint );
520
521 //assigning the input tracks to the fitted vertex through VxTrackAtVertices
522 {
523 if( fittedVertex->vxTrackAtVertexAvailable() ) // TODO: I don't think vxTrackAtVertexAvailable() does the same thing as a null pointer check!
524 {
525 if(!fittedVertex->vxTrackAtVertex().empty())
526 {
527 for(unsigned int i = 0; i <vectorTrk.size(); ++i)
528 {
529
530 LinkToXAODTrackParticle * linkTT = new LinkToXAODTrackParticle;
531 linkTT->setElement(vectorTrk[i]);
532
533 // vxtrackatvertex takes ownership!
534 ( fittedVertex->vxTrackAtVertex() )[i].setOrigTrack(linkTT);
535 }//end of loop for setting orig tracks in.
536 }//end of protection against unsuccessfull updates (no tracks were added)
537 }//end of vector of tracks check
538 }//end of pointer check
539
540 //now set links to xAOD::TrackParticles directly in the xAOD::Vertex
541 unsigned int VTAVsize = fittedVertex->vxTrackAtVertex().size();
542 for (unsigned int i = 0 ; i < VTAVsize ; ++i)
543 {
544 Trk::VxTrackAtVertex* VTAV = &( fittedVertex->vxTrackAtVertex().at(i) );
545 //TODO: Will this pointer really hold 0 if no VxTrackAtVertex is found?
546 if (not VTAV){
547 ATH_MSG_WARNING (" Trying to set link to xAOD::TrackParticle. The VxTrackAtVertex is not found");
548 continue;
549 }
550
551 Trk::ITrackLink* trklink = VTAV->trackOrParticleLink();
552
553 // See if the trklink is to an xAOD::TrackParticle
554 Trk::LinkToXAODTrackParticle* linkToXAODTP = dynamic_cast<Trk::LinkToXAODTrackParticle*>(trklink);
555 if (linkToXAODTP)
556 {
557
558 //Now set the new link to the xAOD vertex
559 fittedVertex->addTrackAtVertex(*linkToXAODTP, VTAV->weight());
560
561 } else {
562 ATH_MSG_WARNING ("Skipping track. Trying to set link to something else than xAOD::TrackParticle. Neutrals not supported.");
563 }
564 } //end of loop
565
566 return fittedVertex;
567
568 }//end of the xAOD constrained fit method
CxxUtils::at
mapped_type at(key_type key) const
Look up an element in the map.

◆ initialize()

StatusCode Trk::FullVertexFitter::initialize ( )
overridevirtual

Definition at line 91 of file FullVertexFitter.cxx.

92 {
93
94
95 if ( m_linFactory.retrieve().isFailure() )
96 {
97 msg(MSG::FATAL)<< "Failed to retrieve tool " << m_linFactory << endmsg;
98 return StatusCode::FAILURE;
99 }
100
101
102 msg(MSG::INFO) << "Retrieved tool " << m_linFactory << endmsg;
103
104
105 return StatusCode::SUCCESS;
106 }

Member Data Documentation

◆ m_linFactory

ToolHandle<Trk::IVertexLinearizedTrackFactory> Trk::FullVertexFitter::m_linFactory
private

Data members to store the results.

Definition at line 180 of file FullVertexFitter.h.

◆ m_maxDchi2PerNdf

double Trk::FullVertexFitter::m_maxDchi2PerNdf
private

Definition at line 176 of file FullVertexFitter.h.

◆ m_maxIterations

unsigned int Trk::FullVertexFitter::m_maxIterations
private

Definition at line 175 of file FullVertexFitter.h.

The documentation for this class was generated from the following files:
Generated on for ATLAS Offline Software by doxygen 1.16.1