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