ATLAS Offline Software
Loading...
Searching...
No Matches
MuonSegmentMomentumFromField.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3*/
4
7#include <sstream>
8#include <iostream>
9#include <vector>
11
18#include <TString.h> // for Form
19
20MuonSegmentMomentumFromField::MuonSegmentMomentumFromField(const std::string& type,const std::string& name,const IInterface*
21 parent):AthAlgTool(type,name,parent)
22{
23 declareInterface<IMuonSegmentMomentumEstimator>(this);
24}
25
27{
28 ATH_MSG_VERBOSE(" MuonSegmentMomentumFromField::Initializing ");
29 ATH_CHECK(AthAlgTool::initialize());
30 ATH_CHECK(m_propagator.retrieve());
31 ATH_CHECK(m_navigator.retrieve());
32 ATH_CHECK(m_idHelperSvc.retrieve());
34 ATH_MSG_VERBOSE("End of Initializing");
35 return StatusCode::SUCCESS;
36}
37
38void MuonSegmentMomentumFromField::fitMomentumVectorSegments( const EventContext& ctx, const std::vector <const Muon::MuonSegment*> & segments, double & signedMomentum ) const
39{
40
43
44 ATH_MSG_VERBOSE(" Executing MuonSegmentMomentumTool fitMomentumVectorSegments ");
45 ATH_MSG_DEBUG(" fitMomentumVectorSegments " << segments.size() << " segments " );
46
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();
50
51 ++it2;
52 double maxintegral=0;
53 while (it2!=it_end){
54 double integral= m_doOld ? fieldIntegralEstimate_old(ctx,*it,*it2) : fieldIntegralEstimate(ctx,*it,*it2);
55 if (std::abs(integral)>maxintegral){
56 maxintegral=std::abs(integral);
57 if( m_doOld ) fitMomentum2Segments_old(ctx, *it, *it2, signedMomentum);
58 else fitMomentum2Segments(ctx, *it, *it2, signedMomentum);
59 }
60 ++it;
61 ++it2;
62 }
63 ATH_MSG_DEBUG( " Estimated signed momentum " << signedMomentum );
64}
65
66double MuonSegmentMomentumFromField::fieldIntegralEstimate( const EventContext& ctx, const Muon::MuonSegment* segment1, const Muon::MuonSegment* segment2) const
67{
68
69 Amg::Vector3D pos1=segment1->globalPosition();
70 Amg::Vector3D pos2=segment2->globalPosition();
71 Amg::Vector3D posdiff=pos2-pos1;
72 if (posdiff.dot(segment1->globalDirection())<0){
73 posdiff=-posdiff;
74 pos1=segment2->globalPosition();
75 pos2=segment1->globalPosition();
76 }
77 Amg::Vector3D point1=pos1+.25*posdiff;
78 Amg::Vector3D point2=pos1+.5*posdiff;
79 Amg::Vector3D point3=pos1+.75*posdiff;
80 Amg::Vector3D field1,field2,field3;
81
84 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
85 if (!fieldCondObj) {
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()));
87 }
88 fieldCondObj->getInitializedCache(fieldCache);
89
90 fieldCache.getField(point1.data(),field1.data()); // field in kiloTesla
91 fieldCache.getField(point2.data(),field2.data());
92 fieldCache.getField(point3.data(),field3.data());
93 ATH_MSG_DEBUG("Mid Point " << Amg::toString(point2) << " field " << Amg::toString(field2) );
94 field1[2]=field2[2]=field3[2]=0;
95 Amg::Vector3D fieldsum=field1+field2+field3;
96 Amg::Vector3D crossvec=posdiff.unit().cross(field1+field2+field3);
97 Amg::Vector2D rphidir(-posdiff.y(),posdiff.x());
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();
102}
103
104void MuonSegmentMomentumFromField::fitMomentum2Segments( const EventContext& ctx, const Muon::MuonSegment* segment1, const Muon::MuonSegment* segment2, double & signedMomentum ) const
105{
106
109
110 ATH_MSG_VERBOSE(" Executing MuonSegmentMomentumTool fitMomentum2Segments ");
111
112 const Muon::MuonSegment* myseg1=segment1,*myseg2=segment2;
113 if (myseg1->globalDirection().dot(myseg2->globalPosition()-myseg1->globalPosition())<0) {
114 myseg1=segment2;
115 myseg2=segment1;
116 }
117 double fieldintegral=fieldIntegralEstimate(ctx, myseg1,myseg2);
118 double theta1=myseg1->globalDirection().theta();
119 double theta2=myseg2->globalDirection().theta();
120 double phi1=myseg1->globalDirection().phi();
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;
126 else {
127 if (theta1>M_PI/2) deltatheta=2*M_PI-theta2-theta1;
128 else deltatheta=-theta1-theta2;
129 }
130 double dist1=-1,dist2=-1;
131 const Trk::RIO_OnTrack *firstphi1=nullptr,*lastphi1=nullptr,*firstphi2=nullptr,*lastphi2=nullptr;
132
133 const Muon::MuonSegment *bestseg=myseg1,*worstseg=myseg2;
134 for (int i=0;i<(int)myseg1->numberOfMeasurementBases();i++){
135 const Trk::RIO_OnTrack *rot=dynamic_cast<const Trk::RIO_OnTrack *>(myseg1->measurement(i));
136 if (!rot){
137 const Trk::CompetingRIOsOnTrack *crot=dynamic_cast<const Trk::CompetingRIOsOnTrack *>(myseg1->measurement(i));
138 if (crot) rot=&crot->rioOnTrack(0);
139 }
140 if (!rot) continue;
141 Identifier id=rot->identify();
142
143 if ((m_idHelperSvc->isRpc(id) && m_idHelperSvc->rpcIdHelper().measuresPhi(id)) || (m_idHelperSvc->isCsc(id) && m_idHelperSvc->cscIdHelper().measuresPhi(id)) || (m_idHelperSvc->isTgc(id) && m_idHelperSvc->tgcIdHelper().isStrip(id))
144 || (m_idHelperSvc->issTgc(id) && m_idHelperSvc->stgcIdHelper().measuresPhi(id) ) ){
145 if (!firstphi1) firstphi1=rot;
146 lastphi1=rot;
147 }
148 }
149 for (int i=0;i<(int)myseg2->numberOfMeasurementBases();i++){
150 const Trk::RIO_OnTrack *rot=dynamic_cast<const Trk::RIO_OnTrack *>(myseg2->measurement(i));
151 if (!rot){
152 const Trk::CompetingRIOsOnTrack *crot=dynamic_cast<const Trk::CompetingRIOsOnTrack *>(myseg2->measurement(i));
153 if (crot) rot=&crot->rioOnTrack(0);
154 }
155 if (!rot) continue;
156 Identifier id=rot->identify();
157 if ((m_idHelperSvc->isRpc(id) && m_idHelperSvc->rpcIdHelper().measuresPhi(id)) || (m_idHelperSvc->isCsc(id) && m_idHelperSvc->cscIdHelper().measuresPhi(id)) || (m_idHelperSvc->isTgc(id) && m_idHelperSvc->tgcIdHelper().isStrip(id))
158 || (m_idHelperSvc->issTgc(id) && m_idHelperSvc->stgcIdHelper().measuresPhi(id) ) ){
159 if (!firstphi2) firstphi2=rot;
160 lastphi2=rot;
161 }
162 }
163 bool flip = false;
164 if (firstphi1) dist1=std::abs((firstphi1->globalPosition()-lastphi1->globalPosition()).dot(myseg1->globalDirection()));
165 if (firstphi2) dist2=std::abs((firstphi2->globalPosition()-lastphi2->globalPosition()).dot(myseg2->globalDirection()));
166 if (dist2>dist1) {
167 flip = true;
168 bestseg=myseg2;
169 worstseg=myseg1;
170 }
171
172 for (int itry=0;itry<2;itry++){
173
174 if(itry==1) {
175// retry with other as best segmt segment
176 if(!flip) {
177 bestseg = myseg2;
178 worstseg=myseg1;
179 } else {
180 bestseg = myseg1;
181 worstseg=myseg2;
182 }
183 }
184
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{};
192 Trk::AtaPlane startpar(bestseg->globalPosition(),bestseg->globalDirection().phi(),
193 bestseg->globalDirection().theta(),1/signedMomentum,bestseg->associatedSurface());
194 auto par=m_propagator->propagateParameters(
195 ctx, startpar,worstseg->associatedSurface(),
196 (bestseg==myseg1) ? Trk::alongMomentum : Trk::oppositeMomentum,false,
198 ATH_MSG_DEBUG("par: " << par.get() << " jac: " << *jac );
199 if (par && jac && (*jac)(1,4)!=0){
200 residual = worstseg->localParameters()[Trk::locY] - par->parameters()[Trk::locY];
201 resi[i] = residual;
202 qoverp[i] = 1/signedMomentum;
203// update qoverp
204 double delta_qoverp=residual/(*jac)(1,4);
205 double der_simple = -10.*(bestseg->globalPosition()-worstseg->globalPosition()).mag()/(.3*fieldintegral);
206 if(bestseg->globalDirection().z()<0) der_simple = - der_simple;
207 if(i>0) {
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;
214 }
215 }
216 } else {
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;
220 }
221 }
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 );
226// protect against too low momenta as propagation will fail
227 signedMomentum = signedMomentum_updated>0? 1000.:-1000;
228 } else {
229// do not change momentum beyond last iteration
230 if(i<3) signedMomentum = signedMomentum_updated;
231 }
232 //delete par;
233 if (std::abs(residual)<1) break;
234 }
235 }
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 ");
239 } else break;
240 }
241}
242
243
244double MuonSegmentMomentumFromField::fieldIntegralEstimate_old( const EventContext& ctx, const Muon::MuonSegment* segment1, const Muon::MuonSegment* segment2) const
245{
246
247 Amg::Vector3D pos1=segment1->globalPosition();
248 Amg::Vector3D pos2=segment2->globalPosition();
249 Amg::Vector3D posdiff=pos2-pos1;
250 if (posdiff.dot(segment1->globalDirection())<0){
251 posdiff=-posdiff;
252 pos1=segment2->globalPosition();
253 pos2=segment1->globalPosition();
254 }
255 Amg::Vector3D point1=pos1+.25*posdiff;
256 Amg::Vector3D point2=pos1+.5*posdiff;
257 Amg::Vector3D point3=pos1+.75*posdiff;
258 Amg::Vector3D field1,field2,field3;
259
260 MagField::AtlasFieldCache fieldCache;
262 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
263 if (!fieldCondObj) {
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()));
265 }
266 fieldCondObj->getInitializedCache(fieldCache);
267
268 fieldCache.getField(point1.data(),field1.data()); // field in kiloTesla
269 fieldCache.getField(point2.data(),field2.data());
270 fieldCache.getField(point3.data(),field3.data());
271 ATH_MSG_DEBUG("Mid Point " << Amg::toString(point2) << " field " << Amg::toString(field2) );
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();
275}
276
277void MuonSegmentMomentumFromField::fitMomentum2Segments_old( const EventContext& ctx, const Muon::MuonSegment* segment1, const Muon::MuonSegment* segment2, double & signedMomentum ) const
278{
279
282
283 ATH_MSG_VERBOSE(" Executing MuonSegmentMomentumTool fitMomentum2Segments ");
284
285 const Muon::MuonSegment* myseg1=segment1,*myseg2=segment2;
286 if (myseg1->globalDirection().dot(myseg2->globalPosition()-myseg1->globalPosition())<0) {
287 myseg1=segment2;
288 myseg2=segment1;
289 }
290 double fieldintegral=fieldIntegralEstimate(ctx, myseg1,myseg2);
291 double theta1=myseg1->globalDirection().theta();
292 double theta2=myseg2->globalDirection().theta();
293 double phi1=myseg1->globalDirection().phi();
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;
299 else {
300 if (theta1>M_PI/2) deltatheta=2*M_PI-theta2-theta1;
301 else deltatheta=-theta1-theta2;
302 }
303 double dist1=-1,dist2=-1;
304 const Trk::RIO_OnTrack *firstphi1=nullptr,*lastphi1=nullptr,*firstphi2=nullptr,*lastphi2=nullptr;
305
306 const Muon::MuonSegment *bestseg=myseg1,*worstseg=myseg2;
307 for (int i=0;i<(int)myseg1->numberOfMeasurementBases();i++){
308 const Trk::RIO_OnTrack *rot=dynamic_cast<const Trk::RIO_OnTrack *>(myseg1->measurement(i));
309 if (!rot){
310 const Trk::CompetingRIOsOnTrack *crot=dynamic_cast<const Trk::CompetingRIOsOnTrack *>(myseg1->measurement(i));
311 if (crot) rot=&crot->rioOnTrack(0);
312 }
313 if (!rot) continue;
314 Identifier id=rot->identify();
315
316 if ((m_idHelperSvc->isRpc(id) && m_idHelperSvc->rpcIdHelper().measuresPhi(id)) || (m_idHelperSvc->isCsc(id) && m_idHelperSvc->cscIdHelper().measuresPhi(id)) ||
317 (m_idHelperSvc->isTgc(id) && m_idHelperSvc->tgcIdHelper().isStrip(id))){
318 if (!firstphi1) firstphi1=rot;
319 lastphi1=rot;
320 }
321 }
322 for (int i=0;i<(int)myseg2->numberOfMeasurementBases();i++){
323 const Trk::RIO_OnTrack *rot=dynamic_cast<const Trk::RIO_OnTrack *>(myseg2->measurement(i));
324 if (!rot){
325 const Trk::CompetingRIOsOnTrack *crot=dynamic_cast<const Trk::CompetingRIOsOnTrack *>(myseg2->measurement(i));
326 if (crot) rot=&crot->rioOnTrack(0);
327 }
328 if (!rot) continue;
329 Identifier id=rot->identify();
330 if ((m_idHelperSvc->isRpc(id) && m_idHelperSvc->rpcIdHelper().measuresPhi(id)) || (m_idHelperSvc->isCsc(id) && m_idHelperSvc->cscIdHelper().measuresPhi(id)) ||
331 (m_idHelperSvc->isTgc(id) && m_idHelperSvc->tgcIdHelper().isStrip(id))){
332 if (!firstphi2) firstphi2=rot;
333 lastphi2=rot;
334 }
335 }
336
337 if (firstphi1) dist1=std::abs((firstphi1->globalPosition()-lastphi1->globalPosition()).dot(myseg1->globalDirection()));
338 if (firstphi2) dist2=std::abs((firstphi2->globalPosition()-lastphi2->globalPosition()).dot(myseg2->globalDirection()));
339 if (dist2>dist1) {
340 bestseg=myseg2;
341 worstseg=myseg1;
342 }
343 signedMomentum =-.3e3*fieldintegral/deltatheta;
344 ATH_MSG_DEBUG("integral: " << fieldintegral << " deltatheta: " << deltatheta << " signedmomentum : " << signedMomentum);
345 for (int i=0;i<3;i++){
346 Trk::AtaPlane startpar(bestseg->globalPosition(),bestseg->globalDirection().phi(),bestseg->globalDirection().theta(),
347 1/signedMomentum,bestseg->associatedSurface());
348 std::optional<Trk::TransportJacobian> jac{};
349 auto par=m_propagator->propagateParameters(
350 ctx,startpar,worstseg->associatedSurface(),
351 (bestseg==myseg1) ? Trk::alongMomentum : Trk::oppositeMomentum,false,
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);
358 // delete par;
359 if (std::abs(residual)<1) break;
360 }
361 }
362}
#define M_PI
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_DEBUG(x)
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
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
This is the common class for 3D segments used in the muon spectrometer.
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.
Definition RIO_OnTrack.h:70
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)
double integral(TH1 *h)
Definition computils.cxx:59
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
@ oppositeMomentum
@ alongMomentum
@ FullField
Field is set to be realistic, but within a given Volume.
@ locY
local cartesian
Definition ParamDefs.h:38
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane