ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/IPHCalignment2/TrackingTools/PatternTools/src/ClosestApproachInRPhi.cc
Revision: 1.1
Committed: Fri Nov 25 16:38:25 2011 UTC (13 years, 5 months ago) by econte
Content type: text/plain
Branch: MAIN
CVS Tags: TBD2011, TBD_2011, HEAD
Error occurred while calculating annotation data.
Log Message:
new IPHC alignment

File Contents

# Content
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 }