171 {
173 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
176 return {};
177 }
178 MagField::AtlasFieldCache fieldCache;
180
181 const double bfield1 =
getBField (mytracks.getFirstPerigee(), fieldCache);
182 const double bfield2 =
getBField (mytracks.getSecondPerigee(), fieldCache);
183
184
185
186 const double phitanpoca1=getphipoca(mytracks.getFirstPerigee());
187 const double phitanpoca2=getphipoca(mytracks.getSecondPerigee());
188
189
190 double abs1 = getRadiusOfCurvature(mytracks.getFirstPerigee(), bfield1);
191 double abs2 = getRadiusOfCurvature(mytracks.getSecondPerigee(),bfield2);
192
193 const double sgnr1 = abs1>0 ? 1. : -1.;
194 const double sgnr2 = abs2>0 ? 1. : -1.;
195
196
197 const std::pair<Amg::Vector3D,Amg::Vector3D> centersofcurv
198 (getCenterOfCurvature(mytracks.getFirstPerigee(),abs1,phitanpoca1),
199 getCenterOfCurvature(mytracks.getSecondPerigee(),abs2,phitanpoca2));
200
201#ifdef TRK2DDISTANCESEEDER_DEBUG
202 ATH_MSG_DEBUG(
" Center of track number 1: " << centersofcurv.first <<
" center of track 2 " << centersofcurv.second );
203#endif
204
205 const double distance2d = dist2d(centersofcurv.first,centersofcurv.second);
206
207 abs1 = std::abs(abs1);
208 abs2 = std::abs(abs2);
209
210
211#ifdef TRK2DDISTANCESEEDER_DEBUG
213 ATH_MSG_DEBUG(
"phitanpoca1: " << phitanpoca1 <<
" phitanpoca2: " << phitanpoca2 );
215#endif
216
217 double phi1 = 0;
218 double phi2 = 0;
219
220
221
222 double phi1poca = 0;
223 double phi2poca = 0;
224
225
226 if (distance2d > (abs1+abs2)) {
227
228#ifdef TRK2DDISTANCESEEDER_DEBUG
230#endif
231
232 phi1 = safeAtan2(centersofcurv.second.y()-centersofcurv.first.y(),
233 centersofcurv.second.x()-centersofcurv.first.x());
234
235 if (sgnr2<0) {
236 phi2=phi1;
237 } else {
238 phi2=oppositeangle(phi1);
239 }
240
241 if (sgnr1<0) {
242 phi1=oppositeangle(phi1);
243 }
244
245
246 phi1=takenearest(phi1,phitanpoca1);
247 phi2=takenearest(phi2,phitanpoca2);
248 }
249 else {
250
251
252
253
254 if (fabs(abs2-abs1)>=distance2d) {
255
256#ifdef TRK2DDISTANCESEEDER_DEBUG
258#endif
259
260
261 if (abs2<abs1) {
262
263#ifdef TRK2DDISTANCESEEDER_DEBUG
265#endif
266 phi1 = safeAtan2(centersofcurv.second.y()-centersofcurv.first.y(),
267 centersofcurv.second.x()-centersofcurv.first.x());
268 phi2=phi1;
269 }
270 else {
271
272#ifdef TRK2DDISTANCESEEDER_DEBUG
274#endif
275
276
277 phi1 = safeAtan2(-(centersofcurv.second.y()-centersofcurv.first.y()),
278 -(centersofcurv.second.x()-centersofcurv.first.x()));
279 phi2=phi1;
280 }
281
282
283 if (sgnr2<0) {
284 phi2=oppositeangle(phi2);
285 }
286
287 if (sgnr1<0) {
288 phi1=oppositeangle(phi1);
289 }
290
291 phi1=takenearest(setphipoca(phi1),setphipoca(phitanpoca1));
292 phi2=takenearest(setphipoca(phi2),setphipoca(phitanpoca2));
293 phi1poca = setphipoca(phi1);
294 phi2poca = setphipoca(phi2);
295 }
296 else {
297
298
299 const double projection2=(square(abs1)-square(abs2)-square(distance2d))/2./distance2d;
300 const double cosinus2=projection2/sgnr2/abs2;
301
302
303#ifdef TRK2DDISTANCESEEDER_DEBUG
304 const double projection1=(square(abs2)-square(abs1)-square(distance2d))/2./distance2d;
305 const double cosinus1=projection1/sgnr1/abs1;
306 ATH_MSG_DEBUG(
"projection1: " << projection1 <<
" cosinus1 " << cosinus1
307 << "projection2: " << projection2 << " cosinus2 " << cosinus2 );
308#endif
309
310 const double phibase2 =
311 safeAtan2(centersofcurv.second.y() - centersofcurv.first.y(),
312 centersofcurv.second.x() - centersofcurv.first.x());
313
314 const double addtophi2=
315 safeAcos(cosinus2);
316
317#ifdef TRK2DDISTANCESEEDER_DEBUG
318 ATH_MSG_DEBUG(
" phibase2 is : " << phibase2 <<
" add to phi " << addtophi2 );
319#endif
320
321 const std::pair<double,double> possiblephiontrack2(phibase2+addtophi2,phibase2-addtophi2);
322
323#ifdef TRK2DDISTANCESEEDER_DEBUG
324 ATH_MSG_DEBUG(
"the two phis are: " << possiblephiontrack2.first <<
" and " << possiblephiontrack2.second );
325#endif
326
327 const std::pair<double,double>
328 possiblecosphitrack1((centersofcurv.second.x()-centersofcurv.first.x()
329 +sgnr2*abs2*std::cos(possiblephiontrack2.first))/sgnr1/abs1,
330 (centersofcurv.second.x()-centersofcurv.first.x()
331 +sgnr2*abs2*std::cos(possiblephiontrack2.second))/sgnr1/abs1);
332
333 const std::pair<double,double>
334 possiblesigntrack1(sgnr1*(centersofcurv.second.y()-centersofcurv.first.y()
335 +sgnr2*abs2*std::sin(possiblephiontrack2.first))>0?1.:-1.,
336 sgnr1*(centersofcurv.second.y()-centersofcurv.first.y()
337 +sgnr2*abs2*std::sin(possiblephiontrack2.second))>0?1.:-1.);
338
339 const std::pair<double,double>
340 possiblephiontrack1(possiblesigntrack1.first*safeAcos(possiblecosphitrack1.first),
341 possiblesigntrack1.second*safeAcos(possiblecosphitrack1.second));
342
343
344 const std::pair<PointOnTrack,PointOnTrack>
345 possiblepointsontrack1(PointOnTrack(mytracks.getFirstPerigee(),
346 takenearest(setphipoca(possiblephiontrack1.first),
347 setphipoca(phitanpoca1))),
348 PointOnTrack(mytracks.getFirstPerigee(),
349 takenearest(setphipoca(possiblephiontrack1.second),
350 setphipoca(phitanpoca1))));
351
352 const std::pair<PointOnTrack,PointOnTrack>
353 possiblepointsontrack2(PointOnTrack(mytracks.getSecondPerigee(),
354 takenearest(setphipoca(possiblephiontrack2.first),
355 setphipoca(phitanpoca2))),
356 PointOnTrack(mytracks.getSecondPerigee(),
357 takenearest(setphipoca(possiblephiontrack2.second),
358 setphipoca(phitanpoca2))));
359
360
361 const std::pair<Amg::Vector3D,Amg::Vector3D>
362 possiblepoints1(getSeedPoint(possiblepointsontrack1.first.getPerigee(),centersofcurv.first,
363 abs1*sgnr1,possiblepointsontrack1.first.getPhiPoint()),
364 getSeedPoint(possiblepointsontrack1.second.getPerigee(),centersofcurv.first,
365 abs1*sgnr1,possiblepointsontrack1.second.getPhiPoint()));
366
367 const std::pair<Amg::Vector3D,Amg::Vector3D>
368 possiblepoints2(getSeedPoint(possiblepointsontrack2.first.getPerigee(),centersofcurv.second,
369 abs2*sgnr2,possiblepointsontrack2.first.getPhiPoint()),
370 getSeedPoint(possiblepointsontrack2.second.getPerigee(),centersofcurv.second,
371 abs2*sgnr2,possiblepointsontrack2.second.getPhiPoint()));
372
373
374#ifdef TRK2DDISTANCESEEDER_DEBUG
375 ATH_MSG_DEBUG(
"Point 1a: x " << possiblepoints1.first.x() <<
" y " << possiblepoints1.first.y() );
376 ATH_MSG_DEBUG(
"Point 2a: x " << possiblepoints2.first.x() <<
" y " << possiblepoints2.first.y() );
377 ATH_MSG_DEBUG(
"Point 1b: x " << possiblepoints1.second.x() <<
" y " << possiblepoints1.second.y() );
378 ATH_MSG_DEBUG(
"Point 2b: x " << possiblepoints2.second.x() <<
" y " << possiblepoints2.second.y() );
379 ATH_MSG_DEBUG(
"Distance between 1a and 2a:" << fabs(possiblepoints1.first.z()-possiblepoints2.first.z()) <<
380 "Distance between 1a and 2a:" << fabs(possiblepoints1.second.z()-possiblepoints2.second.z()) );
381#endif
382
383
384
386 {
387 if (fabs(possiblepoints1.first.z()-possiblepoints2.first.z())<
388 fabs(possiblepoints1.second.z()-possiblepoints2.second.z()))
389 {
390 phi1poca = possiblepointsontrack1.first.getPhiPoint();
391 phi2poca = possiblepointsontrack2.first.getPhiPoint();
392 } else {
393 phi1poca = possiblepointsontrack1.second.getPhiPoint();
394 phi2poca = possiblepointsontrack2.second.getPhiPoint();
395 }
396 }
397 else
398 {
399 if (sqrt(possiblepoints1.first.x()*possiblepoints1.first.x()+possiblepoints1.first.y()*possiblepoints1.first.y())<
400 sqrt(possiblepoints2.first.x()*possiblepoints2.first.x()+possiblepoints2.first.y()*possiblepoints2.first.y()))
401 {
402 phi1poca = possiblepointsontrack1.first.getPhiPoint();
403 phi2poca = possiblepointsontrack2.first.getPhiPoint();
404 } else {
405 phi1poca = possiblepointsontrack1.second.getPhiPoint();
406 phi2poca = possiblepointsontrack2.second.getPhiPoint();
407 }
408 }
409 phi1=goFromPhipocaToPhi(phi1poca);
410 phi2=goFromPhipocaToPhi(phi2poca);
411 }
412 }
413
414 if (twopoints) {
415 *twopoints =
TwoPoints (getSeedPoint(centersofcurv.first,
416 abs1*sgnr1,setphipoca(phi1)),
417 getSeedPoint(centersofcurv.second,
418 abs2*sgnr2,setphipoca(phi2)));
419 }
420
421 return std::make_pair (PointOnTrack (mytracks.getFirstPerigee(), phi1poca),
422 PointOnTrack (mytracks.getSecondPerigee(),phi2poca));
423 }
#define ATH_MSG_ERROR(x,...)
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
double getBField(const Perigee &p, MagField::AtlasFieldCache &cache) const
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCacheCondObjInputKey
std::pair< Amg::Vector3D, Amg::Vector3D > TwoPoints