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