ATLAS Offline Software
Loading...
Searching...
No Matches
Trk::TrkV0VertexFitter Class Reference

This class implements a vertex fitting algorithm optimised for V0 finding. More...

#include <TrkV0VertexFitter.h>

Inheritance diagram for Trk::TrkV0VertexFitter:
Collaboration diagram for Trk::TrkV0VertexFitter:

Public Types

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

Public Member Functions

virtual StatusCode initialize () override
virtual StatusCode finalize () override
 TrkV0VertexFitter (const std::string &t, const std::string &n, const IInterface *p)
virtual ~TrkV0VertexFitter ()
 standard destructor
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const Amg::Vector3D &startingPoint) const override
 Interface for xAOD::TrackParticle with Amg::Vector3D starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const std::vector< const xAOD::NeutralParticle * > &, const Amg::Vector3D &startingPoint) const override
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const xAOD::Vertex &constraint) const override
 Interface for xAOD::TrackParticle with xAOD::Vertex starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const std::vector< const xAOD::NeutralParticle * > &, const xAOD::Vertex &constraint) const override
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk) const
 Fit interface for xAOD::TrackParticle with no starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const Amg::Vector3D &startingPoint) const override
 Interface for Trk::TrackParameters with Amg::Vector3D starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const std::vector< const Trk::NeutralParameters * > &, const Amg::Vector3D &startingPoint) const override
 Interface for Trk::TrackParameters and NeutralParameters with Amg::Vector3D starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const xAOD::Vertex &constraint) const override
 Interface for TrackParameters with xAOD::Vertex starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const std::vector< const Trk::NeutralParameters * > &, const xAOD::Vertex &constraint) const override
 Interface for TrackParameters and NeutralParameters with xAOD::Vertex starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList) const override
 Fit interface for TrackParameters with no starting point.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const std::vector< const Trk::NeutralParameters * > &) const override
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const xAOD::TrackParticle * > &vectorTrk, const std::vector< double > &masses, const double &constraintMass, const xAOD::Vertex *pointingVertex, const Amg::Vector3D &startingPoint) const
 Methods specific for the V0 fitter Method taking a vector of tracks, vector of masses and starting point as arguments.
virtual std::unique_ptr< xAOD::Vertexfit (const EventContext &ctx, const std::vector< const Trk::TrackParameters * > &perigeeList, const std::vector< double > &masses, const double &constraintMass, const xAOD::Vertex *pointingVertex, const Amg::Vector3D &startingPoint) const
 Interface for Trk::TrackParameters with mass and pointing constraints.

Private Attributes

int m_maxIterations
double m_maxDchi2PerNdf
double m_maxR
double m_maxZ
bool m_firstMeas
bool m_deltaR
ToolHandle< Trk::IExtrapolatorm_extrapolator
 Data members to store the results.
SG::ReadCondHandleKey< AtlasFieldCacheCondObjm_fieldCacheCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"}

Detailed Description

This class implements a vertex fitting algorithm optimised for V0 finding.

The algorithm fits the full track information, incorporating the possibility of performing full kinematic constraints at the same time. The full covariance matrix from the fit, including track-track and track-vertex correlations is calculated and returned

Definition at line 39 of file TrkV0VertexFitter.h.

Member Enumeration Documentation

◆ FitError

Constructor & Destructor Documentation

◆ TrkV0VertexFitter()

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

Definition at line 38 of file TrkV0VertexFitter.cxx.

38 : base_class(t,n,p),
41 m_maxR(2000.),
42 m_maxZ(5000.),
43 m_firstMeas(true),
44 m_deltaR(false),
45 m_extrapolator("Trk::Extrapolator/InDetExtrapolator", this)
46 {
47 declareProperty("MaxIterations", m_maxIterations);
48 declareProperty("MaxChi2PerNdf", m_maxDchi2PerNdf);
49 declareProperty("MaxR", m_maxR);
50 declareProperty("MaxZ", m_maxZ);
51 declareProperty("FirstMeasuredPoint", m_firstMeas);
52 declareProperty("Use_deltaR", m_deltaR);
53 declareProperty("Extrapolator", m_extrapolator);
54 declareInterface<IVertexFitter>(this);
55 }
ToolHandle< Trk::IExtrapolator > m_extrapolator
Data members to store the results.

◆ ~TrkV0VertexFitter()

Trk::TrkV0VertexFitter::~TrkV0VertexFitter ( )
virtualdefault

standard destructor

Member Function Documentation

◆ finalize()

StatusCode Trk::TrkV0VertexFitter::finalize ( )
overridevirtual

Definition at line 74 of file TrkV0VertexFitter.cxx.

75 {
76 ATH_MSG_DEBUG( "Finalize successful" );
77 return StatusCode::SUCCESS;
78 }
#define ATH_MSG_DEBUG(x)

◆ fit() [1/13]

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

Fit interface for TrackParameters with no starting point.

Interface for Trk::TrackParameters with no starting point.

(0,0,0) will be assumed.

(0,0,0) will be assumed

Definition at line 223 of file TrkV0VertexFitter.cxx.

225 {
226 Amg::Vector3D tmpVtx;
227 tmpVtx.setZero();
228 return fit(ctx, originalPerigees, tmpVtx);
229 }
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 Amg::Vector3D starting point.
Eigen::Matrix< double, 3, 1 > Vector3D

◆ fit() [2/13]

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

Interface for Trk::TrackParameters with Amg::Vector3D starting point.

Definition at line 200 of file TrkV0VertexFitter.cxx.

203 {
204 std::vector<double> masses;
205 double constraintMass = -9999.;
206 xAOD::Vertex * pointingVertex = nullptr;
207 return fit(ctx, originalPerigees, masses, constraintMass, pointingVertex, firstStartingPoint);
208 }
Vertex_v1 Vertex
Define the latest version of the vertex class.

◆ fit() [3/13]

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

Definition at line 145 of file TrkV0VertexFitter.h.

150 {
151 msg(MSG::WARNING) << "TrkV0VertexFitter::fit(fit(const std::vector<const "
152 "Trk::TrackParameters*>&,const std::vector<const "
153 "Trk::NeutralParameters*>&) ignoring neutrals"
154 << endmsg;
155 return fit(ctx, perigeeList);
156 };
#define endmsg
MsgStream & msg
Definition testRead.cxx:32

◆ fit() [4/13]

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

Interface for Trk::TrackParameters and NeutralParameters with Amg::Vector3D starting point.

Definition at line 103 of file TrkV0VertexFitter.h.

108 {
109 msg(MSG::WARNING)
110 << "TrkV0VertexFitter::fit(fit(const std::vector<const "
111 "Trk::TrackParameters*>&,const std::vector<const "
112 "Trk::NeutralParameters*>&,const Amg::Vector3D&) ignoring neutrals"
113 << endmsg;
114 return fit(ctx, perigeeList, startingPoint);
115 };

◆ fit() [5/13]

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

Interface for TrackParameters and NeutralParameters with xAOD::Vertex starting point.

Definition at line 125 of file TrkV0VertexFitter.h.

130 {
131 msg(MSG::WARNING)
132 << "TrkV0VertexFitter::fit(fit(const std::vector<const "
133 "Trk::TrackParameters*>&,const std::vector<const "
134 "Trk::NeutralParameters*>&,const xAOD::Vertex&) ignoring neutrals"
135 << endmsg;
136 return fit(ctx, perigeeList, constraint);
137 };

◆ fit() [6/13]

std::unique_ptr< xAOD::Vertex > Trk::TrkV0VertexFitter::fit ( const EventContext & ctx,
const std::vector< const Trk::TrackParameters * > & perigeeList,
const std::vector< double > & masses,
const double & constraintMass,
const xAOD::Vertex * pointingVertex,
const Amg::Vector3D & startingPoint ) const
virtual

Interface for Trk::TrackParameters with mass and pointing constraints.

Definition at line 232 of file TrkV0VertexFitter.cxx.

238 {
239 if ( originalPerigees.empty() )
240 {
241 ATH_MSG_DEBUG("No tracks to fit in this event.");
242 return nullptr;
243 }
244
245 // Initialisation of variables
246 bool pointingConstraint = false;
247 bool massConstraint = false;
248 if(constraintMass > -100.) massConstraint = true;
249 bool conversion = false;
250 if(constraintMass == 0. && originalPerigees.size() == 2) conversion = true;
251 double x_point=0., y_point=0., z_point=0.;
252 AmgSymMatrix(3) pointingVertexCov; pointingVertexCov.setIdentity();
253 if (pointingVertex != nullptr) {
254 if (pointingVertex->covariancePosition().trace() != 0.) {
255 pointingConstraint = true;
256 Amg::Vector3D pv = pointingVertex->position();
257 x_point = pv.x();
258 y_point = pv.y();
259 z_point = pv.z();
260 pointingVertexCov = pointingVertex->covariancePosition().inverse();
261 }
262 }
263
264 if (msgLvl(MSG::DEBUG)) {
265 msg(MSG::DEBUG) << "massConstraint " << massConstraint << " pointingConstraint " << pointingConstraint << " conversion " << conversion << endmsg;
266 msg(MSG::DEBUG) << "V0Fitter called with: " << endmsg;
267 if (massConstraint && !masses.empty()) msg(MSG::DEBUG) << "mass constraint, V0Mass = " << constraintMass << " particle masses " << masses << endmsg;
268 if (pointingConstraint) msg(MSG::DEBUG) << "pointing constraint, x = " << x_point << " y = " << y_point << " z = " << z_point << endmsg;
269 }
270
271 bool restartFit = true;
272 double chi2 = 2000000000000.;
273 unsigned int nTrk = originalPerigees.size(); // Number of tracks to fit
274 unsigned int nMeas = 5*nTrk; // Number of measurements
275 unsigned int nVert = 1; // Number of vertices
276
277 unsigned int nCnst = 2*nTrk; // Number of constraint equations
278 unsigned int nPntC = 2; // Contribution from pointing constraint in 2D
279 unsigned int nMass = 1; // Contribution from mass constraint
280
281 if (massConstraint) {
282 nCnst = nCnst + nMass;
283 }
284 if (pointingConstraint) {
285 nCnst = nCnst + nPntC;
286 nMeas = nMeas + 3;
287 nVert = nVert + 1;
288 }
289
290 unsigned int nPar = 5*nTrk + 3*nVert; // Number of parameters
291 int ndf = nMeas - (nPar - nCnst); // Number of degrees of freedom
292 if (ndf < 0) {ndf = 1;}
293
294 unsigned int dim = nCnst; //
295 unsigned int n_dim = nMeas; //
296
297 ATH_MSG_DEBUG("ndf " << ndf << " n_dim " << n_dim << " dim " << dim);
298
299 std::vector<V0FitterTrack> v0FitterTracks;
300
301 Amg::VectorX Y_vec(n_dim); Y_vec.setZero();
302 Amg::VectorX Y0_vec(n_dim); Y0_vec.setZero();
303 Amg::Vector3D A_vec; A_vec.setZero();
304
305 Amg::MatrixX Wmeas_mat(n_dim,n_dim); Wmeas_mat.setZero();
306 Amg::MatrixX Wmeas0_mat(n_dim,n_dim); Wmeas0_mat.setZero();
307 Amg::MatrixX Bjac_mat(dim,n_dim); Bjac_mat.setZero();
308 Amg::MatrixX Ajac_mat(dim,3); Ajac_mat.setZero();
309 Amg::MatrixX C11_mat(n_dim,n_dim); C11_mat.setZero();
310 Amg::MatrixX C22_mat(3,3); C22_mat.setZero();
311 Amg::MatrixX C21_mat(3,n_dim); C21_mat.setZero();
312 Amg::MatrixX C31_mat(dim,n_dim); C31_mat.setZero();
313 Amg::MatrixX C32_mat(dim,3); C32_mat.setZero();
314 Amg::MatrixX Wb_mat(dim,dim); Wb_mat.setZero();
315 Amg::MatrixX Btemp_mat(dim,n_dim); Btemp_mat.setZero();
316 Amg::MatrixX Atemp_mat(dim,3); Atemp_mat.setZero();
317 Amg::VectorX DeltaY_vec(n_dim); DeltaY_vec.setZero();
318 Amg::Vector3D DeltaA_vec; DeltaA_vec.setZero();
319 Amg::VectorX DeltaY0_vec(n_dim); DeltaY0_vec.setZero();
320 Amg::VectorX F_vec(dim); F_vec.setZero();
321 Amg::VectorX C_vec(dim); C_vec.setZero();
322 Amg::VectorX C_cor_vec(dim); C_cor_vec.setZero();
323 Amg::MatrixX V_mat(nPar,nPar); V_mat.setZero();
324 Amg::MatrixX Chi_vec(1,n_dim); Chi_vec.setZero();
325 AmgSymMatrix(1) Chi_mat; Chi_mat.setZero();
326 Amg::MatrixX ChiItr_vec(1,n_dim); ChiItr_vec.setZero();
327 AmgSymMatrix(1) ChiItr_mat; ChiItr_mat.setZero();
328 Amg::VectorX F_fac_vec(dim); F_fac_vec.setZero();
329
330 const Amg::Vector3D * globalPosition = &(firstStartingPoint);
331 ATH_MSG_DEBUG("globalPosition of starting point: " << (*globalPosition)[0] << ", " << (*globalPosition)[1] << ", " << (*globalPosition)[2]);
332
333 if (globalPosition->perp() > m_maxR && globalPosition->z() > m_maxZ) return nullptr;
334
335 SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, ctx};
336 if (!readHandle.isValid()) {
337 std::string msg = "Failed to retrieve magmnetic field conditions data ";
338 msg += m_fieldCacheCondObjInputKey.key();
339 throw std::runtime_error(msg);
340 }
341 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
342
343 MagField::AtlasFieldCache fieldCache;
344 fieldCondObj->getInitializedCache (fieldCache);
345
346 // magnetic field
347 double BField[3];
348 fieldCache.getField(globalPosition->data(),BField);
349 double B_z = BField[2]*299.792; // should be in GeV/mm
350 if (B_z == 0. || std::isnan(B_z)) {
351 ATH_MSG_DEBUG("Could not find a magnetic field different from zero: very very strange");
352 B_z = 0.60407; // Value in GeV/mm (ATLAS units)
353 } else {
354 ATH_MSG_VERBOSE("Magnetic field projection of z axis in the perigee position is: " << B_z << " GeV/mm ");
355 }
356// double B_z = 1.998*0.3;
357
358
359 v0FitterTracks.clear();
360 Trk::PerigeeSurface perigeeSurface(*globalPosition);
361 // Extrapolate the perigees to the startpoint of the fit
362 for (const Trk::TrackParameters* chargeParameters : originalPerigees)
363 {
364 if (chargeParameters != nullptr)
365 {
366 // Correct material changes
367 const Amg::Vector3D gMomentum = chargeParameters->momentum();
368 const Amg::Vector3D gDirection = chargeParameters->position() - *globalPosition;
369 const double extrapolationDirection = gMomentum.dot( gDirection );
371 if(extrapolationDirection > 0) mode = Trk::addNoise;
372 std::unique_ptr<const Trk::Perigee> extrapolatedPerigee(nullptr);
373
374 std::unique_ptr<const Trk::TrackParameters> tmp =
375 std::abs(chargeParameters->position().z()) > m_maxZ ? nullptr :
376 m_extrapolator->extrapolate(ctx,
377 *chargeParameters,
378 perigeeSurface,
380 true,
381 Trk::pion,
382 mode);
383
384 //if of right type we want to pass ownership
385 if (tmp && tmp->associatedSurface().type() == Trk::SurfaceType::Perigee) {
386 extrapolatedPerigee.reset(static_cast<const Trk::Perigee*>(tmp.release()));
387 }
388
389 if (extrapolatedPerigee == nullptr) {
390 ATH_MSG_DEBUG("Perigee was not extrapolated! Taking original one!");
391 const Trk::Perigee* tmpPerigee = dynamic_cast<const Trk::Perigee*>(chargeParameters);
392 if (tmpPerigee!=nullptr) extrapolatedPerigee = std::make_unique<Trk::Perigee>(*tmpPerigee);
393 else return nullptr;
394 }
395
396 // store track parameters at starting point
397 V0FitterTrack locV0FitterTrack;
398 locV0FitterTrack.TrkPar[0] = extrapolatedPerigee->parameters()[Trk::d0];
399 locV0FitterTrack.TrkPar[1] = extrapolatedPerigee->parameters()[Trk::z0];
400 locV0FitterTrack.TrkPar[2] = extrapolatedPerigee->parameters()[Trk::phi];
401 locV0FitterTrack.TrkPar[3] = extrapolatedPerigee->parameters()[Trk::theta];
402 locV0FitterTrack.TrkPar[4] = extrapolatedPerigee->parameters()[Trk::qOverP];
403 locV0FitterTrack.Wi_mat = extrapolatedPerigee->covariance()->inverse().eval();
404 locV0FitterTrack.originalPerigee = chargeParameters;
405 v0FitterTracks.push_back(locV0FitterTrack);
406 } else {
407 ATH_MSG_DEBUG("Track parameters are not charged tracks ... fit aborted");
408 return nullptr;
409 }
410 }
411
412 // Iterate fits until the fit criteria are met, or the number of max iterations is reached
413 double chi2New=0.; double chi2Old=chi2;
414 double sumConstr=0.;
415 bool onConstr = false;
416 Amg::Vector3D frameOrigin = firstStartingPoint;
417 Amg::Vector3D frameOriginItr = firstStartingPoint;
418 for (int itr=0; itr < m_maxIterations; ++itr)
419 {
420 ATH_MSG_DEBUG("Iteration number: " << itr);
421 if (!restartFit) chi2Old = chi2New;
422 chi2New = 0.;
423
424 if (restartFit)
425 {
426 // ===> loop over tracks
427 std::vector<V0FitterTrack>::iterator PTIter;
428 int i=0;
429 for (PTIter = v0FitterTracks.begin(); PTIter != v0FitterTracks.end() ; ++PTIter)
430 {
431 V0FitterTrack locP((*PTIter));
432 Wmeas0_mat.block<5,5>(5*i,5*i) = locP.Wi_mat;
433 Wmeas_mat.block<5,5>(5*i,5*i) = locP.Wi_mat;
434 for (int j=0; j<5; ++j) {
435 Y0_vec(j+5*i) = locP.TrkPar[j];
436 }
437 ++i;
438 }
439 if(pointingConstraint) {
440 Y0_vec(5*nTrk + 0) = x_point;
441 Y0_vec(5*nTrk + 1) = y_point;
442 Y0_vec(5*nTrk + 2) = z_point;
443 Wmeas0_mat.block<3,3>(5*nTrk,5*nTrk) = pointingVertexCov;
444 Wmeas_mat.block<3,3>(5*nTrk,5*nTrk) = pointingVertexCov;
445 }
446 Wmeas_mat = Wmeas_mat.inverse();
447 }
448
449 Y_vec = Y0_vec + DeltaY_vec;
450 A_vec = DeltaA_vec;
451
452 // check theta and phi ranges
453 for (unsigned int i=0; i<nTrk; ++i)
454 {
455 if ( fabs ( Y_vec(2+5*i) ) > 100. || fabs ( Y_vec(3+5*i) ) > 100. ) { return nullptr; }
456 while ( fabs ( Y_vec(2+5*i) ) > M_PI ) Y_vec(2+5*i) += ( Y_vec(2+5*i) > 0 ) ? -2*M_PI : 2*M_PI;
457 while ( Y_vec(3+5*i) > 2*M_PI ) Y_vec(3+5*i) -= 2*M_PI;
458 while ( Y_vec(3+5*i) < -M_PI ) Y_vec(3+5*i) += M_PI;
459 if ( Y_vec(3+5*i) > M_PI )
460 {
461 Y_vec(3+5*i) = 2*M_PI - Y_vec(3+5*i);
462 if ( Y_vec(2+5*i) >= 0 ) Y_vec(2+5*i) += ( Y_vec(2+5*i) >0 ) ? -M_PI : M_PI;
463 }
464 if ( Y_vec(3+5*i) < 0.0 )
465 {
466 Y_vec(3+5*i) = - Y_vec(3+5*i);
467 if ( Y_vec(2+5*i) >= 0 ) Y_vec(2+5*i) += ( Y_vec(2+5*i) >0 ) ? -M_PI : M_PI;
468 }
469 }
470
471 double SigE=0., SigPx=0., SigPy=0., SigPz=0., Px=0., Py=0., Pz=0.;
472 Amg::VectorX rho(nTrk), Phi(nTrk), charge(nTrk);
473 rho.setZero(); Phi.setZero(); charge.setZero();
474 Amg::VectorX d0Cor(nTrk), d0Fac(nTrk), xcphiplusysphi(nTrk), xsphiminusycphi(nTrk);
475 d0Cor.setZero(); d0Fac.setZero(); xcphiplusysphi.setZero(); xsphiminusycphi.setZero();
476 AmgVector(2) conv_sign;
477 conv_sign[0] = -1; conv_sign[1] = 1;
478 for (unsigned int i=0; i<nTrk; ++i)
479 {
480 charge[i] = (Y_vec(4+5*i) < 0.) ? -1. : 1.;
481 rho[i] = sin(Y_vec(3+5*i))/(B_z*Y_vec(4+5*i));
482 xcphiplusysphi[i] = A_vec(0)*cos(Y_vec(2+5*i))+A_vec(1)*sin(Y_vec(2+5*i));
483 xsphiminusycphi[i] = A_vec(0)*sin(Y_vec(2+5*i))-A_vec(1)*cos(Y_vec(2+5*i));
484 if(fabs(-xcphiplusysphi[i]/rho[i]) > 1.) return nullptr;
485 d0Cor[i] = 0.5*asin(-xcphiplusysphi[i]/rho[i]);
486 double d0Facsq = 1. - xcphiplusysphi[i]*xcphiplusysphi[i]/(rho[i]*rho[i]);
487 d0Fac[i] = (d0Facsq>0.) ? sqrt(d0Facsq) : 0;
488 Phi[i] = Y_vec(2+5*i) + 2.*d0Cor[i];
489
490 if(massConstraint && !masses.empty() && masses[i] != 0.){
491 SigE += sqrt(1./(Y_vec(4+5*i)*Y_vec(4+5*i)) + masses[i]*masses[i]);
492 SigPx += sin(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
493 SigPy += sin(Y_vec(3+5*i))*sin(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
494 SigPz += cos(Y_vec(3+5*i))*charge[i]/Y_vec(4+5*i);
495 }
496 Px += sin(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
497 Py += sin(Y_vec(3+5*i))*sin(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
498 Pz += cos(Y_vec(3+5*i))*charge[i]/Y_vec(4+5*i);
499 }
500
501 double FMass=0., dFMassdxs=0., dFMassdys=0., dFMassdzs=0.;
502 double FPxy=0., dFPxydxs=0., dFPxydys=0., dFPxydzs=0., dFPxydxp=0., dFPxydyp=0., dFPxydzp=0.;
503 double FPxz=0., dFPxzdxs=0., dFPxzdys=0., dFPxzdzs=0., dFPxzdxp=0., dFPxzdyp=0., dFPxzdzp=0.;
504 Amg::VectorX Fxy(nTrk), Fxz(nTrk), dFMassdPhi(nTrk);
505 Fxy.setZero(); Fxz.setZero(); dFMassdPhi.setZero();
506 Amg::VectorX drhodtheta(nTrk), drhodqOverP(nTrk), csplusbc(nTrk), ccminusbs(nTrk);
507 drhodtheta.setZero(); drhodqOverP.setZero(); csplusbc.setZero(); ccminusbs.setZero();
508 Amg::VectorX dFxydd0(nTrk), dFxydz0(nTrk), dFxydphi(nTrk), dFxydtheta(nTrk), dFxydqOverP(nTrk);
509 dFxydd0.setZero(); dFxydz0.setZero(); dFxydphi.setZero(); dFxydtheta.setZero(); dFxydqOverP.setZero();
510 Amg::VectorX dFxydxs(nTrk), dFxydys(nTrk), dFxydzs(nTrk);
511 dFxydxs.setZero(); dFxydys.setZero(); dFxydzs.setZero();
512 Amg::VectorX dFxzdd0(nTrk), dFxzdz0(nTrk), dFxzdphi(nTrk), dFxzdtheta(nTrk), dFxzdqOverP(nTrk);
513 dFxzdd0.setZero(); dFxzdz0.setZero(); dFxzdphi.setZero(); dFxzdtheta.setZero(); dFxzdqOverP.setZero();
514 Amg::VectorX dFxzdxs(nTrk), dFxzdys(nTrk), dFxzdzs(nTrk);
515 dFxzdxs.setZero(); dFxzdys.setZero(); dFxzdzs.setZero();
516 Amg::VectorX dFMassdd0(nTrk), dFMassdz0(nTrk), dFMassdphi(nTrk), dFMassdtheta(nTrk), dFMassdqOverP(nTrk);
517 dFMassdd0.setZero(); dFMassdz0.setZero(); dFMassdphi.setZero(); dFMassdtheta.setZero(); dFMassdqOverP.setZero();
518 Amg::VectorX dFPxydd0(nTrk), dFPxydz0(nTrk), dFPxydphi(nTrk), dFPxydtheta(nTrk), dFPxydqOverP(nTrk);
519 dFPxydd0.setZero(); dFPxydz0.setZero(); dFPxydphi.setZero(); dFPxydtheta.setZero(); dFPxydqOverP.setZero();
520 Amg::VectorX dFPxzdd0(nTrk), dFPxzdz0(nTrk), dFPxzdphi(nTrk), dFPxzdtheta(nTrk), dFPxzdqOverP(nTrk);
521 dFPxzdd0.setZero(); dFPxzdz0.setZero(); dFPxzdphi.setZero(); dFPxzdtheta.setZero(); dFPxzdqOverP.setZero();
522 Amg::VectorX dPhidd0(nTrk), dPhidz0(nTrk), dPhidphi0(nTrk), dPhidtheta(nTrk), dPhidqOverP(nTrk);
523 dPhidd0.setZero(); dPhidz0.setZero(); dPhidphi0.setZero(); dPhidtheta.setZero(); dPhidqOverP.setZero();
524 Amg::VectorX dPhidxs(nTrk), dPhidys(nTrk), dPhidzs(nTrk);
525 dPhidxs.setZero(); dPhidys.setZero(); dPhidzs.setZero();
526 //
527 // constraint equations for V0vertex fitter
528 //
529 // FMass = mass vertex constraint
530 //
531 if (conversion) {
532 FMass = Phi[1] - Phi[0];
533 } else {
534 FMass = constraintMass*constraintMass - SigE*SigE + SigPx*SigPx + SigPy*SigPy + SigPz*SigPz;
535 }
536 //
537 // FPxy = pointing constraint in xy
538 //
539 FPxy = Px*(frameOriginItr[1] - y_point) - Py*(frameOriginItr[0]- x_point);
540 //
541 // FPxz = pointing constraint in xz
542 //
543 FPxz = Px*(frameOriginItr[2] - z_point) - Pz*(frameOriginItr[0]- x_point);
544
545 for (unsigned int i=0; i<nTrk; ++i)
546 {
547 //
548 // Fxy = vertex constraint in xy plane (one for each track)
549 //
550 Fxy[i] = Y_vec(0+5*i) + xsphiminusycphi[i] - 2.*rho[i]*sin(d0Cor[i])*sin(d0Cor[i]);
551 //
552 // Fxz = vertex constraint in xz plane (one for each track)
553 //
554 Fxz[i] = Y_vec(1+5*i) - A_vec(2) - rho[i]*2.*d0Cor[i]/tan(Y_vec(3+5*i));
555 //
556 // derivatives
557 //
558 drhodtheta[i] = cos(Y_vec(3+5*i))/(B_z*Y_vec(4+5*i));
559 drhodqOverP[i] = -sin(Y_vec(3+5*i))/(B_z*Y_vec(4+5*i)*Y_vec(4+5*i));
560
561 dFxydd0[i] = 1.;
562 dFxydphi[i] = xcphiplusysphi[i]*(1. + xsphiminusycphi[i]/(d0Fac[i]*rho[i]));
563 dFxydtheta[i] = (xcphiplusysphi[i]*xcphiplusysphi[i]/(d0Fac[i]*rho[i]*rho[i])-2.*sin(d0Cor[i])*sin(d0Cor[i]))*drhodtheta[i];
564 dFxydqOverP[i] = (xcphiplusysphi[i]*xcphiplusysphi[i]/(d0Fac[i]*rho[i]*rho[i])-2.*sin(d0Cor[i])*sin(d0Cor[i]))*drhodqOverP[i];
565 dFxydxs[i] = sin(Y_vec(2+5*i)) - cos(Y_vec(2+5*i))*xcphiplusysphi[i]/(d0Fac[i]*rho[i]);
566 dFxydys[i] = -cos(Y_vec(2+5*i)) - sin(Y_vec(2+5*i))*xcphiplusysphi[i]/(d0Fac[i]*rho[i]);
567
568 dFxzdz0[i] = 1.;
569 dFxzdphi[i] = -xsphiminusycphi[i]/(d0Fac[i]*tan(Y_vec(3+5*i)));
570 dFxzdtheta[i] = -((xcphiplusysphi[i]/(d0Fac[i]*rho[i]) + 2.*d0Cor[i])*tan(Y_vec(3+5*i))*drhodtheta[i] -
571 rho[i]*2.*d0Cor[i]/(cos(Y_vec(3+5*i))*cos(Y_vec(3+5*i))))/(tan(Y_vec(3+5*i))*tan(Y_vec(3+5*i)));
572 dFxzdqOverP[i] = -(xcphiplusysphi[i]/(d0Fac[i]*rho[i]) + 2.*d0Cor[i])*drhodqOverP[i]/tan(Y_vec(3+5*i));
573 dFxzdxs[i] = cos(Y_vec(2+5*i))/(d0Fac[i]*tan(Y_vec(3+5*i)));
574 dFxzdys[i] = sin(Y_vec(2+5*i))/(d0Fac[i]*tan(Y_vec(3+5*i)));
575 dFxzdzs[i] = -1.;
576
577 dPhidphi0[i] = 1. + xsphiminusycphi[i]/(d0Fac[i]*rho[i]);
578 dPhidtheta[i] = xcphiplusysphi[i]*drhodtheta[i]/(d0Fac[i]*rho[i]*rho[i]);
579 dPhidqOverP[i] = xcphiplusysphi[i]*drhodqOverP[i]/(d0Fac[i]*rho[i]*rho[i]);
580 dPhidxs[i] = -cos(Y_vec(2+5*i))/(d0Fac[i]*rho[i]);
581 dPhidys[i] = -sin(Y_vec(2+5*i))/(d0Fac[i]*rho[i]);
582
583 if (massConstraint && !masses.empty() && masses[i] != 0.){
584 if (conversion) {
585 dFMassdphi[i] = conv_sign[i]*dPhidphi0[i];
586 dFMassdtheta[i] = conv_sign[i]*dPhidtheta[i];
587 dFMassdqOverP[i] = conv_sign[i]*dPhidqOverP[i];
588 dFMassdxs += conv_sign[i]*dPhidxs[i];
589 dFMassdys += conv_sign[i]*dPhidys[i];
590 } else {
591 csplusbc[i] = SigPy*sin(Y_vec(2+5*i))+SigPx*cos(Y_vec(2+5*i));
592 ccminusbs[i] = SigPy*cos(Y_vec(2+5*i))-SigPx*sin(Y_vec(2+5*i));
593 dFMassdphi[i] = 2.*sin(Y_vec(3+5*i))*ccminusbs[i]*charge[i]/Y_vec(4+5*i);
594 dFMassdtheta[i] = 2.*(cos(Y_vec(3+5*i))*csplusbc[i] - sin(Y_vec(3+5*i))*SigPz)*charge[i]/Y_vec(4+5*i);
595 dFMassdqOverP[i] = 2.*SigE/(sqrt(1./(Y_vec(4+5*i)*Y_vec(4+5*i)) + masses[i]*masses[i])*Y_vec(4+5*i)*Y_vec(4+5*i)*Y_vec(4+5*i)) -
596 2.*charge[i]*(sin(Y_vec(3+5*i))*csplusbc[i] + cos(Y_vec(3+5*i))*SigPz)/(Y_vec(4+5*i)*Y_vec(4+5*i));
597 }
598 }
599
600 if (pointingConstraint){
601 dFPxydphi[i] = -sin(Y_vec(3+5*i))*(sin(Y_vec(2+5*i))*(frameOriginItr[1]-y_point)+cos(Y_vec(2+5*i))*(frameOriginItr[0]-x_point))*charge[i]/Y_vec(4+5*i);
602 dFPxydtheta[i] = cos(Y_vec(3+5*i))*(cos(Y_vec(2+5*i))*(frameOriginItr[1]-y_point)-sin(Y_vec(2+5*i))*(frameOriginItr[0]-x_point))*charge[i]/Y_vec(4+5*i);
603 dFPxydqOverP[i] = -sin(Y_vec(3+5*i))*(cos(Y_vec(2+5*i))*(frameOriginItr[1]-y_point)-sin(Y_vec(2+5*i))*(frameOriginItr[0]-x_point))*charge[i]/(Y_vec(4+5*i)*Y_vec(4+5*i));
604 dFPxydxs += -sin(Y_vec(3+5*i))*sin(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
605 dFPxydys += sin(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
606 dFPxydxp += sin(Y_vec(3+5*i))*sin(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
607 dFPxydyp += -sin(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
608
609 dFPxzdphi[i] = -sin(Y_vec(3+5*i))*sin(Y_vec(2+5*i))*(frameOriginItr[2]-z_point)*charge[i]/Y_vec(4+5*i);
610 dFPxzdtheta[i] = cos(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*(frameOriginItr[2]-z_point)*charge[i]/Y_vec(4+5*i)
611 +sin(Y_vec(3+5*i))*(frameOriginItr[0]-x_point)*charge[i]/Y_vec(4+5*i);
612 dFPxzdqOverP[i] = -sin(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*(frameOriginItr[2]-z_point)*charge[i]/(Y_vec(4+5*i)*Y_vec(4+5*i))
613 +cos(Y_vec(3+5*i))*(frameOriginItr[0]-x_point)*charge[i]/(Y_vec(4+5*i)*Y_vec(4+5*i));
614 dFPxzdxs += -cos(Y_vec(3+5*i))*charge[i]/Y_vec(4+5*i);
615 dFPxzdzs += sin(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
616 dFPxzdxp += cos(Y_vec(3+5*i))*charge[i]/Y_vec(4+5*i);
617 dFPxzdzp += -sin(Y_vec(3+5*i))*cos(Y_vec(2+5*i))*charge[i]/Y_vec(4+5*i);
618 }
619
620 // fill vector of constraints
621 F_vec[i] = -Fxy[i];
622 F_vec[i+nTrk] = -Fxz[i];
623 F_fac_vec[i] = 1.;
624 F_fac_vec[i+nTrk] = 1.;
625 }
626 if(massConstraint) F_vec(2*nTrk+0) = -FMass;
627 //if(massConstraint) F_fac_vec(2*nTrk+0) = 1.;
628 if(massConstraint) F_fac_vec(2*nTrk+0) = 0.000001;
629 if(pointingConstraint) {
630 if(massConstraint) {
631 F_vec(2*nTrk+1) = -FPxy;
632 F_vec(2*nTrk+2) = -FPxz;
633 F_fac_vec(2*nTrk+1) = 0.000001;
634 F_fac_vec(2*nTrk+2) = 0.000001;
635 } else {
636 F_vec(2*nTrk+0) = -FPxy;
637 F_vec(2*nTrk+1) = -FPxz;
638 F_fac_vec(2*nTrk+0) = 0.000001;
639 F_fac_vec(2*nTrk+1) = 0.000001;
640 }
641 }
642
643 sumConstr = 0.;
644 for (unsigned int i=0; i<dim; ++i)
645 {
646 sumConstr += F_fac_vec[i]*fabs(F_vec[i]);
647 }
648 if ( std::isnan(sumConstr) ) { return nullptr; }
649 if (sumConstr < 0.001) { onConstr = true; }
650 ATH_MSG_DEBUG("sumConstr " << sumConstr);
651
652 for (unsigned int i=0; i<nTrk; ++i)
653 {
654 Bjac_mat(i,0+5*i) = dFxydd0(i);
655 Bjac_mat(i,1+5*i) = dFxydz0(i);
656 Bjac_mat(i,2+5*i) = dFxydphi(i);
657 Bjac_mat(i,3+5*i) = dFxydtheta(i);
658 Bjac_mat(i,4+5*i) = dFxydqOverP(i);
659 Bjac_mat(i+nTrk,0+5*i) = dFxzdd0(i);
660 Bjac_mat(i+nTrk,1+5*i) = dFxzdz0(i);
661 Bjac_mat(i+nTrk,2+5*i) = dFxzdphi(i);
662 Bjac_mat(i+nTrk,3+5*i) = dFxzdtheta(i);
663 Bjac_mat(i+nTrk,4+5*i) = dFxzdqOverP(i);
664 if(massConstraint) {
665 Bjac_mat(2*nTrk,0+5*i) = dFMassdd0(i);
666 Bjac_mat(2*nTrk,1+5*i) = dFMassdz0(i);
667 Bjac_mat(2*nTrk,2+5*i) = dFMassdphi(i);
668 Bjac_mat(2*nTrk,3+5*i) = dFMassdtheta(i);
669 Bjac_mat(2*nTrk,4+5*i) = dFMassdqOverP(i);
670 }
671 if(pointingConstraint) {
672 if(massConstraint) {
673 Bjac_mat(2*nTrk+1,0+5*i) = dFPxydd0(i);
674 Bjac_mat(2*nTrk+1,1+5*i) = dFPxydz0(i);
675 Bjac_mat(2*nTrk+1,2+5*i) = dFPxydphi(i);
676 Bjac_mat(2*nTrk+1,3+5*i) = dFPxydtheta(i);
677 Bjac_mat(2*nTrk+1,4+5*i) = dFPxydqOverP(i);
678 Bjac_mat(2*nTrk+1,5*nTrk) = dFPxydxp;
679 Bjac_mat(2*nTrk+1,5*nTrk+1) = dFPxydyp;
680 Bjac_mat(2*nTrk+1,5*nTrk+2) = dFPxydzp;
681 Bjac_mat(2*nTrk+2,0+5*i) = dFPxzdd0(i);
682 Bjac_mat(2*nTrk+2,1+5*i) = dFPxzdz0(i);
683 Bjac_mat(2*nTrk+2,2+5*i) = dFPxzdphi(i);
684 Bjac_mat(2*nTrk+2,3+5*i) = dFPxzdtheta(i);
685 Bjac_mat(2*nTrk+2,4+5*i) = dFPxzdqOverP(i);
686 Bjac_mat(2*nTrk+2,5*nTrk) = dFPxzdxp;
687 Bjac_mat(2*nTrk+2,5*nTrk+1) = dFPxzdyp;
688 Bjac_mat(2*nTrk+2,5*nTrk+2) = dFPxzdzp;
689 } else {
690 Bjac_mat(2*nTrk+0,0+5*i) = dFPxydd0(i);
691 Bjac_mat(2*nTrk+0,1+5*i) = dFPxydz0(i);
692 Bjac_mat(2*nTrk+0,2+5*i) = dFPxydphi(i);
693 Bjac_mat(2*nTrk+0,3+5*i) = dFPxydtheta(i);
694 Bjac_mat(2*nTrk+0,4+5*i) = dFPxydqOverP(i);
695 Bjac_mat(2*nTrk+0,5*nTrk) = dFPxydxp;
696 Bjac_mat(2*nTrk+0,5*nTrk+1) = dFPxydyp;
697 Bjac_mat(2*nTrk+0,5*nTrk+2) = dFPxydzp;
698 Bjac_mat(2*nTrk+1,0+5*i) = dFPxzdd0(i);
699 Bjac_mat(2*nTrk+1,1+5*i) = dFPxzdz0(i);
700 Bjac_mat(2*nTrk+1,2+5*i) = dFPxzdphi(i);
701 Bjac_mat(2*nTrk+1,3+5*i) = dFPxzdtheta(i);
702 Bjac_mat(2*nTrk+1,4+5*i) = dFPxzdqOverP(i);
703 Bjac_mat(2*nTrk+1,5*nTrk) = dFPxzdxp;
704 Bjac_mat(2*nTrk+1,5*nTrk+1) = dFPxzdyp;
705 Bjac_mat(2*nTrk+1,5*nTrk+2) = dFPxzdzp;
706 }
707 }
708
709 Ajac_mat(i,0) = dFxydxs(i);
710 Ajac_mat(i,1) = dFxydys(i);
711 Ajac_mat(i,2) = dFxydzs(i);
712 Ajac_mat(i+nTrk,0) = dFxzdxs(i);
713 Ajac_mat(i+nTrk,1) = dFxzdys(i);
714 Ajac_mat(i+nTrk,2) = dFxzdzs(i);
715 if(massConstraint) {
716 Ajac_mat(2*nTrk,0) = dFMassdxs;
717 Ajac_mat(2*nTrk,1) = dFMassdys;
718 Ajac_mat(2*nTrk,2) = dFMassdzs;
719 }
720 if(pointingConstraint) {
721 if(massConstraint) {
722 Ajac_mat(2*nTrk+1,0) = dFPxydxs;
723 Ajac_mat(2*nTrk+1,1) = dFPxydys;
724 Ajac_mat(2*nTrk+1,2) = dFPxydzs;
725 Ajac_mat(2*nTrk+2,0) = dFPxzdxs;
726 Ajac_mat(2*nTrk+2,1) = dFPxzdys;
727 Ajac_mat(2*nTrk+2,2) = dFPxzdzs;
728 } else {
729 Ajac_mat(2*nTrk+0,0) = dFPxydxs;
730 Ajac_mat(2*nTrk+0,1) = dFPxydys;
731 Ajac_mat(2*nTrk+0,2) = dFPxydzs;
732 Ajac_mat(2*nTrk+1,0) = dFPxzdxs;
733 Ajac_mat(2*nTrk+1,1) = dFPxzdys;
734 Ajac_mat(2*nTrk+1,2) = dFPxzdzs;
735 }
736 }
737 }
738
739 Wb_mat = Wmeas_mat.similarity(Bjac_mat) ;
740 Wb_mat = Wb_mat.inverse();
741
742 C22_mat = Wb_mat.similarity(Ajac_mat.transpose());
743 C22_mat = C22_mat.inverse();
744
745 Btemp_mat = Wb_mat * Bjac_mat * Wmeas_mat;
746 Atemp_mat = Wb_mat * Ajac_mat;
747
748 C21_mat = - C22_mat * Ajac_mat.transpose() * Btemp_mat;
749 C32_mat = Atemp_mat * C22_mat;
750 C31_mat = Btemp_mat + Atemp_mat * C21_mat;
751 Amg::MatrixX mat_prod_1 = Wmeas_mat * Bjac_mat.transpose();
752 Amg::MatrixX mat_prod_2 = Wmeas_mat * Bjac_mat.transpose() * Wb_mat * Ajac_mat;
753 C11_mat = Wmeas_mat - Wb_mat.similarity( mat_prod_1 ) + C22_mat.similarity( mat_prod_2 );
754
755 C_cor_vec = Ajac_mat*DeltaA_vec + Bjac_mat*DeltaY_vec;
756 C_vec = C_cor_vec + F_vec;
757
758 DeltaY_vec = C31_mat.transpose()*C_vec;
759 DeltaA_vec = C32_mat.transpose()*C_vec;
760
761 for (unsigned int i=0; i<n_dim; ++i)
762 {
763 ChiItr_vec(0,i) = DeltaY_vec(i);
764 }
765 ChiItr_mat = Wmeas0_mat.similarity( ChiItr_vec );
766 chi2New = ChiItr_mat(0,0);
767
768 // current vertex position in global coordinates
769 frameOriginItr[0] += DeltaA_vec(0);
770 frameOriginItr[1] += DeltaA_vec(1);
771 frameOriginItr[2] += DeltaA_vec(2);
772 if (msgLvl(MSG::DEBUG)) {
773 msg(MSG::DEBUG) << "New vertex, global coordinates: " << frameOriginItr.transpose() << endmsg;
774 msg(MSG::DEBUG) << "chi2Old: " << chi2Old << " chi2New: " << chi2New << " fabs(chi2Old-chi2New): " << fabs(chi2Old-chi2New) << endmsg;
775 }
776
777 const Amg::Vector3D * globalPositionItr = &frameOriginItr;
778 if (globalPositionItr->perp() > m_maxR && globalPositionItr->z() > m_maxZ) return nullptr;
779
780 if (onConstr && fabs(chi2Old-chi2New) < 0.1) { break; }
781
782 double BFieldItr[3];
783 fieldCache.getField(globalPositionItr->data(),BFieldItr);
784 double B_z_new = BFieldItr[2]*299.792; // should be in GeV/mm
785 if (B_z_new == 0. || std::isnan(B_z_new)) {
786 ATH_MSG_DEBUG("Using old B_z");
787 B_z_new = B_z;
788 }
789
790 restartFit = false;
791 double deltaR = sqrt(DeltaA_vec(0)*DeltaA_vec(0)+DeltaA_vec(1)*DeltaA_vec(1)+DeltaA_vec(2)*DeltaA_vec(2));
792 double deltaB_z = fabs(B_z-B_z_new)/B_z;
793 bool changeBz = false;
794
795 if (m_deltaR) {
796 if (deltaR > 5. && itr < m_maxIterations-1) changeBz = true;
797 } else {
798 if (deltaB_z > 0.000001 && itr < m_maxIterations-1) changeBz = true;
799 }
800
801 if (changeBz) {
802 B_z = B_z_new;
803
804 v0FitterTracks.clear();
805 Trk::PerigeeSurface perigeeSurfaceItr(*globalPositionItr);
806 // Extrapolate the perigees to the new startpoint of the fit
807 for (const Trk::TrackParameters* chargeParameters : originalPerigees)
808 {
809 if (chargeParameters != nullptr)
810 {
811 // Correct material changes
812 const Amg::Vector3D gMomentum = chargeParameters->momentum();
813 const Amg::Vector3D gDirection = chargeParameters->position() - *globalPositionItr;
814 const double extrapolationDirection = gMomentum .dot( gDirection );
816 if(extrapolationDirection > 0) mode = Trk::addNoise;
817 std::unique_ptr<const Trk::Perigee> extrapolatedPerigee(nullptr);
818
819 std::unique_ptr<const Trk::TrackParameters> tmp =
820 std::abs(chargeParameters->position().z()) > m_maxZ ? nullptr :
821 m_extrapolator->extrapolate(ctx,
822 *chargeParameters,
823 perigeeSurfaceItr,
825 true,
826 Trk::pion,
827 mode);
828
829 // if of right type we want to pass ownership
830 if (tmp && tmp->associatedSurface().type() == Trk::SurfaceType::Perigee) {
831 extrapolatedPerigee.reset(
832 static_cast<const Trk::Perigee*>(tmp.release()));
833 }
834
835 if (extrapolatedPerigee == nullptr) {
836 ATH_MSG_DEBUG("Perigee was not extrapolated! Taking original one!");
837 const Trk::Perigee* tmpPerigee = dynamic_cast<const Trk::Perigee*>(chargeParameters);
838 if (tmpPerigee!=nullptr) extrapolatedPerigee = std::make_unique<Trk::Perigee>(*tmpPerigee);
839 else return nullptr;
840 }
841
842 // store track parameters at new starting point
843 V0FitterTrack locV0FitterTrack;
844 locV0FitterTrack.TrkPar[0] = extrapolatedPerigee->parameters()[Trk::d0];
845 locV0FitterTrack.TrkPar[1] = extrapolatedPerigee->parameters()[Trk::z0];
846 locV0FitterTrack.TrkPar[2] = extrapolatedPerigee->parameters()[Trk::phi];
847 locV0FitterTrack.TrkPar[3] = extrapolatedPerigee->parameters()[Trk::theta];
848 locV0FitterTrack.TrkPar[4] = extrapolatedPerigee->parameters()[Trk::qOverP];
849 locV0FitterTrack.Wi_mat = extrapolatedPerigee->covariance()->inverse().eval();
850 locV0FitterTrack.originalPerigee = chargeParameters;
851 v0FitterTracks.push_back(locV0FitterTrack);
852 } else {
853 ATH_MSG_DEBUG("Track parameters are not charged tracks ... fit aborted");
854 return nullptr;
855 }
856 }
857 frameOrigin = frameOriginItr;
858 Y0_vec *= 0.;
859 Y_vec *= 0.;
860 A_vec *= 0.;
861 DeltaY_vec *= 0.;
862 DeltaA_vec *= 0.;
863 chi2Old = 2000000000000.;
864 chi2New = 0.;
865 sumConstr = 0.;
866 onConstr = false;
867 restartFit = true;
868 }
869
870 //if (onConstr && fabs(chi2Old-chi2New) < 0.1) { break; }
871
872 } // end of iteration
873
874 frameOrigin[0] += DeltaA_vec(0);
875 frameOrigin[1] += DeltaA_vec(1);
876 frameOrigin[2] += DeltaA_vec(2);
877 if ( std::isnan(frameOrigin[0]) || std::isnan(frameOrigin[1]) || std::isnan(frameOrigin[2]) ) return nullptr;
878
879 Y_vec = Y0_vec + DeltaY_vec;
880
881 // check theta and phi ranges
882 for (unsigned int i=0; i<nTrk; ++i)
883 {
884 if ( fabs ( Y_vec(2+5*i) ) > 100. || fabs ( Y_vec(3+5*i) ) > 100. ) { return nullptr; }
885 while ( fabs ( Y_vec(2+5*i) ) > M_PI ) Y_vec(2+5*i) += ( Y_vec(2+5*i) > 0 ) ? -2*M_PI : 2*M_PI;
886 while ( Y_vec(3+5*i) > 2*M_PI ) Y_vec(3+5*i) -= 2*M_PI;
887 while ( Y_vec(3+5*i) < -M_PI ) Y_vec(3+5*i) += M_PI;
888 if ( Y_vec(3+5*i) > M_PI )
889 {
890 Y_vec(3+5*i) = 2*M_PI - Y_vec(3+5*i);
891 if ( Y_vec(2+5*i) >= 0 ) Y_vec(2+5*i) += ( Y_vec(2+5*i) >0 ) ? -M_PI : M_PI;
892 }
893 if ( Y_vec(3+5*i) < 0.0 )
894 {
895 Y_vec(3+5*i) = - Y_vec(3+5*i);
896 if ( Y_vec(2+5*i) >= 0 ) Y_vec(2+5*i) += ( Y_vec(2+5*i) >0 ) ? -M_PI : M_PI;
897 }
898 }
899
900 for (unsigned int i=0; i<n_dim; ++i)
901 {
902 Chi_vec(0,i) = DeltaY_vec(i);
903 }
904 Chi_mat = Wmeas0_mat.similarity( Chi_vec );
905 chi2 = Chi_mat(0,0);
906
907 V_mat.setZero();
908 V_mat.block(0,0,n_dim,n_dim) = C11_mat;
909 V_mat.block<3,3>(n_dim,n_dim) = C22_mat;
910 V_mat.block(n_dim,0,3,n_dim) = C21_mat;
911 V_mat.block(0,n_dim,n_dim,3) = C21_mat.transpose();
912
913 // ===> loop over tracks
914 std::vector<V0FitterTrack>::iterator BTIter;
915 int iRP=0;
916 for (BTIter = v0FitterTracks.begin(); BTIter != v0FitterTracks.end() ; ++BTIter)
917 {
918 // chi2 per track
919 AmgSymMatrix(5) covTrk = Wmeas0_mat.block<5,5>(5*iRP,5*iRP);
920 AmgVector(5) chi_vec; chi_vec.setZero();
921 for (unsigned int i=0; i<5; ++i) chi_vec(i) = DeltaY_vec(i+5*iRP);
922 double chi2Trk = chi_vec.dot(covTrk*chi_vec);
923 (*BTIter).chi2=chi2Trk;
924 iRP++;
925 }
926
927 // Store the vertex
928 auto vx = std::make_unique<xAOD::Vertex>();
929 vx->makePrivateStore();
930 vx->setPosition (frameOrigin);
931 vx->setCovariancePosition (C22_mat);
932 vx->setFitQuality(chi2,static_cast<float>(ndf));
933 vx->setVertexType(xAOD::VxType::V0Vtx);
934
935 // Store the tracks at vertex
936 std::vector<VxTrackAtVertex> & tracksAtVertex = vx->vxTrackAtVertex(); tracksAtVertex.clear();
937 Amg::Vector3D Vertex(frameOrigin[0],frameOrigin[1],frameOrigin[2]);
938 const Trk::PerigeeSurface Surface(Vertex);
939 Trk::Perigee * refittedPerigee(nullptr);
940 unsigned int iterf=0;
941 std::vector<V0FitterTrack>::iterator BTIterf;
942 for (BTIterf = v0FitterTracks.begin(); BTIterf != v0FitterTracks.end() ; ++BTIterf)
943 {
944 AmgSymMatrix(5) CovMtxP = V_mat.block<5,5>(5*iterf, 5*iterf);
945 refittedPerigee = new Trk::Perigee (Y_vec(0+5*iterf),Y_vec(1+5*iterf),Y_vec(2+5*iterf),Y_vec(3+5*iterf),Y_vec(4+5*iterf),
946 Surface, std::move(CovMtxP));
947 tracksAtVertex.emplace_back((*BTIterf).chi2, refittedPerigee, (*BTIterf).originalPerigee);
948 iterf++;
949 }
950
951 // Full Covariance Matrix
952 unsigned int sfcmv = nPar*(nPar+1)/2;
953 std::vector<float> floatErrMtx(sfcmv,0.);
954 unsigned int ipnt = 0;
955 for (unsigned int i=0; i<nPar; ++i) {
956 for (unsigned int j=0; j<i+1; ++j) {
957 floatErrMtx[ipnt++]=V_mat(i,j);
958 }
959 }
960 vx->setCovariance(floatErrMtx);
961
962 return vx;
963 }
#define M_PI
Scalar perp() const
perp method - perpendicular length
Scalar deltaR(const MatrixBase< Derived > &vec) const
#define ATH_MSG_VERBOSE(x)
double charge(const T &p)
Definition AtlasPID.h:997
#define AmgSymMatrix(dim)
#define AmgVector(rows)
boost::graph_traits< boost::adjacency_list< boost::vecS, boost::vecS, boost::bidirectionalS > >::vertex_descriptor Vertex
@ Phi
Definition RPCdef.h:8
if(pathvar)
void clear()
Empty the pool.
Eigen::Matrix< double, 3, 1 > Vector3D
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCacheCondObjInputKey
double chi2(TH1 *h0, TH1 *h1)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)
@ anyDirection
@ V0Vtx
Vertex from V0 Decay.
Definition VertexType.h:31
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
@ z
global position (cartesian)
Definition ParamDefs.h:57
@ theta
Definition ParamDefs.h:66
@ qOverP
perigee
Definition ParamDefs.h:67
@ phi
Definition ParamDefs.h:75
@ d0
Definition ParamDefs.h:63
@ z0
Definition ParamDefs.h:64
MaterialUpdateMode
This is a steering enum to force the material update it can be: (1) addNoise (-1) removeNoise Second ...
ParametersBase< TrackParametersDim, Charged > TrackParameters

◆ fit() [7/13]

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

Interface for TrackParameters with xAOD::Vertex starting point.

Interface for Trk::TrackParameters with xAOD::Vertex starting point.

Definition at line 211 of file TrkV0VertexFitter.cxx.

214 {
215 std::vector<double> masses;
216 double constraintMass = -9999.;
217 xAOD::Vertex * pointingVertex = nullptr;
218 const Amg::Vector3D& startingPoint = firstStartingPoint.position();
219 return fit(ctx, originalPerigees, masses, constraintMass, pointingVertex, startingPoint);
220 }
const Amg::Vector3D & position() const
Returns the 3-pos.

◆ fit() [8/13]

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

Fit interface for xAOD::TrackParticle with no starting point.

Interface for xAOD::TrackParticle with no starting point.

(0,0,0) will be assumed

Definition at line 105 of file TrkV0VertexFitter.cxx.

107 {
108 Amg::Vector3D tmpVtx;
109 tmpVtx.setZero();
110 return fit(ctx, vectorTrk, tmpVtx);
111 }

◆ fit() [9/13]

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

Interface for xAOD::TrackParticle with Amg::Vector3D starting point.

Definition at line 82 of file TrkV0VertexFitter.cxx.

85 {
86 std::vector<double> masses;
87 double constraintMass = -9999.;
88 xAOD::Vertex * pointingVertex = nullptr;
89 return fit(ctx, vectorTrk, masses, constraintMass, pointingVertex, firstStartingPoint);
90 }

◆ fit() [10/13]

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

Definition at line 58 of file TrkV0VertexFitter.h.

62 {
63 msg(MSG::WARNING)
64 << "TrkV0VertexFitter::fit(fit(const std::vector<const "
65 "TrackParticle*>&,const std::vector<const "
66 "Trk::NeutralParticle*>&,const Amg::Vector3D&) ignoring neutrals"
67 << endmsg;
68 return fit(ctx, vectorTrk, startingPoint);
69 };

◆ fit() [11/13]

virtual std::unique_ptr< xAOD::Vertex > Trk::TrkV0VertexFitter::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 76 of file TrkV0VertexFitter.h.

80 {
81 msg(MSG::WARNING)
82 << "TrkV0VertexFitter::fit(fit(const std::vector<const "
83 "TrackParticle*>&,const std::vector<const "
84 "Trk::NeutralParticle*>&,const xAOD::Vertex&) ignoring neutrals"
85 << endmsg;
86 return fit(ctx, vectorTrk, constraint);
87 };

◆ fit() [12/13]

std::unique_ptr< xAOD::Vertex > Trk::TrkV0VertexFitter::fit ( const EventContext & ctx,
const std::vector< const xAOD::TrackParticle * > & vectorTrk,
const std::vector< double > & masses,
const double & constraintMass,
const xAOD::Vertex * pointingVertex,
const Amg::Vector3D & startingPoint ) const
virtual

Methods specific for the V0 fitter Method taking a vector of tracks, vector of masses and starting point as arguments.

Interface for xAOD::TrackParticle with mass and pointing constraints.

Action: Fits a vertex out of initial set of tracks and applies a mass constraint to the particle reconstructed at the vertex and a pointing constraint: in the case the pointingVertex is a zero pointer, no pointing constraint applied.

Definition at line 114 of file TrkV0VertexFitter.cxx.

120 {
121 std::vector<const Trk::TrackParameters*> measuredPerigees;
122 std::vector<const Trk::TrackParameters*> measuredPerigees_delete;
123 for (const xAOD::TrackParticle* p : vectorTrk)
124 {
125 if (m_firstMeas) {
126 unsigned int indexFMP;
127 if (p->indexOfParameterAtPosition(indexFMP, xAOD::FirstMeasurement)) {
128 measuredPerigees.push_back(new CurvilinearParameters(p->curvilinearParameters(indexFMP)));
129 measuredPerigees_delete.push_back(measuredPerigees.back());
130 ATH_MSG_DEBUG("first measurement on track exists");
131 ATH_MSG_DEBUG("first measurement " << p->curvilinearParameters(indexFMP));
132 ATH_MSG_DEBUG("first measurement covariance " << *(p->curvilinearParameters(indexFMP)).covariance());
133 } else {
134 Amg::Transform3D CylTrf;
135 CylTrf.setIdentity();
136 Trk::CylinderSurface estimationCylinder(CylTrf, p->radiusOfFirstHit(), 10e10);
137 const Trk::TrackParameters* chargeParameters = &p->perigeeParameters();
139
140 const Trk::TrackParameters* extrapolatedPerigee =
141 std::abs(chargeParameters->position().z()) > m_maxZ ? nullptr :
142 m_extrapolator->extrapolate(ctx,
143 *chargeParameters,
144 estimationCylinder,
146 true,
147 Trk::pion,
148 mode).release();
149
150 if (extrapolatedPerigee != nullptr) {
151 ATH_MSG_DEBUG("extrapolated to first measurement");
152 measuredPerigees.push_back (extrapolatedPerigee);
153 measuredPerigees_delete.push_back (extrapolatedPerigee);
154 } else {
155
156 extrapolatedPerigee =
157 std::abs(chargeParameters->position().z()) > m_maxZ ? nullptr :
158 m_extrapolator->extrapolateDirectly(ctx,
159 *chargeParameters,
160 estimationCylinder,
162 true,
163 Trk::pion).release();
164
165 if (extrapolatedPerigee != nullptr) {
166 ATH_MSG_DEBUG( "extrapolated (direct) to first measurement");
167 measuredPerigees.push_back (extrapolatedPerigee);
168 measuredPerigees_delete.push_back (extrapolatedPerigee);
169 } else {
170 ATH_MSG_DEBUG("Failed to extrapolate to the first measurement on track, using Perigee parameters");
171 measuredPerigees.push_back (&p->perigeeParameters());
172 }
173 }
174 }
175 } else {
176 measuredPerigees.push_back (&p->perigeeParameters());
177 }
178 }
179
180 std::unique_ptr<xAOD::Vertex> fittedVxCandidate = fit(ctx, measuredPerigees, masses, constraintMass, pointingVertex, firstStartingPoint);
181
182 // assign the used tracks to the V0Candidate
183 if (fittedVxCandidate) {
184 for (const xAOD::TrackParticle* p : vectorTrk)
185 {
186 ElementLink<xAOD::TrackParticleContainer> el;
187 el.setElement(p);
188 fittedVxCandidate->addTrackAtVertex (el);
189 }
190 }
191
192 for (const auto *ptr : measuredPerigees_delete){ delete ptr; }
193
194 return fittedVxCandidate;
195 }
const Amg::Vector3D & position() const
Access method for the position.
Eigen::Affine3d Transform3D
@ alongMomentum
CurvilinearParametersT< TrackParametersDim, Charged, PlaneSurface > CurvilinearParameters
void * ptr(T *p)
Definition SGImplSvc.cxx:74
TrackParticle_v1 TrackParticle
Reference the current persistent version:
@ FirstMeasurement
Parameter defined at the position of the 1st measurement.

◆ fit() [13/13]

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

Interface for xAOD::TrackParticle with xAOD::Vertex starting point.

Definition at line 93 of file TrkV0VertexFitter.cxx.

96 {
97 std::vector<double> masses;
98 double constraintMass = -9999.;
99 xAOD::Vertex * pointingVertex = nullptr;
100 const Amg::Vector3D& startingPoint = firstStartingPoint.position();
101 return fit(ctx, vectorTrk, masses, constraintMass, pointingVertex, startingPoint);
102 }

◆ initialize()

StatusCode Trk::TrkV0VertexFitter::initialize ( )
overridevirtual

Definition at line 59 of file TrkV0VertexFitter.cxx.

60 {
61 if ( m_extrapolator.retrieve().isFailure() ) {
62 ATH_MSG_FATAL("Failed to retrieve tool " << m_extrapolator);
63 return StatusCode::FAILURE;
64 }
65 ATH_MSG_DEBUG( "Retrieved tool " << m_extrapolator );
66
67
69
70 ATH_MSG_DEBUG( "Initialize successful");
71 return StatusCode::SUCCESS;
72 }
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_FATAL(x)

Member Data Documentation

◆ m_deltaR

bool Trk::TrkV0VertexFitter::m_deltaR
private

Definition at line 188 of file TrkV0VertexFitter.h.

◆ m_extrapolator

ToolHandle< Trk::IExtrapolator > Trk::TrkV0VertexFitter::m_extrapolator
private

Data members to store the results.

Definition at line 192 of file TrkV0VertexFitter.h.

◆ m_fieldCacheCondObjInputKey

SG::ReadCondHandleKey<AtlasFieldCacheCondObj> Trk::TrkV0VertexFitter::m_fieldCacheCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"}
private

Definition at line 193 of file TrkV0VertexFitter.h.

194{this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};

◆ m_firstMeas

bool Trk::TrkV0VertexFitter::m_firstMeas
private

Definition at line 187 of file TrkV0VertexFitter.h.

◆ m_maxDchi2PerNdf

double Trk::TrkV0VertexFitter::m_maxDchi2PerNdf
private

Definition at line 184 of file TrkV0VertexFitter.h.

◆ m_maxIterations

int Trk::TrkV0VertexFitter::m_maxIterations
private

Definition at line 183 of file TrkV0VertexFitter.h.

◆ m_maxR

double Trk::TrkV0VertexFitter::m_maxR
private

Definition at line 185 of file TrkV0VertexFitter.h.

◆ m_maxZ

double Trk::TrkV0VertexFitter::m_maxZ
private

Definition at line 186 of file TrkV0VertexFitter.h.


The documentation for this class was generated from the following files: