1 |
#include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
|
2 |
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
|
3 |
#include "MagneticField/Engine/interface/MagneticField.h"
|
4 |
#include "FWCore/Utilities/interface/Exception.h"
|
5 |
|
6 |
using namespace std;
|
7 |
|
8 |
bool ClosestApproachInRPhi::calculate(const TrajectoryStateOnSurface & sta,
|
9 |
const TrajectoryStateOnSurface & stb)
|
10 |
{
|
11 |
TrackCharge chargeA = sta.charge(); TrackCharge chargeB = stb.charge();
|
12 |
GlobalVector momentumA = sta.globalMomentum();
|
13 |
GlobalVector momentumB = stb.globalMomentum();
|
14 |
GlobalPoint positionA = sta.globalPosition();
|
15 |
GlobalPoint positionB = stb.globalPosition();
|
16 |
paramA = sta.globalParameters();
|
17 |
paramB = stb.globalParameters();
|
18 |
// compute magnetic field ONCE
|
19 |
bz = sta.freeState()->parameters().magneticField().inTesla(positionA).z() * 2.99792458e-3;
|
20 |
|
21 |
return compute(chargeA, momentumA, positionA, chargeB, momentumB, positionB);
|
22 |
|
23 |
}
|
24 |
|
25 |
|
26 |
bool ClosestApproachInRPhi::calculate(const FreeTrajectoryState & sta,
|
27 |
const FreeTrajectoryState & stb)
|
28 |
{
|
29 |
TrackCharge chargeA = sta.charge(); TrackCharge chargeB = stb.charge();
|
30 |
GlobalVector momentumA = sta.momentum();
|
31 |
GlobalVector momentumB = stb.momentum();
|
32 |
GlobalPoint positionA = sta.position();
|
33 |
GlobalPoint positionB = stb.position();
|
34 |
paramA = sta.parameters();
|
35 |
paramB = stb.parameters();
|
36 |
// compute magnetic field ONCE
|
37 |
bz = sta.parameters().magneticField().inTesla(positionA).z() * 2.99792458e-3;
|
38 |
|
39 |
return compute(chargeA, momentumA, positionA, chargeB, momentumB, positionB);
|
40 |
|
41 |
}
|
42 |
|
43 |
pair<GlobalPoint, GlobalPoint> ClosestApproachInRPhi::points() const
|
44 |
{
|
45 |
if (!status_)
|
46 |
throw cms::Exception("TrackingTools/PatternTools","ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
|
47 |
return pair<GlobalPoint, GlobalPoint> (posA, posB);
|
48 |
}
|
49 |
|
50 |
|
51 |
GlobalPoint
|
52 |
ClosestApproachInRPhi::crossingPoint() const
|
53 |
{
|
54 |
if (!status_)
|
55 |
throw cms::Exception("TrackingTools/PatternTools","ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
|
56 |
return GlobalPoint(0.5*(posA.basicVector() + posB.basicVector()));
|
57 |
|
58 |
}
|
59 |
|
60 |
|
61 |
float ClosestApproachInRPhi::distance() const
|
62 |
{
|
63 |
if (!status_)
|
64 |
throw cms::Exception("TrackingTools/PatternTools","ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
|
65 |
return (posB - posA).mag();
|
66 |
}
|
67 |
|
68 |
|
69 |
bool ClosestApproachInRPhi::compute(const TrackCharge & chargeA,
|
70 |
const GlobalVector & momentumA,
|
71 |
const GlobalPoint & positionA,
|
72 |
const TrackCharge & chargeB,
|
73 |
const GlobalVector & momentumB,
|
74 |
const GlobalPoint & positionB)
|
75 |
{
|
76 |
|
77 |
|
78 |
// centres and radii of track circles
|
79 |
double xca, yca, ra;
|
80 |
circleParameters(chargeA, momentumA, positionA, xca, yca, ra, bz);
|
81 |
double xcb, ycb, rb;
|
82 |
circleParameters(chargeB, momentumB, positionB, xcb, ycb, rb, bz);
|
83 |
|
84 |
// points of closest approach in transverse plane
|
85 |
double xg1, yg1, xg2, yg2;
|
86 |
int flag = transverseCoord(xca, yca, ra, xcb, ycb, rb, xg1, yg1, xg2, yg2);
|
87 |
if (flag == 0) {
|
88 |
status_ = false;
|
89 |
return false;
|
90 |
}
|
91 |
|
92 |
double xga, yga, zga, xgb, ygb, zgb;
|
93 |
|
94 |
if (flag == 1) {
|
95 |
// two crossing points on each track in transverse plane
|
96 |
// select point for which z-coordinates on the 2 tracks are the closest
|
97 |
double za1 = zCoord(momentumA, positionA, ra, xca, yca, xg1, yg1);
|
98 |
double zb1 = zCoord(momentumB, positionB, rb, xcb, ycb, xg1, yg1);
|
99 |
double za2 = zCoord(momentumA, positionA, ra, xca, yca, xg2, yg2);
|
100 |
double zb2 = zCoord(momentumB, positionB, rb, xcb, ycb, xg2, yg2);
|
101 |
|
102 |
if (abs(zb1 - za1) < abs(zb2 - za2)) {
|
103 |
xga = xg1; yga = yg1; zga = za1; zgb = zb1;
|
104 |
}
|
105 |
else {
|
106 |
xga = xg2; yga = yg2; zga = za2; zgb = zb2;
|
107 |
}
|
108 |
xgb = xga; ygb = yga;
|
109 |
}
|
110 |
else {
|
111 |
// one point of closest approach on each track in transverse plane
|
112 |
xga = xg1; yga = yg1;
|
113 |
zga = zCoord(momentumA, positionA, ra, xca, yca, xga, yga);
|
114 |
xgb = xg2; ygb = yg2;
|
115 |
zgb = zCoord(momentumB, positionB, rb, xcb, ycb, xgb, ygb);
|
116 |
}
|
117 |
|
118 |
posA = GlobalPoint(xga, yga, zga);
|
119 |
posB = GlobalPoint(xgb, ygb, zgb);
|
120 |
status_ = true;
|
121 |
return true;
|
122 |
}
|
123 |
|
124 |
pair <GlobalTrajectoryParameters, GlobalTrajectoryParameters>
|
125 |
ClosestApproachInRPhi::trajectoryParameters () const
|
126 |
{
|
127 |
if (!status_)
|
128 |
throw cms::Exception("TrackingTools/PatternTools","ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
|
129 |
pair <GlobalTrajectoryParameters, GlobalTrajectoryParameters>
|
130 |
ret ( newTrajectory( posA, paramA, bz),
|
131 |
newTrajectory( posB, paramB, bz) );
|
132 |
return ret;
|
133 |
}
|
134 |
|
135 |
GlobalTrajectoryParameters
|
136 |
ClosestApproachInRPhi::newTrajectory( const GlobalPoint & newpt, const GlobalTrajectoryParameters & oldgtp, double bz )
|
137 |
{
|
138 |
// First we need the centers of the circles.
|
139 |
double qob = oldgtp.charge()/bz;
|
140 |
double xc = oldgtp.position().x() + qob * oldgtp.momentum().y();
|
141 |
double yc = oldgtp.position().y() - qob * oldgtp.momentum().x();
|
142 |
|
143 |
// and of course....
|
144 |
double npx = (newpt.y()-yc)*(bz/oldgtp.charge());
|
145 |
double npy = (xc-newpt.x())*(bz/oldgtp.charge());
|
146 |
|
147 |
/*
|
148 |
* old code: slow and wrong
|
149 |
*
|
150 |
// now we do a translation, move the center of circle to (0,0,0).
|
151 |
double dx1 = oldgtp.position().x() - xc;
|
152 |
double dy1 = oldgtp.position().y() - yc;
|
153 |
double dx2 = newpt.x() - xc;
|
154 |
double dy2 = newpt.y() - yc;
|
155 |
|
156 |
// now for the angles:
|
157 |
double cosphi = ( dx1 * dx2 + dy1 * dy2 ) /
|
158 |
( sqrt ( dx1 * dx1 + dy1 * dy1 ) * sqrt ( dx2 * dx2 + dy2 * dy2 ));
|
159 |
double sinphi = - oldgtp.charge() * sqrt ( 1 - cosphi * cosphi );
|
160 |
|
161 |
// Finally, the new momenta:
|
162 |
double px = cosphi * oldgtp.momentum().x() - sinphi * oldgtp.momentum().y();
|
163 |
double py = sinphi * oldgtp.momentum().x() + cosphi * oldgtp.momentum().y();
|
164 |
|
165 |
std::cout << px-npx << " " << py-npy << ", " << oldgtp.charge() << std::endl;
|
166 |
*/
|
167 |
|
168 |
GlobalVector vta ( npx, npy, oldgtp.momentum().z() );
|
169 |
GlobalTrajectoryParameters gta( newpt , vta , oldgtp.charge(), &(oldgtp.magneticField()) );
|
170 |
return gta;
|
171 |
}
|
172 |
|
173 |
void
|
174 |
ClosestApproachInRPhi::circleParameters(const TrackCharge& charge,
|
175 |
const GlobalVector& momentum,
|
176 |
const GlobalPoint& position,
|
177 |
double& xc, double& yc, double& r,
|
178 |
double bz)
|
179 |
{
|
180 |
|
181 |
// compute radius of circle
|
182 |
/** temporary code, to be replaced by call to curvature() when bug
|
183 |
* is fixed.
|
184 |
*/
|
185 |
// double bz = MagneticField::inInverseGeV(position).z();
|
186 |
|
187 |
// signed_r directed towards circle center, along F_Lorentz = q*v X B
|
188 |
double qob = charge/bz;
|
189 |
double signed_r = qob*momentum.transverse();
|
190 |
r = abs(signed_r);
|
191 |
/** end of temporary code
|
192 |
*/
|
193 |
|
194 |
// compute centre of circle
|
195 |
// double phi = momentum.phi();
|
196 |
// xc = signed_r*sin(phi) + position.x();
|
197 |
// yc = -signed_r*cos(phi) + position.y();
|
198 |
xc = position.x() + qob * momentum.y();
|
199 |
yc = position.y() - qob * momentum.x();
|
200 |
|
201 |
}
|
202 |
|
203 |
|
204 |
int
|
205 |
ClosestApproachInRPhi::transverseCoord(double cxa, double cya, double ra,
|
206 |
double cxb, double cyb, double rb,
|
207 |
double & xg1, double & yg1,
|
208 |
double & xg2, double & yg2)
|
209 |
{
|
210 |
int flag = 0;
|
211 |
double x1, y1, x2, y2;
|
212 |
|
213 |
// new reference frame with origin in (cxa, cya) and x-axis
|
214 |
// directed from (cxa, cya) to (cxb, cyb)
|
215 |
|
216 |
double d_ab = sqrt((cxb - cxa)*(cxb - cxa) + (cyb - cya)*(cyb - cya));
|
217 |
if (d_ab == 0) { // concentric circles
|
218 |
return 0;
|
219 |
}
|
220 |
// elements of rotation matrix
|
221 |
double u = (cxb - cxa) / d_ab;
|
222 |
double v = (cyb - cya) / d_ab;
|
223 |
|
224 |
// conditions for circle intersection
|
225 |
if (d_ab <= ra + rb && d_ab >= abs(rb - ra)) {
|
226 |
|
227 |
// circles cross each other
|
228 |
flag = 1;
|
229 |
|
230 |
// triangle (ra, rb, d_ab)
|
231 |
double cosphi = (ra*ra - rb*rb + d_ab*d_ab) / (2*ra*d_ab);
|
232 |
double sinphi2 = 1. - cosphi*cosphi;
|
233 |
if (sinphi2 < 0.) { sinphi2 = 0.; cosphi = 1.; }
|
234 |
|
235 |
// intersection points in new frame
|
236 |
double sinphi = sqrt(sinphi2);
|
237 |
x1 = ra*cosphi; y1 = ra*sinphi; x2 = x1; y2 = -y1;
|
238 |
}
|
239 |
else if (d_ab > ra + rb) {
|
240 |
|
241 |
// circles are external to each other
|
242 |
flag = 2;
|
243 |
|
244 |
// points of closest approach in new frame
|
245 |
// are on line between 2 centers
|
246 |
x1 = ra; y1 = 0; x2 = d_ab - rb; y2 = 0;
|
247 |
}
|
248 |
else if (d_ab < abs(rb - ra)) {
|
249 |
|
250 |
// circles are inside each other
|
251 |
flag = 2;
|
252 |
|
253 |
// points of closest approach in new frame are on line between 2 centers
|
254 |
// choose 2 closest points
|
255 |
double sign = 1.;
|
256 |
if (ra <= rb) sign = -1.;
|
257 |
x1 = sign*ra; y1 = 0; x2 = d_ab + sign*rb; y2 = 0;
|
258 |
}
|
259 |
else {
|
260 |
return 0;
|
261 |
}
|
262 |
|
263 |
// intersection points in global frame, transverse plane
|
264 |
xg1 = u*x1 - v*y1 + cxa; yg1 = v*x1 + u*y1 + cya;
|
265 |
xg2 = u*x2 - v*y2 + cxa; yg2 = v*x2 + u*y2 + cya;
|
266 |
|
267 |
return flag;
|
268 |
}
|
269 |
|
270 |
|
271 |
double
|
272 |
ClosestApproachInRPhi::zCoord(const GlobalVector& mom,
|
273 |
const GlobalPoint& pos,
|
274 |
double r, double xc, double yc,
|
275 |
double xg, double yg)
|
276 |
{
|
277 |
|
278 |
// starting point
|
279 |
double x = pos.x(); double y = pos.y(); double z = pos.z();
|
280 |
|
281 |
double px = mom.x(); double py = mom.y(); double pz = mom.z();
|
282 |
|
283 |
// rotation angle phi from starting point to crossing point (absolute value)
|
284 |
// -- compute sin(phi/2) if phi smaller than pi/4,
|
285 |
// -- cos(phi) if phi larger than pi/4
|
286 |
double phi = 0.;
|
287 |
double sinHalfPhi = sqrt((x-xg)*(x-xg) + (y-yg)*(y-yg))/(2*r);
|
288 |
if (sinHalfPhi < 0.383) { // sin(pi/8)
|
289 |
phi = 2*asin(sinHalfPhi);
|
290 |
}
|
291 |
else {
|
292 |
double cosPhi = ((x-xc)*(xg-xc) + (y-yc)*(yg-yc))/(r*r);
|
293 |
if (std::abs(cosPhi) > 1) cosPhi = (cosPhi > 0 ? 1 : -1);
|
294 |
phi = abs(acos(cosPhi));
|
295 |
}
|
296 |
// -- sign of phi
|
297 |
double signPhi = ((x - xc)*(yg - yc) - (xg - xc)*(y - yc) > 0) ? 1. : -1.;
|
298 |
|
299 |
// sign of track angular momentum
|
300 |
// if rotation is along angular momentum, delta z is along pz
|
301 |
double signOmega = ((x - xc)*py - (y - yc)*px > 0) ? 1. : -1.;
|
302 |
|
303 |
// delta z
|
304 |
// -- |dz| = |cos(theta) * path along helix|
|
305 |
// = |cos(theta) * arc length along circle / sin(theta)|
|
306 |
double dz = signPhi*signOmega*(pz/mom.transverse())*phi*r;
|
307 |
|
308 |
return z + dz;
|
309 |
}
|