ATLAS Offline Software
Loading...
Searching...
No Matches
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
10namespace {
11
12double 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
32namespace 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
42double 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 */
vkalMagCnvCst
#define vkalMagCnvCst
Definition CommonPars.h:23
diff
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
Definition Jet.cxx:631
angle
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
Definition TRTDetectorFactory_Full.cxx:71
cross_ref
#define cross_ref(a_1, a_2)
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::vkPerigeeToP
void vkPerigeeToP(const double *perig3, double *pp, double effectiveBMAG)
Definition cfMomentum.cxx:15
Trk::vkvFastV
double vkvFastV(double *p1, double *p2, const double *vRef, double dbmag, double *out)
Definition VKvFast.cxx:42
Generated on for ATLAS Offline Software by doxygen 1.14.0