ATLAS Offline Software
VKvFast.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "TrkVKalVrtCore/VKvFast.h"
6 #include "TrkVKalVrtCore/CommonPars.h"
7 #include "TrkVKalVrtCore/cfMomentum.h"
8 #include <cmath>
9 
10 namespace {
11 
12 double vkvang_(double *crs, double *centr, const double *r__, const double *xd,
13  const double *yd) {
14  double ret_val, cosf, sinf, cosp, sinp;
15 
16  /* Parameter adjustments */
17  --centr;
18  --crs;
19 
20  /* Function Body */
21  sinf = (crs[1] - centr[1]) / *r__;
22  cosf = -(crs[2] - centr[2]) / *r__;
23  sinp = *xd * sinf - *yd * cosf;
24  /* vector prod. */
25  cosp = *xd * cosf + *yd * sinf;
26  /* scalar prod. */
27  ret_val = atan2(sinp, cosp);
28  return ret_val;
29 }
30 } // namespace
31 
32 namespace Trk {
33 
34 /* ================================================================== */
35 /* Fast estimation of the crossing point of 2 charged tracks */
36 /* p1(5),p2(5) - perigee parameters of tracks */
37 /* dBMAG - magnetic field */
38 /* out(3) - vertex position estimation */
39 /* Author: V.Kostioukhine */
40 /* ================================================================== */
41 
42 double vkvFastV(double *p1, double *p2, const double* vRef, double dbmag, double *out)
43 {
44  double d__1;
45  double diff, coef, cent1[2], cent2[2], angle, cross[6];/* was [3][2] */
46  double r1, r2, z1, z2, a01[3], a02[3], dc, an[2];
47  double ar1, ar2, xd1, xd2, yd1, yd2, dz1, dz2, pt1, pt2, pz1, pz2, det;
48  double loc_bmag;
49 
50 
51 #define cross_ref(a_1,a_2) cross[(a_2)*3 + (a_1) - 4]
52 
53 /* ---- */
54  /* Parameter adjustments */
55  --out;
56  --p2;
57  --p1;
58 
59  /* Function Body */
60  loc_bmag = dbmag;
61  if (dbmag <= .1) {loc_bmag = 0.1;} // To avoid ZERO magnetic field
62 /* ---- */
63  double psum[3]={0.,0.,0.}; double ptmp[3]={0.,0.,0.};
64  vkPerigeeToP(&p1[3],ptmp,loc_bmag); psum[0]+=ptmp[0];psum[1]+=ptmp[1];psum[2]+=ptmp[2];
65  vkPerigeeToP(&p2[3],ptmp,loc_bmag); psum[0]+=ptmp[0];psum[1]+=ptmp[1];psum[2]+=ptmp[2];
66 /*------*/
67  r1 = 1. / p1[5];
68  r2 = 1. / p2[5];
69  ar1 = std::abs(r1);
70  ar2 = std::abs(r2);
71 // pt1 = ar1 * loc_bmag * .00299979; // For GeV and CM
72 // pt2 = ar2 * loc_bmag * .00299979;
73 // pt1 = ar1 * loc_bmag * .299979; // For MeV and MM
74 // pt2 = ar2 * loc_bmag * .299979;
75  pt1 = ar1 * loc_bmag * vkalMagCnvCst; // from CommonPars.h
76  pt2 = ar2 * loc_bmag * vkalMagCnvCst;
77  pz1 = pt1 / tan(p1[3]);
78  pz2 = pt2 / tan(p2[3]);
79 /* --- */
80  xd1 = cos(p1[4]);
81 /* Need for calculation of Z of crossing point */
82  xd2 = cos(p2[4]);
83  yd1 = sin(p1[4]);
84  yd2 = sin(p2[4]);
85 /* --- Point of minimal approach to the beam */
86  a01[0] = p1[1] * yd1;
87  a01[1] = -p1[1] * xd1;
88  a01[2] = p1[2];
89  a02[0] = p2[1] * yd2;
90  a02[1] = -p2[1] * xd2;
91  a02[2] = p2[2];
92 /* -- coordinates of centres of trajectories */
93  cent1[0] = a01[0] - r1 * yd1;
94  cent1[1] = a01[1] + r1 * xd1;
95  cent2[0] = a02[0] - r2 * yd2;
96  cent2[1] = a02[1] + r2 * xd2;
97 /* -- */
98  an[0] = cent2[0] - cent1[0];
99  an[1] = cent2[1] - cent1[1];
100  dc = sqrt(an[0]*an[0] + an[1]*an[1]);
101  an[0] /= dc;
102  an[1] /= dc;
103 /* -- Single point*/
104  if (dc <= std::abs(ar1 - ar2)) {
105  if (ar1 < ar2) {
106  cross_ref(1, 1) = ar1 / dc * (cent1[0] - cent2[0]) + cent1[0];
107  cross_ref(1, 2) = ar2 / dc * (cent1[0] - cent2[0]) + cent2[0];
108  cross_ref(2, 1) = ar1 / dc * (cent1[1] - cent2[1]) + cent1[1];
109  cross_ref(2, 2) = ar2 / dc * (cent1[1] - cent2[1]) + cent2[1];
110  } else {
111  cross_ref(1, 1) = ar1 / dc * (cent2[0] - cent1[0]) + cent1[0];
112  cross_ref(1, 2) = ar2 / dc * (cent2[0] - cent1[0]) + cent2[0];
113  cross_ref(2, 1) = ar1 / dc * (cent2[1] - cent1[1]) + cent1[1];
114  cross_ref(2, 2) = ar2 / dc * (cent2[1] - cent1[1]) + cent2[1];
115  }
116  angle = vkvang_(&cross_ref(1, 1), cent1, &r1, &xd1, &yd1);
117  z1 = pz1 * r1 * angle / pt1 + a01[2];
118  angle = vkvang_(&cross_ref(1, 2), cent2, &r2, &xd2, &yd2);
119  z2 = pz2 * r2 * angle / pt2 + a02[2];
120  //diffz = std::abs(z2 - z1);
121  out[1] = (cross_ref(1, 1) + cross_ref(1, 2)) / 2.;
122  out[2] = (cross_ref(2, 1) + cross_ref(2, 2)) / 2.;
123  out[3] = (z1 + z2) / 2.;
124  return std::abs(z1 - z2);
125  }
126 /* ---------------------------------------------------------- */
127  double bestZdiff=0.;
128  diff = dc - (ar1 + ar2);
129  if (diff < 0.) {
130 /* 2 distinct points, A1 & A2. */
131  coef = (r1*r1 + dc*dc - r2*r2) / (dc * 2);
132  if (r1*r1 - coef*coef > 0.) {
133  det = sqrt(r1 * r1 - coef * coef);
134  } else {
135  det = 0.;
136  }
137  cross_ref(1, 1) = cent1[0] + coef * an[0] + det * an[1];
138  cross_ref(2, 1) = cent1[1] + coef * an[1] - det * an[0];
139  cross_ref(3, 1) = 0.;
140  cross_ref(1, 2) = cent1[0] + coef * an[0] - det * an[1];
141  cross_ref(2, 2) = cent1[1] + coef * an[1] + det * an[0];
142  cross_ref(3, 2) = 0.;
143 /* -First cross */
144  angle = vkvang_(&cross_ref(1, 1), cent1, &r1, &xd1, &yd1);
145  z1 = pz1 * r1 * angle / pt1 + a01[2];
146  angle = vkvang_(&cross_ref(1, 1), cent2, &r2, &xd2, &yd2);
147  z2 = pz2 * r2 * angle / pt2 + a02[2];
148  dz1 = std::abs(z2 - z1);
149  cross_ref(3, 1) = (z1 + z2) / 2.;
150 /* -Second cross */
151  angle = vkvang_(&cross_ref(1, 2), cent1, &r1, &xd1, &yd1);
152  z1 = pz1 * r1 * angle / pt1 + a01[2];
153  angle = vkvang_(&cross_ref(1, 2), cent2, &r2, &xd2, &yd2);
154  z2 = pz2 * r2 * angle / pt2 + a02[2];
155  dz2 = std::abs(z2 - z1);
156  cross_ref(3, 2) = (z1 + z2) / 2.;
157 /* if reterence vertex is present -> check 2D mom direction -> must be avay from it*/
158  if(vRef && std::abs(dz1-dz2)<5.){
159  double dir1=ptmp[0]*(cross_ref(1,1)-vRef[0])+ptmp[1]*(cross_ref(2,1)-vRef[1]);
160  double dir2=ptmp[0]*(cross_ref(1,2)-vRef[0])+ptmp[1]*(cross_ref(2,2)-vRef[1]);
161  if(dir1*dir2<0){ // points are on different sides of ref. vertex
162  if(dir1<0) dz1+=1000.;
163  if(dir2<0) dz2+=1000.;
164  }
165  }
166 /* - choice of best dz */
167 
168  if (dz1 < dz2) {
169  out[1] = cross_ref(1, 1);
170  out[2] = cross_ref(2, 1);
171  out[3] = cross_ref(3, 1);
172  bestZdiff=dz1;
173  } else {
174  out[1] = cross_ref(1, 2);
175  out[2] = cross_ref(2, 2);
176  out[3] = cross_ref(3, 2);
177  bestZdiff=dz2;
178  }
179  } else {
180  cross_ref(1, 1) = (ar1 * cent2[0] + ar2 * cent1[0]) / (ar1 + ar2);
181  cross_ref(2, 1) = (ar1 * cent2[1] + ar2 * cent1[1]) / (ar1 + ar2);
182  out[1] = cross_ref(1, 1);
183  out[2] = cross_ref(2, 1);
184  out[3] = 0.;
185  d__1 = r1 * dc / (ar1 + ar2);
186  angle = vkvang_(cross, cent1, &d__1, &xd1, &yd1);
187  z1 = pz1 * r1 * angle / pt1 + a01[2];
188  d__1 = r2 * dc / (ar1 + ar2);
189  angle = vkvang_(cross, cent2, &d__1, &xd2, &yd2);
190  z2 = pz2 * r2 * angle / pt2 + a02[2];
191  out[3] = (z1 + z2) / 2.;
192  bestZdiff = std::abs(z1 - z2);
193  }
194  return bestZdiff;
195 }
196 #undef cross_ref
197 
198 
199 
200 } /* End of VKalVrtCore namespace */
CommonPars.h
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
VKvFast.h
vkalMagCnvCst
#define vkalMagCnvCst
Definition: CommonPars.h:23
angle
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
Definition: TRTDetectorFactory_Full.cxx:73
WritePulseShapeToCool.det
det
Definition: WritePulseShapeToCool.py:204
Trk::vkvFastV
double vkvFastV(double *p1, double *p2, const double *vRef, double dbmag, double *out)
Definition: VKvFast.cxx:42
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
plotBeamSpotCompare.xd
xd
Definition: plotBeamSpotCompare.py:220
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
cross_ref
#define cross_ref(a_1, a_2)
Trk::vkPerigeeToP
void vkPerigeeToP(const double *perig3, double *pp, double BMAG)
Definition: cfMomentum.cxx:15
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
WriteCalibToCool.coef
coef
Definition: WriteCalibToCool.py:582
cfMomentum.h
Generated on Fri Jun 28 2024 21:33:10 for ATLAS Offline Software by   doxygen 1.8.18