23 declareInterface<IMuonSegmentMomentumEstimator>(
this);
35 return StatusCode::SUCCESS;
44 ATH_MSG_VERBOSE(
" Executing MuonSegmentMomentumTool fitMomentumVectorSegments ");
45 ATH_MSG_DEBUG(
" fitMomentumVectorSegments " << segments.size() <<
" segments " );
47 std::vector<const Muon::MuonSegment*>::const_iterator it = segments.begin();
48 std::vector<const Muon::MuonSegment*>::const_iterator it2 = segments.begin();
49 std::vector<const Muon::MuonSegment*>::const_iterator it_end = segments.end();
63 ATH_MSG_DEBUG(
" Estimated signed momentum " << signedMomentum );
86 throw std::runtime_error(Form(
"File: %s, Line: %d\nMuonSegmentMomentumFromField::fieldIntegralEstimate() - Failed to retrieve AtlasFieldCacheCondObj with key %s", __FILE__, __LINE__, (
m_fieldCondObjInputKey.key()).c_str()));
90 fieldCache.
getField(point1.data(),field1.data());
91 fieldCache.
getField(point2.data(),field2.data());
92 fieldCache.
getField(point3.data(),field3.data());
94 field1[2]=field2[2]=field3[2]=0;
96 Amg::Vector3D crossvec=posdiff.unit().cross(field1+field2+field3);
98 double averagelcrossB=crossvec.mag()/3;
99 if (rphidir.x()*fieldsum.x()+rphidir.y()*fieldsum.y()<0) averagelcrossB=-averagelcrossB;
100 ATH_MSG_DEBUG(
"field integral " << averagelcrossB <<
"dist " << posdiff.mag() <<
" tot " << averagelcrossB*posdiff.mag());
101 return averagelcrossB*posdiff.mag();
110 ATH_MSG_VERBOSE(
" Executing MuonSegmentMomentumTool fitMomentum2Segments ");
119 double theta2=myseg2->globalDirection().theta();
121 double phi2=myseg2->globalDirection().phi();
122 double deltaphi=std::abs(phi2-phi1);
123 double deltatheta=99999999;
124 if (std::abs(deltaphi-2*
M_PI)<deltaphi) deltaphi=std::abs(deltaphi-2*
M_PI);
125 if (deltaphi<
M_PI/2) deltatheta=theta2-theta1;
127 if (theta1>
M_PI/2) deltatheta=2*
M_PI-theta2-theta1;
128 else deltatheta=-theta1-theta2;
130 double dist1=-1,dist2=-1;
131 const Trk::RIO_OnTrack *firstphi1=
nullptr,*lastphi1=
nullptr,*firstphi2=
nullptr,*lastphi2=
nullptr;
145 if (!firstphi1) firstphi1=rot;
149 for (
int i=0;i<(int)myseg2->numberOfMeasurementBases();i++){
159 if (!firstphi2) firstphi2=rot;
165 if (firstphi2) dist2=std::abs((firstphi2->globalPosition()-lastphi2->globalPosition()).dot(myseg2->globalDirection()));
172 for (
int itry=0;itry<2;itry++){
185 signedMomentum =-.3e3*fieldintegral/deltatheta;
186 if(std::abs(signedMomentum)<1000.) signedMomentum = 1e6;
187 ATH_MSG_DEBUG(
"integral: " << fieldintegral <<
" deltatheta: " << deltatheta <<
" signedmomentum : " << signedMomentum);
188 double residual = 9999.;
189 double resi[4],qoverp[4];
190 for (
int i=0;i<4;i++){
191 std::optional<Trk::TransportJacobian> jac{};
199 if (par && jac && (*jac)(1,4)!=0){
202 qoverp[i] = 1/signedMomentum;
204 double delta_qoverp=residual/(*jac)(1,4);
205 double der_simple = -10.*(bestseg->
globalPosition()-worstseg->globalPosition()).mag()/(.3*fieldintegral);
208 if(qoverp[i]!=qoverp[i-1]) {
209 double derivative = -(resi[i]-resi[i-1])/(qoverp[i]-qoverp[i-1]);
210 ATH_MSG_DEBUG(
" numerical derivative " << derivative <<
" derivative from track " << (*jac)(1,4) <<
" der_simple " << der_simple);
211 if(std::abs(derivative)>std::abs((*jac)(1,4))) {
212 ATH_MSG_DEBUG(
" use numerical derivative " << derivative <<
" derivative from track " << (*jac)(1,4));
213 delta_qoverp = residual/derivative;
217 if(std::abs(der_simple)>std::abs((*jac)(1,4))) {
218 ATH_MSG_DEBUG(
" use simple numerical derivative " << der_simple <<
" derivative from track " << (*jac)(1,4));
219 delta_qoverp = residual/der_simple;
222 ATH_MSG_DEBUG(
"residual: " << residual <<
" jac " << (*jac)(1,4) <<
" signedmomentum: " << signedMomentum <<
" delta_qoverp " << delta_qoverp);
223 double signedMomentum_updated = signedMomentum/(1+signedMomentum*delta_qoverp);
224 if(std::abs(signedMomentum_updated)<1000.) {
225 ATH_MSG_DEBUG(
"Too low signed momentum " << signedMomentum_updated );
227 signedMomentum = signedMomentum_updated>0? 1000.:-1000;
230 if(i<3) signedMomentum = signedMomentum_updated;
233 if (std::abs(residual)<1)
break;
236 if(std::abs(residual)>10.) {
237 ATH_MSG_DEBUG(
"NOT converged residual: " << residual <<
" itry " << itry);
238 if(itry==1)
ATH_MSG_DEBUG(
"NOT converged residual after two trials ");
264 throw std::runtime_error(Form(
"File: %s, Line: %d\nMuonSegmentMomentumFromField::fieldIntegralEstimate_old() - Failed to retrieve AtlasFieldCacheCondObj with key %s", __FILE__, __LINE__, (
m_fieldCondObjInputKey.key()).c_str()));
268 fieldCache.
getField(point1.data(),field1.data());
269 fieldCache.
getField(point2.data(),field2.data());
270 fieldCache.
getField(point3.data(),field3.data());
272 double averageBcrossl=(field1.cross(posdiff.unit()).mag()+field2.cross(posdiff.unit()).mag()+field3.cross(posdiff.unit()).mag())/3;
273 ATH_MSG_DEBUG(
"field integral " << averageBcrossl <<
"dist " << posdiff.mag() <<
" tot " << averageBcrossl*posdiff.mag());
274 return averageBcrossl*posdiff.mag();
283 ATH_MSG_VERBOSE(
" Executing MuonSegmentMomentumTool fitMomentum2Segments ");
292 double theta2=myseg2->globalDirection().theta();
294 double phi2=myseg2->globalDirection().phi();
295 double deltaphi=std::abs(phi2-phi1);
296 double deltatheta=99999999;
297 if (std::abs(deltaphi-2*
M_PI)<deltaphi) deltaphi=std::abs(deltaphi-2*
M_PI);
298 if (deltaphi<
M_PI/2) deltatheta=theta2-theta1;
300 if (theta1>
M_PI/2) deltatheta=2*
M_PI-theta2-theta1;
301 else deltatheta=-theta1-theta2;
303 double dist1=-1,dist2=-1;
304 const Trk::RIO_OnTrack *firstphi1=
nullptr,*lastphi1=
nullptr,*firstphi2=
nullptr,*lastphi2=
nullptr;
318 if (!firstphi1) firstphi1=rot;
322 for (
int i=0;i<(int)myseg2->numberOfMeasurementBases();i++){
332 if (!firstphi2) firstphi2=rot;
338 if (firstphi2) dist2=std::abs((firstphi2->globalPosition()-lastphi2->globalPosition()).dot(myseg2->globalDirection()));
343 signedMomentum =-.3e3*fieldintegral/deltatheta;
344 ATH_MSG_DEBUG(
"integral: " << fieldintegral <<
" deltatheta: " << deltatheta <<
" signedmomentum : " << signedMomentum);
345 for (
int i=0;i<3;i++){
348 std::optional<Trk::TransportJacobian> jac{};
353 if (par && jac && (*jac)(1,4)!=0){
354 double residual=worstseg->localParameters()[
Trk::locY] - par->parameters()[
Trk::locY];
355 double delta_qoverp=residual/(*jac)(1,4);
356 signedMomentum=1/(1/signedMomentum+delta_qoverp);
357 ATH_MSG_DEBUG(
"residual: " << residual <<
" jac " << (*jac)(1,4) <<
" dp " << delta_qoverp <<
" signedmomentum: " << signedMomentum);
359 if (std::abs(residual)<1)
break;
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
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,...
double fieldIntegralEstimate(const EventContext &ctx, const Muon::MuonSegment *segment1, const Muon::MuonSegment *segment2) const
double fieldIntegralEstimate_old(const EventContext &ctx, const Muon::MuonSegment *segment1, const Muon::MuonSegment *segment2) const
virtual void fitMomentum2Segments_old(const EventContext &ctx, const Muon::MuonSegment *segment1, const Muon::MuonSegment *segment2, double &signedMomentum) const
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCondObjInputKey
virtual void fitMomentum2Segments(const EventContext &ctx, const Muon::MuonSegment *segment1, const Muon::MuonSegment *segment2, double &signedMomentum) const
fits a momentum to 2 segments
ToolHandle< Trk::IPropagator > m_propagator
virtual void fitMomentumVectorSegments(const EventContext &, const std::vector< const Muon::MuonSegment * > &, double &signedMomentum) const
fits a momentum to a vector of segments
virtual StatusCode initialize()
to initiate private members
MuonSegmentMomentumFromField(const std::string &, const std::string &, const IInterface *)
constructor
ToolHandle< Trk::INavigator > m_navigator
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Gaudi::Property< bool > m_doOld
This is the common class for 3D segments used in the muon spectrometer.
const Amg::Vector3D & globalDirection() const
global direction
virtual const Amg::Vector3D & globalPosition() const override final
global position
virtual const Trk::PlaneSurface & associatedSurface() const override final
returns the surface for the local to global transformation
Base class for all CompetingRIOsOnTack implementations, extends the common MeasurementBase.
virtual const RIO_OnTrack & rioOnTrack(unsigned int) const =0
returns the RIO_OnTrack (also known as ROT) objects depending on the integer.
magnetic field properties to steer the behavior of the extrapolation
virtual const S & associatedSurface() const override final
Access to the Surface method.
Class to handle RIO On Tracks ROT) for InDet and Muons, it inherits from the common MeasurementBase.
Identifier identify() const
return the identifier -extends MeasurementBase
virtual const Amg::Vector3D & globalPosition() const override=0
Interface method to get the global Position.
const MeasurementBase * measurement(unsigned int) const
returns the Trk::MeasurementBase objects depending on the integer
unsigned int numberOfMeasurementBases() const
Return the number of contained Trk::MeasurementBase (s)
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
@ FullField
Field is set to be realistic, but within a given Volume.
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane