Estimate signed momentum for two segments by fitting 2 segments to one approximate track model
105{
106
109
110 ATH_MSG_VERBOSE(
" Executing MuonSegmentMomentumTool fitMomentum2Segments ");
111
112 const Muon::MuonSegment* myseg1=segment1,*myseg2=segment2;
114 myseg1=segment2;
115 myseg2=segment1;
116 }
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;
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;
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));
139 }
140 if (!rot) continue;
142
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));
154 }
155 if (!rot) continue;
159 if (!firstphi2) firstphi2=rot;
160 lastphi2=rot;
161 }
162 }
163 bool flip = false;
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
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);
189 double resi[4],qoverp[4];
190 for (
int i=0;
i<4;
i++){
191 std::optional<Trk::TransportJacobian> jac{};
195 ctx, startpar,worstseg->associatedSurface(),
199 if (par && jac && (*jac)(1,4)!=0){
202 qoverp[
i] = 1/signedMomentum;
203
204 double delta_qoverp=
residual/(*jac)(1,4);
205 double der_simple = -10.*(bestseg->
globalPosition()-worstseg->globalPosition()).mag()/(.3*fieldintegral);
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));
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));
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
227 signedMomentum = signedMomentum_updated>0? 1000.:-1000;
228 } else {
229
230 if(i<3) signedMomentum = signedMomentum_updated;
231 }
232
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}
#define ATH_MSG_VERBOSE(x)
double fieldIntegralEstimate(const EventContext &ctx, const Muon::MuonSegment *segment1, const Muon::MuonSegment *segment2) const
ToolHandle< Trk::IPropagator > m_propagator
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
virtual const Trk::PlaneSurface & associatedSurface() const override final
returns the surface for the local to global transformation
virtual const RIO_OnTrack & rioOnTrack(unsigned int) const =0
returns the RIO_OnTrack (also known as ROT) objects depending on the integer.
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)
constexpr double derivative(const double x)
Evaluates the d-th derivative of the n-th Legendre polynomial at x.
@ FullField
Field is set to be realistic, but within a given Volume.
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane