1 package org.lcsim.recon.tracking.trfcylplane;
2
3 import org.lcsim.recon.tracking.trfbase.Propagator;
4 import org.lcsim.recon.tracking.trfbase.PropDirected;
5 import org.lcsim.recon.tracking.trfbase.PropDir;
6 import org.lcsim.recon.tracking.trfbase.TrackVector;
7 import org.lcsim.recon.tracking.trfbase.TrackDerivative;
8 import org.lcsim.recon.tracking.trfzp.PropZZ;
9 import org.lcsim.recon.tracking.trfbase.PropStat;
10 import org.lcsim.recon.tracking.trfutil.TRFMath;
11 import org.lcsim.recon.tracking.trfutil.Assert;
12 import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
13 import org.lcsim.recon.tracking.trfzp.SurfZPlane;
14 import org.lcsim.recon.tracking.trfbase.Surface;
15 import org.lcsim.recon.tracking.trfbase.VTrack;
16 import org.lcsim.recon.tracking.trfbase.ETrack;
17
18 /**
19 * Propagates tracks from an XYPlane to a ZPlane in a constant magnetic field.
20 *<p>
21 * Propagation will fail if either the origin is not an XYPlane
22 * or destination is not a ZPlane.
23 * Propagator works incorrectly for tracks with very small curvatures.
24 *<p>
25 *@author Norman A. Graf
26 *@version 1.0
27 *
28 */
29 public class PropXYZ extends PropDirected
30 {
31
32 // attributes
33
34 private double _bfac;
35 private PropZZ _propzz;
36
37 // Assign track parameter indices.
38
39 private static final int IV = SurfXYPlane.IV;
40 private static final int IZ = SurfXYPlane.IZ;
41 private static final int IDVDU = SurfXYPlane.IDVDU;
42 private static final int IDZDU = SurfXYPlane.IDZDU;
43 private static final int IQP_XY = SurfXYPlane.IQP;
44
45 private static final int IX = SurfZPlane.IX;
46 private static final int IY = SurfZPlane.IY;
47 private static final int IDXDZ = SurfZPlane.IDXDZ;
48 private static final int IDYDZ = SurfZPlane.IDYDZ;
49 private static final int IQP_Z = SurfZPlane.IQP;
50
51 // static methods
52
53
54 /**
55 *Return a String representation of the class' type name.
56 *Included for completeness with the C++ version.
57 *
58 * @return A String representation of the class' type name.
59 */
60 public static String typeName()
61 { return "PropXYZ";
62 }
63
64 /**
65 *Return a String representation of the class' type name.
66 *Included for completeness with the C++ version.
67 *
68 * @return A String representation of the class' type name.
69 */
70 public static String staticType()
71 { return typeName();
72 }
73
74 /**
75 *Construct an instance from a constant solenoidal magnetic field in Tesla.
76 *
77 * @param bfield The magnetic field strength in Tesla.
78 */
79 public PropXYZ(double bfield)
80 {
81 _bfac = TRFMath.BFAC*bfield;
82 _propzz = new PropZZ(bfield);
83 }
84
85
86 /**
87 *Clone an instance.
88 *
89 * @return A Clone of this instance.
90 */
91 public Propagator newPropagator( )
92 {
93 return new PropXYZ( bField() );
94 }
95
96 /**
97 *Propagate a track without error in the specified direction
98 *and return the derivative matrix in deriv.
99 *
100 * The track parameters for a cylinder are:
101 * phi z alpha tan(lambda) curvature
102 *
103 * @param trv The VTrack to propagate.
104 * @param srf The Surface to which to propagate.
105 * @param dir The direction in which to propagate.
106 * @param deriv The track derivatives to update at the surface srf.
107 * @return The propagation status.
108 */
109 public PropStat vecDirProp( VTrack trv, Surface srf,
110 PropDir dir,TrackDerivative deriv)
111 {
112 PropStat pstat = new PropStat();
113 Propagator.reduceDirection(dir);
114
115 // Check destination is a ZPlane.
116 Assert.assertTrue( srf.pureType().equals(SurfZPlane.staticType()) );
117 if ( !srf.pureType( ).equals(SurfZPlane.staticType()) )
118 return pstat;
119
120 // Fetch phi of the destination plane
121
122 VTrack trv0 = trv; // we wish to modify this track!
123 TrackDerivative deriv1 = new TrackDerivative();
124 TrackDerivative deriv2 = new TrackDerivative();
125
126 if ( deriv != null ) pstat = vecTransformXYZ( _bfac, trv0, dir, deriv1 );
127 else pstat = vecTransformXYZ( _bfac, trv0, dir, deriv);
128
129 if ( ! pstat.success() ) return pstat;
130
131 if ( deriv != null ) pstat = _propzz.vecDirProp(trv0,srf,dir, deriv2);
132 else pstat = _propzz.vecDirProp(trv0,srf,dir, deriv);
133
134 if ( pstat.success() )
135 {
136 //trv = trv0;
137 if ( deriv != null ) deriv.set(deriv2.times(deriv1));
138 }
139
140 return pstat;
141 }
142
143 /**
144 *Propagate a track without error in the specified direction.
145 *
146 * The track parameters for a cylinder are:
147 * phi z alpha tan(lambda) curvature
148 *
149 * @param trv The VTrack to propagate.
150 * @param srf The Surface to which to propagate.
151 * @param dir The direction in which to propagate.
152 * @return The propagation status.
153 */
154 public PropStat vecDirProp( VTrack trv, Surface srf,
155 PropDir dir)
156 {
157 TrackDerivative deriv = null;
158 return vecDirProp(trv, srf, dir, deriv);
159
160 }
161
162 /**
163 *Return the strength of the magnetic field in Tesla.
164 *
165 * @return The strength of the magnetic field in Tesla.
166 */
167 public double bField()
168 {
169 return _bfac/TRFMath.BFAC;
170 }
171
172
173 /**
174 *Return a String representation of the class' type name.
175 *Included for completeness with the C++ version.
176 *
177 * @return A String representation of the class' type name.
178 */
179 public String type()
180 { return staticType();
181 }
182
183 /**
184 *output stream
185 * @return A String representation of this instance.
186 */
187 public String toString()
188 {
189 return "XYPlane-ZPlane propagation with constant "
190 + bField() + " Tesla field";
191
192 }
193
194
195
196 //**********************************************************************
197
198 // Private function to propagate a track without error
199 // The corresponding track parameters are:
200 // On ZPlane:
201 // z (cm) is fixed
202 // 0 - x (cm)
203 // 1 - y (cm)
204 // 2 - dx/dz
205 // 3 - dy/dz
206 // 4 - q/p p is momentum of a track, q is its charge
207 // On XYPlane:
208 // u (cm) is fixed
209 // 0 - v (cm)
210 // 1 - z (cm)
211 // 2 - dv/du
212 // 3 - dz/du
213 // 4 - q/p p is momentum of a track, q is its charge
214 // If deriv is nonzero, return the derivative matrix there.
215
216 private PropStat
217 vecTransformXYZ( double B, VTrack trv,
218 PropDir dir,
219 TrackDerivative deriv)
220 {
221
222 // construct return status
223 PropStat pstat = new PropStat();
224
225 // fetch the originating surface and vector
226 Surface srf1 = trv.surface();
227 // TrackVector vec1 = trv.get_vector();
228
229 // Check origin is a XYPlane.
230 Assert.assertTrue( srf1.pureType().equals(SurfXYPlane.staticType()) );
231 if ( !srf1.pureType( ).equals(SurfXYPlane.staticType()) )
232 return pstat;
233 SurfXYPlane sxyp1 = (SurfXYPlane) srf1;
234
235 // Fetch the dist and phi of the plane and the starting track vector.
236 int iphi = SurfXYPlane.NORMPHI;
237 int idist = SurfXYPlane.DISTNORM;
238 double sphi = sxyp1.parameter(iphi);
239 double u = sxyp1.parameter(idist);
240 TrackVector vec = trv.vector();
241 double b1 = vec.get(IV); // v
242 double b2 = vec.get(IZ); // z
243 double b3 = vec.get(IDVDU); // dv/du
244 double b4 = vec.get(IDZDU); // dz/du
245 double b5 = vec.get(IQP_XY); // q/p
246
247 // check that dz != 0
248 if(b4 == 0.) return pstat;
249
250 double zpos = b2;
251 double cos_sphi = Math.cos(sphi);
252 double sin_sphi = Math.sin(sphi);
253
254 double a1 = u*cos_sphi - b1*sin_sphi;
255 double a2 = b1*cos_sphi + u*sin_sphi;
256 double a3 = cos_sphi/b4 - b3/b4*sin_sphi;
257 double a4 = b3/b4*cos_sphi + sin_sphi/b4;
258 double a5 = b5;
259
260
261 // if dz/du*du/dt > 0 and forward : dz>0; dz/du*du/dt < 0 dz<0
262 // if dz/du*du/dt < 0 and backward : dz>0; dz/du*du/dt > 0 dz<0
263
264 int sign_du = 0;
265 if(trv.isForward()) sign_du = 1;
266 if(trv.isBackward()) sign_du = -1;
267 if(sign_du == 0)
268 {
269 System.out.println("PropXYZ._vec_propagate: Unknown direction of a track ");
270 System.exit(1);
271 }
272
273 int sign_dz = 0;
274 if(b4*sign_du > 0) sign_dz = 1;
275 if(b4*sign_du < 0) sign_dz = -1;
276
277 vec.set(IX , a1);
278 vec.set(IY , a2);
279 vec.set(IDXDZ , a3);
280 vec.set(IDYDZ , a4);
281 vec.set(IQP_Z , a5);
282
283 // Update trv
284 SurfZPlane szp = new SurfZPlane(zpos);
285 trv.setSurface(szp.newPureSurface());
286 trv.setVector(vec);
287
288 // set new direction of the track
289 if(sign_dz == 1) trv.setForward();
290 if(sign_dz == -1) trv.setBackward();
291
292 // Set the return status.
293 pstat.setSame();
294
295 // exit now if user did not ask for error matrix.
296 if ( deriv == null ) return pstat;
297
298 double a34_hat = sign_dz*Math.sqrt(a3*a3 + a4*a4);
299 double a34_hat2 = a34_hat*a34_hat;
300
301 double Rcos_phi = a3/Math.sqrt(1.+a34_hat2)*sign_dz;
302 double Rsin_phi = a4/Math.sqrt(1.+a34_hat2)*sign_dz;
303
304 // ddphi / db
305
306 double ddphi_db2 = -B*a5*Math.sqrt(a34_hat2+1.)*sign_dz;
307 double ddphi_db2_nob = -Math.sqrt(a34_hat2+1.)*sign_dz;
308
309 // da3 / db
310
311 double da3_db3 = -sin_sphi/b4;
312 double da3_db4 = (-cos_sphi + b3*sin_sphi)/(b4*b4);
313
314 // da4 / db
315
316 double da4_db3 = cos_sphi/b4;
317 double da4_db4 = (- sin_sphi - b3*cos_sphi)/(b4*b4);
318
319 // da1/ db
320
321 double da1n_db1 = -sin_sphi;
322 double da1n_db2 = Rcos_phi*ddphi_db2_nob;
323 double da1n_db3 = 0.;
324 double da1n_db4 = 0.;
325 double da1n_db5 = 0.;
326
327 // da2/ db
328
329 double da2n_db1 = cos_sphi;
330 double da2n_db2 = Rsin_phi*ddphi_db2_nob;
331 double da2n_db3 = 0.;
332 double da2n_db4 = 0.;
333 double da2n_db5 = 0.;
334
335 // da3/ db
336
337 double da3n_db1 = 0.;
338 double da3n_db2 = -a4*ddphi_db2;
339 double da3n_db3 = da3_db3;
340 double da3n_db4 = da3_db4;
341 double da3n_db5 = 0.;
342
343 // da4/ db
344
345 double da4n_db1 = 0.;
346 double da4n_db2 = a3*ddphi_db2;
347 double da4n_db3 = da4_db3;
348 double da4n_db4 = da4_db4;
349 double da4n_db5 = 0.;
350
351 // da5n
352
353 double da5n_db1 = 0.;
354 double da5n_db2 = 0.;
355 double da5n_db3 = 0.;
356 double da5n_db4 = 0.;
357 double da5n_db5 = 1.;
358
359 deriv.set(IX,IV , da1n_db1);
360 deriv.set(IX,IZ , da1n_db2);
361 deriv.set(IX,IDVDU , da1n_db3);
362 deriv.set(IX,IDZDU , da1n_db4);
363 deriv.set(IX,IQP_XY , da1n_db5);
364 deriv.set(IY,IV , da2n_db1);
365 deriv.set(IY,IZ , da2n_db2);
366 deriv.set(IY,IDVDU , da2n_db3);
367 deriv.set(IY,IDZDU , da2n_db4);
368 deriv.set(IY,IQP_XY , da2n_db5);
369 deriv.set(IDXDZ,IV , da3n_db1);
370 deriv.set(IDXDZ,IZ , da3n_db2);
371 deriv.set(IDXDZ,IDVDU , da3n_db3);
372 deriv.set(IDXDZ,IDZDU , da3n_db4);
373 deriv.set(IDXDZ,IQP_XY , da3n_db5);
374 deriv.set(IDYDZ,IV , da4n_db1);
375 deriv.set(IDYDZ,IZ , da4n_db2);
376 deriv.set(IDYDZ,IDVDU , da4n_db3);
377 deriv.set(IDYDZ,IDZDU , da4n_db4);
378 deriv.set(IDYDZ,IQP_XY , da4n_db5);
379 deriv.set(IQP_Z,IV , da5n_db1);
380 deriv.set(IQP_Z,IZ , da5n_db2);
381 deriv.set(IQP_Z,IDVDU , da5n_db3);
382 deriv.set(IQP_Z,IDZDU , da5n_db4);
383 deriv.set(IQP_Z,IQP_XY , da5n_db5);
384
385 return pstat;
386 }
387
388 }
389