35 #include "CLHEP/Vector/ThreeVector.h"
43 CLHEP::HepMatrix getPhiThetaQOverPToPxPyPzJacobian(
double qOverP,
double theta,
double phi) {
46 const double sinPhi(
sin(
phi));
48 const double cosPhi(
cos(
phi));
49 const double absQoP(fabs(
qOverP));
51 transform[0][0]=-1./absQoP*sinTheta*sinPhi;
52 transform[0][1]=1./absQoP*sinTheta*cosPhi;
54 transform[1][0]=1./absQoP*cosTheta*cosPhi;
55 transform[1][1]=1./absQoP*cosTheta*sinPhi;
71 m_linearizedTrackFactoryIsAvailable(false)
74 declareInterface< InDetJetFitterUtils >(
this) ;
89 return StatusCode::FAILURE;
101 if (
sc.isFailure()) {
103 return StatusCode::FAILURE;
107 return StatusCode::SUCCESS;
116 const AmgMatrix(5,3) &
A = linTrack->positionJacobian();
117 const AmgMatrix(5,3) &
B = linTrack->momentumJacobian();
118 const AmgSymMatrix(5) & trkParametersWeight = linTrack->expectedWeightAtPCA();
119 AmgSymMatrix(3) Sm = B.transpose() * trkParametersWeight * B;
122 if(Sm.determinant()==0){
123 msg(MSG::WARNING) <<
"Inversion of S matrix fails in track refit" <<
endmsg;
124 msg(MSG::WARNING) <<
" This track is returned not refitted" <<
endmsg;
125 throw std::string(
"Inversion of S matrix fails in track parameters refit");
127 AmgMatrix(3,3) posMomentumCovariance = -vrt_cov *
A.transpose() * trkParametersWeight * B *Sm;
141 msg(MSG::ERROR) <<
" No LinearizedTrackFactory defined. Cannot calculate compatibility. 0 compatibility returned" <<
endmsg;
142 return std::pair<double,double>(0,0);
149 vertexPosition[0]=
vertex.position()[0];
150 vertexPosition[1]=
vertex.position()[1];
151 vertexPosition[2]=
vertex.position()[2];
157 const AmgSymMatrix(5) & ExpectedCovariance=myLinearizedTrack->expectedCovarianceAtPCA();
159 AmgSymMatrix(2) weightReduced=ExpectedCovariance.block<2,2>(0,0);
163 weightReduced+=errorVertexReduced;
165 weightReduced.
inverse().eval();
168 if(weightReduced.determinant()==0)
170 msg(MSG::WARNING) <<
" Problem inverting cov matrix in compatibility method" <<
endmsg;
174 Amg::Vector2D paramsReduced((myLinearizedTrack->expectedParametersAtPCA())[0],(myLinearizedTrack->expectedParametersAtPCA())[1]);
176 double returnv2= paramsReduced.transpose() * weightReduced.inverse() * paramsReduced;
178 delete myLinearizedTrack;
179 myLinearizedTrack=
nullptr;
181 return std::pair<double,double>(returnv2,0);
190 msg(MSG::ERROR) <<
" No LinearizedTrackFactory defined. Cannot calculate compatibility. 0 compatibility returned" <<
endmsg;
191 return std::pair<double,double>(0,0);
198 vertexPosition[1]=
vertex.position()[1];
199 vertexPosition[2]=
vertex.position()[2];
204 const AmgSymMatrix(5) & ExpectedCovariance=myLinearizedTrack->expectedCovarianceAtPCA();
206 AmgSymMatrix(2) weightReduced=ExpectedCovariance.block<2,2>(0,0);
210 weightReduced+=errorVertexReduced;
212 weightReduced.
inverse().eval();
214 if(weightReduced.determinant()==0)
216 msg(MSG::WARNING) <<
" Problem inverting cov matrix in compatibility method" <<
endmsg;
218 Amg::Vector2D paramsReduced((myLinearizedTrack->expectedParametersAtPCA())[0],(myLinearizedTrack->expectedParametersAtPCA())[1]);
220 double returnv2= paramsReduced.transpose() * weightReduced * paramsReduced;
223 delete myLinearizedTrack;
224 myLinearizedTrack=
nullptr;
226 return std::pair<double,double>(returnv2,0);
242 double sign=(jetMomentum.cross(trackMom)).
dot(trackMom.cross(primaryPos-trackPos));
250 std::vector<Trk::VxTrackAtVertex*>::const_iterator vtxIter=myVxCandidate.
vxTrackAtVertex()->begin();
265 if (firstTrackPerigee==
nullptr ||secondTrackPerigee==
nullptr)
267 ATH_MSG_DEBUG(
"No Perigee in one of the two tracks at vertex. No sensible charge returned." );
274 fabs(firstTrackPerigee->parameters()[
Trk::qOverP])+
276 fabs(secondTrackPerigee->parameters()[
Trk::qOverP])+0.5 )
285 std::vector<Trk::VxTrackAtVertex>::const_iterator vtxIter=myVxCandidate.
vxTrackAtVertex().begin();
300 if (firstTrackPerigee==
nullptr ||secondTrackPerigee==
nullptr)
302 ATH_MSG_DEBUG(
"No Perigee in one of the two tracks at vertex. No sensible charge returned." );
309 fabs(firstTrackPerigee->parameters()[
Trk::qOverP])+
311 fabs(secondTrackPerigee->parameters()[
Trk::qOverP])+0.5 )
317 double highestMomMass,
318 double lowestMomMass)
const
323 std::vector<Trk::VxTrackAtVertex*>::const_iterator vtxIter=myVxCandidate.
vxTrackAtVertex()->begin();
338 if (firstTrackPerigee==
nullptr ||secondTrackPerigee==
nullptr)
340 ATH_MSG_DEBUG(
"No Perigee in one of the two tracks at vertex. No sensible mass returned." );
348 CLHEP::HepLorentzVector first4Mom;
349 CLHEP::HepLorentzVector second4Mom;
351 if (firstMomentum.mag2()>secondMomentum.mag2())
353 first4Mom=CLHEP::HepLorentzVector(firstMomentum.x(),firstMomentum.y(),firstMomentum.z(),TMath::Sqrt(highestMomMass*highestMomMass+firstMomentum.mag()*firstMomentum.mag()));
354 second4Mom=CLHEP::HepLorentzVector(secondMomentum.x(),secondMomentum.y(),secondMomentum.z(),TMath::Sqrt(lowestMomMass*lowestMomMass+secondMomentum.mag()*secondMomentum.mag()));
358 first4Mom=CLHEP::HepLorentzVector(firstMomentum.x(),firstMomentum.y(),firstMomentum.z(),TMath::Sqrt(lowestMomMass*highestMomMass+firstMomentum.mag()*firstMomentum.mag()));
359 second4Mom=CLHEP::HepLorentzVector(secondMomentum.x(),secondMomentum.y(),secondMomentum.z(),TMath::Sqrt(highestMomMass*lowestMomMass+secondMomentum.mag()*secondMomentum.mag()));
362 return (first4Mom+second4Mom).mag();
369 double highestMomMass,
370 double lowestMomMass)
const
375 std::vector<Trk::VxTrackAtVertex>::const_iterator vtxIter=myVxCandidate.
vxTrackAtVertex().begin();
390 if (firstTrackPerigee==
nullptr ||secondTrackPerigee==
nullptr)
392 ATH_MSG_DEBUG(
"No Perigee in one of the two tracks at vertex. No sensible mass returned." );
400 CLHEP::HepLorentzVector first4Mom;
401 CLHEP::HepLorentzVector second4Mom;
403 if (firstMomentum.mag2()>secondMomentum.mag2())
405 first4Mom=CLHEP::HepLorentzVector(firstMomentum.x(),firstMomentum.y(),firstMomentum.z(),TMath::Sqrt(highestMomMass*highestMomMass+firstMomentum.mag()*firstMomentum.mag()));
406 second4Mom=CLHEP::HepLorentzVector(secondMomentum.x(),secondMomentum.y(),secondMomentum.z(),TMath::Sqrt(lowestMomMass*lowestMomMass+secondMomentum.mag()*secondMomentum.mag()));
410 first4Mom=CLHEP::HepLorentzVector(firstMomentum.x(),firstMomentum.y(),firstMomentum.z(),TMath::Sqrt(lowestMomMass*highestMomMass+firstMomentum.mag()*firstMomentum.mag()));
411 second4Mom=CLHEP::HepLorentzVector(secondMomentum.x(),secondMomentum.y(),secondMomentum.z(),TMath::Sqrt(highestMomMass*lowestMomMass+secondMomentum.mag()*secondMomentum.mag()));
414 return (first4Mom+second4Mom).mag();
428 second.covariancePosition();
431 sumErrorsThenInverted.
inverse().eval();
438 if (sumErrorsThenInverted.determinant()>0 &&
distance>0)
440 double temp=differenceUnit.transpose() * sumErrorsThenInverted * differenceUnit;
443 error=1./std::sqrt(temp );
447 ATH_MSG_DEBUG(
"The significance of the distance to the PV is negative or zero definite: " << temp );
453 if (sumErrorsThenInverted.determinant()<=0)
455 ATH_MSG_DEBUG(
"Sum of cov matrices of PV + single vertex fit is zero or negative. Error on distance is returned as 1000mm." );
468 const EventContext& ctx = Gaudi::Hive::currentContext();
471 msg(MSG::ERROR) <<
"Cannot perform requested extrapolation. No extrapolator defined...Returning 0 compatibility..." <<
endmsg;
472 return std::pair<double,double>(0,0);
479 mySurface).release();
480 if (newMeasPerigee==
nullptr)
482 msg(MSG::WARNING) <<
" Extrapolation failed. Wrong d0 and z0 returned " <<
endmsg;
483 return std::pair<double,double>
484 (trackPerigee.parameters()[
Trk::d0],
488 double IPd0=newMeasPerigee->parameters()[
Trk::d0];
489 double IPz0=newMeasPerigee->parameters()[
Trk::z0]*
492 delete newMeasPerigee;
493 newMeasPerigee=
nullptr;
495 return std::pair<double,double>(IPd0,IPz0);
505 ATH_MSG_ERROR(
"Cannot perform requested extrapolation. No extrapolator defined...Returning 0 compatibility..." );
506 return std::pair<double,double>(0,0);
511 vertexPosition[0]=
vertex.position()[0];
513 vertexPosition[2]=
vertex.position()[2];
515 const AmgSymMatrix(5) & ExpectedCovariance=myLinearizedTrack->expectedCovarianceAtPCA();
516 AmgSymMatrix(2) weightReduced=ExpectedCovariance.block<2,2>(0,0);
520 double IPd0Sig=(myLinearizedTrack->expectedParametersAtPCA())[0]/sqrt( weightReduced(0,0) );
521 double IPz0Sig=(myLinearizedTrack->expectedParametersAtPCA())[1]/sqrt( weightReduced(1,1) );
532 myLinearizedTrack=
nullptr;
534 return std::pair<
double,
double>(IPd0Sig,IPz0Sig);
561 const std::vector<const Trk::ITrackLink*> & vectorOfTracks)
564 std::vector<const Trk::ITrackLink*>::const_iterator vectorOfTracksBegin=vectorOfTracks.begin();
565 std::vector<const Trk::ITrackLink*>::const_iterator vectorOfTracksEnd=vectorOfTracks.end();
567 for (std::vector<const Trk::ITrackLink*>::const_iterator vectorOfTracksIter=vectorOfTracksBegin;
568 vectorOfTracksIter!=vectorOfTracksEnd;++vectorOfTracksIter)
570 if (*vectorOfTracksIter==trackToCheck)
579 const std::vector<const Trk::LinkToTrackParticleBase*> & vectorOfTracks)
582 std::vector<const Trk::LinkToTrackParticleBase*>::const_iterator vectorOfTracksBegin=vectorOfTracks.begin();
583 std::vector<const Trk::LinkToTrackParticleBase*>::const_iterator vectorOfTracksEnd=vectorOfTracks.end();
585 for (std::vector<const Trk::LinkToTrackParticleBase*>::const_iterator vectorOfTracksIter=vectorOfTracksBegin;
586 vectorOfTracksIter!=vectorOfTracksEnd;++vectorOfTracksIter)
588 if (*vectorOfTracksIter==trackToCheck)
597 const std::vector<const Trk::LinkToTrackParticleBase*> & )
649 const std::vector<const xAOD::Vertex*> & vectorOfCandidates)
652 std::vector<const xAOD::Vertex*>::const_iterator vectorOfCandidatesBegin=vectorOfCandidates.begin();
653 std::vector<const xAOD::Vertex*>::const_iterator vectorOfCandidatesEnd=vectorOfCandidates.end();
655 for (std::vector<const xAOD::Vertex*>::const_iterator vectorOfCandidatesIter=vectorOfCandidatesBegin;
656 vectorOfCandidatesIter!=vectorOfCandidatesEnd;
657 ++vectorOfCandidatesIter)
659 if (*vectorOfCandidatesIter==vertexToCheck)
668 const std::vector<const xAOD::Vertex*> & vectorOfVxCandidates)
671 std::vector<const xAOD::Vertex*>::const_iterator verticesToVetoBegin=vectorOfVxCandidates.begin();
672 std::vector<const xAOD::Vertex*>::const_iterator verticesToVetoEnd=vectorOfVxCandidates.end();
674 for (std::vector<const xAOD::Vertex*>::const_iterator verticesToVetoIter=verticesToVetoBegin;
675 verticesToVetoIter!=verticesToVetoEnd;++verticesToVetoIter)
684 if (trackToCheck==linkToTP1||
685 trackToCheck==linkToTP2)
699 const double s_pion=139.57018;
702 CLHEP::HepLorentzVector massVector(0,0,0,0);
704 const std::vector<Trk::VxTrackAtVertex*> & tracksOfVertex=myVxVertexOnJetAxis.
getTracksAtVertex();
705 std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackBegin=tracksOfVertex.begin();
706 std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackEnd=tracksOfVertex.end();
707 for (std::vector<Trk::VxTrackAtVertex*>::const_iterator clustersOfTrackIter=clustersOfTrackBegin;
708 clustersOfTrackIter!=clustersOfTrackEnd;
709 ++clustersOfTrackIter)
711 if (
dynamic_cast<const Trk::Perigee*
>((*clustersOfTrackIter)->perigeeAtVertex())!=
nullptr)
716 massVector+=CLHEP::HepLorentzVector(mytrack.x(),mytrack.y(),mytrack.z(),TMath::Sqrt(s_pion*s_pion+mytrack.mag()*mytrack.mag()));