1 package org.lcsim.util.swim;
2
3 import hep.physics.vec.Hep3Vector;
4
5 import org.lcsim.event.LCIOParameters;
6 import org.lcsim.event.Track;
7 import org.lcsim.spacegeom.CartesianPoint;
8 import org.lcsim.spacegeom.CartesianVector;
9 import org.lcsim.spacegeom.SpacePoint;
10 import org.lcsim.spacegeom.SpaceVector;
11
12 import static java.lang.Math.atan;
13 import static java.lang.Math.sqrt;
14 import static org.lcsim.constants.Constants.fieldConversion;
15 import static org.lcsim.event.LCIOParameters.ParameterName.phi0;
16 import static org.lcsim.event.LCIOParameters.ParameterName.omega;
17 import static org.lcsim.event.LCIOParameters.ParameterName.tanLambda;
18
19
20
21
22
23
24
25
26
27
28
29 public class HelixSwimmer {
30 protected double field;
31 protected Trajectory _trajectory;
32 protected SpaceVector _momentum;
33 protected double _charge;
34
35
36
37
38
39
40
41
42 public HelixSwimmer(double B) {
43 field = B * fieldConversion;
44 }
45
46
47
48
49
50
51
52
53
54
55
56 public void setTrack(Hep3Vector p0, SpacePoint r0, int iq) {
57 SpaceVector p = new CartesianVector(p0.v());
58 double phi = Math.atan2(p.y(), p.x());
59 double lambda = Math.atan2(p.z(), p.rxy());
60
61 if (iq != 0 && field != 0) {
62 double radius = p.rxy() / (iq * field);
63 _trajectory = new Helix(r0, radius, phi, lambda);
64 } else {
65 _trajectory = new Line(r0, phi, lambda);
66 }
67 _momentum = p;
68 _charge = iq;
69 }
70
71
72
73
74
75
76
77
78
79
80
81 public void setTrack(Hep3Vector p0, SpacePoint r0, double q) {
82 SpaceVector p = new CartesianVector(p0.v());
83 double phi = Math.atan2(p.y(), p.x());
84 double lambda = Math.atan2(p.z(), p.rxy());
85
86 if (q != 0 && field != 0) {
87 double radius = p.rxy() / (q * field);
88 _trajectory = new Helix(r0, radius, phi, lambda);
89 } else {
90 _trajectory = new Line(r0, phi, lambda);
91 }
92 _momentum = p;
93 _charge = q;
94 }
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109 @Deprecated
110 public void setTrack(Hep3Vector p, Hep3Vector r0, int iq) {
111 double pt = sqrt(p.x() * p.x() + p.y() * p.y());
112 double phi = Math.atan2(p.y(), p.x());
113 double lambda = Math.atan2(p.z(), pt);
114
115 if (iq != 0 && field != 0) {
116 double radius = pt / (iq * field);
117 _trajectory = new Helix(r0, radius, phi, lambda);
118 } else {
119 _trajectory = new Line(r0, phi, lambda);
120 }
121 _momentum = new CartesianVector(p.v());
122 _charge = iq;
123 }
124
125
126
127
128
129
130
131
132
133 public void setTrack(Track t) {
134 double pt = sqrt(t.getPX() * t.getPX() + t.getPY() * t.getPY());
135 LCIOParameters parameters = new LCIOParameters(t.getTrackParameters(), pt);
136
137 SpacePoint ref = new CartesianPoint(t.getReferencePoint());
138 SpacePoint origin = LCIOParameters.Parameters2Position(parameters, ref);
139 _trajectory = new Helix(origin, 1 / parameters.get(omega), parameters
140 .get(phi0), atan(parameters.get(tanLambda)));
141
142 _momentum = new CartesianVector(LCIOParameters.Parameters2Momentum(parameters).v());
143 }
144
145
146
147
148
149
150
151
152 @Deprecated
153 public SpacePoint getPointAtDistance(double alpha) {
154 return getPointAtLength(alpha);
155 }
156
157
158
159
160
161
162
163
164 public SpacePoint getPointAtLength(double alpha) {
165 if (_trajectory == null) {
166 throw new RuntimeException("Trajectory not set");
167 }
168 return _trajectory.getPointAtDistance(alpha);
169 }
170
171 public double getDistanceToRadius(double r) {
172 if (_trajectory == null) {
173 throw new RuntimeException("Trajectory not set");
174 }
175 return _trajectory.getDistanceToInfiniteCylinder(r);
176 }
177 public double getDistanceToPolyhedra(double r, int nsides)
178 {
179 if (_trajectory == null) {
180 throw new RuntimeException("Trajectory not set");
181 }
182 double s = _trajectory.getDistanceToInfiniteCylinder(r);
183 if(Double.isNaN(s))return s;
184 return ((Helix) (_trajectory)).getDistanceToPolyhedra(r,nsides);
185 }
186
187 public double getDistanceToZ(double z) {
188 if (_trajectory == null) {
189 throw new RuntimeException("Trajectory not set");
190 }
191 double result = _trajectory.getDistanceToZPlane(z);
192 if (result < 0) {
193 result = _trajectory.getDistanceToZPlane(-z);
194 }
195 return result;
196 }
197
198 public double getDistanceToCylinder(double r, double z) {
199 double x1 = getDistanceToRadius(r);
200 double x2 = getDistanceToZ(z);
201 return Double.isNaN(x1) ? x2 : Math.min(x1, x2);
202 }
203
204
205
206
207
208
209
210
211
212 public double getTrackLengthToPoint(Hep3Vector point) {
213 return _trajectory.getDistanceToPoint(point);
214 }
215
216 public double getDistanceToPoint(Hep3Vector point) {
217 return getTrackLengthToPoint(point);
218 }
219
220
221
222
223
224
225
226
227 public SpaceVector getMomentumAtLength(double alpha) {
228
229 SpaceVector unitDirection = _trajectory.getUnitTangentAtLength(alpha);
230 double magnitude = _momentum.rxy();
231
232 return VectorArithmetic.multiply(unitDirection, magnitude);
233 }
234
235 public Trajectory getTrajectory() {
236 return _trajectory;
237 }
238 }