1 package org.lcsim.recon.tracking.trfcylplane;
2
3
4 import org.lcsim.recon.tracking.trfbase.Propagator;
5 import org.lcsim.recon.tracking.trfbase.PropDirected;
6 import org.lcsim.recon.tracking.trfbase.PropDir;
7 import org.lcsim.recon.tracking.trfbase.TrackVector;
8 import org.lcsim.recon.tracking.trfbase.TrackDerivative;
9 import org.lcsim.recon.tracking.trfxyp.PropXYXY;
10 import org.lcsim.recon.tracking.trfbase.PropStat;
11 import org.lcsim.recon.tracking.trfutil.TRFMath;
12 import org.lcsim.recon.tracking.trfutil.Assert;
13 import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
14 import org.lcsim.recon.tracking.trfzp.SurfZPlane;
15 import org.lcsim.recon.tracking.trfbase.Surface;
16 import org.lcsim.recon.tracking.trfbase.VTrack;
17 import org.lcsim.recon.tracking.trfbase.ETrack;
18
19 /**
20 * Propagates tracks from a ZPlane to an XYPlane in a constant magnetic field.
21 *<p>
22 * Propagation will fail if either the origin is not a ZPlane
23 * or destination is not a XYPlane.
24 * Propagator works incorrectly for tracks with very small curvatures
25 *<p>
26 *@author Norman A. Graf
27 *@version 1.0
28 *
29 */
30
31 public class PropZXY extends PropDirected
32 {
33
34 // Assign track parameter indices.
35
36 private static final int IV = SurfXYPlane.IV;
37 private static final int IZ = SurfXYPlane.IZ;
38 private static final int IDVDU = SurfXYPlane.IDVDU;
39 private static final int IDZDU = SurfXYPlane.IDZDU;
40 private static final int IQP_XY = SurfXYPlane.IQP;
41
42 private static final int IX = SurfZPlane.IX;
43 private static final int IY = SurfZPlane.IY;
44 private static final int IDXDZ = SurfZPlane.IDXDZ;
45 private static final int IDYDZ = SurfZPlane.IDYDZ;
46 private static final int IQP_Z = SurfZPlane.IQP;
47
48
49 // attributes
50
51 private double _bfac;
52 private PropXYXY _propxyxy;
53
54 // static methods
55
56 // Return the type name.
57 /**
58 *Return a String representation of the class' type name.
59 *Included for completeness with the C++ version.
60 *
61 * @return A String representation of the class' type name.
62 */
63 public static String typeName()
64 { return "PropZXY";
65 }
66
67 // Return the type.
68 /**
69 *Return a String representation of the class' type name.
70 *Included for completeness with the C++ version.
71 *
72 * @return A String representation of the class' type name.
73 */
74 public static String staticType()
75 { return typeName();
76 }
77
78
79
80 // constructor from magnetic field in Tesla.
81 /**
82 *Construct an instance from a constant solenoidal magnetic field in Tesla.
83 *
84 * @param bfield The magnetic field strength in Tesla.
85 */
86 public PropZXY(double bfield)
87 {
88 _bfac = TRFMath.BFAC*bfield;
89 _propxyxy = new PropXYXY(bfield);
90 }
91
92
93 // Clone.
94 /**
95 *Clone an instance.
96 *
97 * @return A Clone of this instance.
98 */
99 public Propagator newPropagator( )
100 {
101 return new PropZXY(bField());
102
103 }
104
105 /**
106 *Propagate a track without error in the specified direction
107 *and return the derivative matrix in deriv.
108 *
109 * The track parameters for a cylinder are:
110 * phi z alpha tan(lambda) curvature
111 *
112 * @param trv The VTrack to propagate.
113 * @param srf The Surface to which to propagate.
114 * @param dir The direction in which to propagate.
115 * @param deriv The track derivatives to update at the surface srf.
116 * @return The propagation status.
117 */
118 public PropStat vecDirProp( VTrack trv, Surface srf,
119 PropDir dir, TrackDerivative deriv )
120 {
121
122 PropStat pstat = new PropStat();
123 Propagator.reduceDirection(dir);
124
125 // Check destination is a XYPlane.
126 Assert.assertTrue( srf.pureType().equals(SurfXYPlane.staticType()) );
127 if ( !srf.pureType( ).equals(SurfXYPlane.staticType()) )
128 return pstat;
129 SurfXYPlane sxyp2 = ( SurfXYPlane ) srf;
130
131 // Fetch phi of the destination plane
132
133 int iphi = SurfXYPlane.NORMPHI;
134 double phi_n = sxyp2.parameter(iphi);
135 VTrack trv0 = trv; // we want to change this vector!
136 TrackDerivative deriv1 = new TrackDerivative();
137 TrackDerivative deriv2 = new TrackDerivative();
138 TrackDerivative lderiv = deriv;
139 if ( deriv != null ) pstat = vecTransformZXY( _bfac, trv0, phi_n, dir,deriv1 );
140 else pstat = vecTransformZXY( _bfac, trv0, phi_n, dir, lderiv );
141 if ( ! pstat.success() ) return pstat;
142
143 if ( deriv != null ) pstat = _propxyxy.vecDirProp(trv0,srf,dir,deriv2);
144 else pstat = _propxyxy.vecDirProp(trv0,srf,dir, lderiv);
145 if ( pstat.success() )
146 {
147 trv = trv0;
148 if ( deriv != null ) deriv.set(deriv2.times(deriv1));
149 }
150 return pstat;
151 }
152
153 /**
154 *Propagate a track without error in the specified direction.
155 *
156 * The track parameters for a cylinder are:
157 * phi z alpha tan(lambda) curvature
158 *
159 * @param trv The VTrack to propagate.
160 * @param srf The Surface to which to propagate.
161 * @param dir The direction in which to propagate.
162 * @return The propagation status.
163 */
164 public PropStat vecDirProp( VTrack trv, Surface srf,
165 PropDir dir)
166 {
167 TrackDerivative deriv = null;
168 return vecDirProp(trv, srf, dir, deriv);
169
170 }
171
172 /**
173 *Return the strength of the magnetic field in Tesla.
174 *
175 * @return The strength of the magnetic field in Tesla.
176 */
177 public double bField()
178 {
179 return _bfac/TRFMath.BFAC;
180 }
181
182 // Return the type.
183 /**
184 *Return a String representation of the class' type name.
185 *Included for completeness with the C++ version.
186 *
187 * @return A String representation of the class' type name.
188 */
189 public String type()
190 { return staticType();
191 }
192
193 /**
194 *output stream
195 * @return A String representation of this instance.
196 */
197 public String toString()
198 {
199 return "ZPlane-XYPlane propagation with constant "
200 + bField() + " Tesla field";
201
202 }
203
204 //**********************************************************************
205
206 // Private function to propagate a track without error
207 // The corresponding track parameters are:
208 // On ZPlane:
209 // z (cm) is fixed
210 // 0 - x (cm)
211 // 1 - y (cm)
212 // 2 - dx/dz
213 // 3 - dy/dz
214 // 4 - q/p p is momentum of a track, q is its charge
215 // On XYPlane:
216 // u (cm) is fixed
217 // 0 - v (cm)
218 // 1 - z (cm)
219 // 2 - dv/du
220 // 3 - dz/du
221 // 4 - q/p p is momentum of a track, q is its charge
222 // If pderiv is nonzero, return the derivative matrix there.
223
224 private PropStat
225 vecTransformZXY( double B, VTrack trv, double phi_n,
226 PropDir dir,
227 TrackDerivative deriv )
228 {
229
230 // construct return status
231 PropStat pstat = new PropStat();
232
233 // fetch the originating surface and vector
234 Surface srf1 = trv.surface();
235 // TrackVector vec1 = trv.get_vector();
236
237 // Check origin is a ZPlane.
238 Assert.assertTrue( srf1.pureType().equals(SurfZPlane.staticType()) );
239 if ( !srf1.pureType( ).equals(SurfZPlane.staticType()) )
240 return pstat;
241 SurfZPlane szp1 = ( SurfZPlane ) srf1;
242
243
244 // Fetch the zpos's of the planes and the starting track vector.
245 int izpos = SurfZPlane.ZPOS;
246
247 double zpos_0 = szp1.parameter(izpos);
248
249 TrackVector vec = trv.vector();
250 double a1 = vec.get(IX); // x
251 double a2 = vec.get(IY); // y
252 double a3 = vec.get(IDXDZ); // dx/dz
253 double a4 = vec.get(IDYDZ); // dy/dz
254 double a5 = vec.get(IQP_Z); // q/p
255
256 int sign_dz = 0;
257 if(trv.isForward()) sign_dz = 1;
258 if(trv.isBackward()) sign_dz = -1;
259 if(sign_dz == 0)
260 {
261 System.out.println("PropZXY._vec_propagate: Unknown direction of a track ");
262 System.exit(1);
263 }
264
265 double sinphi_n = Math.sin(phi_n);
266 double cosphi_n = Math.cos(phi_n);
267
268 double u = a1*cosphi_n + a2*sinphi_n;
269 if( u < 0. )
270 {
271 u = - u;
272 sinphi_n = - sinphi_n;
273 cosphi_n = - cosphi_n;
274 phi_n = phi_n + Math.PI;
275 if( phi_n >= TRFMath.TWOPI) phi_n -= TRFMath.TWOPI;
276 }
277 // check if du == 0 ( that is track moves parallel to the destination plane )
278 double du_dz = a3*cosphi_n + a4*sinphi_n;
279 if(du_dz == 0.) return pstat;
280
281 double a_hat_phin = 1./du_dz;
282 double a_hat_phin2 = a_hat_phin*a_hat_phin;
283
284 // double u = a1*cosphi_n + a2*sinphi_n;
285 double b1 = -a1*sinphi_n + a2*cosphi_n;
286 double b2 = zpos_0;
287 double b3 = (a4*cosphi_n - a3*sinphi_n)*a_hat_phin;
288 double b4 = a_hat_phin;
289 double b5 = a5;
290
291 int sign_du = 0;
292 if(b4*sign_dz > 0) sign_du = 1;
293 if(b4*sign_dz < 0) sign_du = -1;
294
295 vec.set(IV , b1);
296 vec.set(IZ , b2);
297 vec.set(IDVDU , b3);
298 vec.set(IDZDU , b4);
299 vec.set(IQP_XY , b5);
300
301 // Update trv
302 SurfXYPlane sxyp = new SurfXYPlane(u,phi_n);
303 trv.setSurface(sxyp.newPureSurface());
304 trv.setVector(vec);
305
306 // set new direction of the track
307 if(sign_du == 1) trv.setForward();
308 if(sign_du == -1) trv.setBackward();
309
310 // Set the return status.
311 pstat.setSame();
312
313 // exit now if user did not ask for error matrix.
314 if ( deriv == null ) return pstat;
315
316 double b34_hat = Math.sqrt(1 + b3*b3 + b4*b4);
317 double rsinphi = (b5*B)*sign_du*b34_hat;
318
319 // du_da
320
321 double du_da1 = cosphi_n;
322 double du_da2 = sinphi_n;
323
324 // db1_da
325
326 double db1_da1 = - sinphi_n;
327 double db1_da2 = cosphi_n;
328
329 // db3_da
330
331 double db3_da3 = -a4*a_hat_phin2;
332 double db3_da4 = a3*a_hat_phin2;
333
334 // db4_da
335
336 double db4_da3 = - cosphi_n*a_hat_phin2;
337 double db4_da4 = - sinphi_n*a_hat_phin2;
338
339
340 // db3_n_db
341
342 double db3_n_du = -(1. + b3*b3)*rsinphi;
343
344 // db4_n_db
345
346 double db4_n_du = -b4*b3*rsinphi;
347
348 // db5_n_db
349
350
351 // db1_n_da
352
353 double db1_n_da1 = -b3 * du_da1 + db1_da1;
354 double db1_n_da2 = -b3 * du_da2 + db1_da2;
355 double db1_n_da3 = 0.;
356 double db1_n_da4 = 0.;
357 double db1_n_da5 = 0.;
358
359 // db2_n_da
360
361 double db2_n_da1 = -b4 * du_da1;
362 double db2_n_da2 = -b4 * du_da2;
363 double db2_n_da3 = 0.;
364 double db2_n_da4 = 0.;
365 double db2_n_da5 = 0.;
366
367 // db3_n_da
368
369 double db3_n_da1 = db3_n_du * du_da1;
370 double db3_n_da2 = db3_n_du * du_da2;
371 double db3_n_da3 = db3_da3;
372 double db3_n_da4 = db3_da4;
373 double db3_n_da5 = 0.;
374
375 // db4_n_da
376
377 double db4_n_da1 = db4_n_du * du_da1;
378 double db4_n_da2 = db4_n_du * du_da2;
379 double db4_n_da3 = db4_da3;
380 double db4_n_da4 = db4_da4;
381 double db4_n_da5 = 0.;
382
383 // db5_n_da
384
385 double db5_n_da1 = 0.;
386 double db5_n_da2 = 0.;
387 double db5_n_da3 = 0.;
388 double db5_n_da4 = 0.;
389 double db5_n_da5 = 1.;
390
391 deriv.set(IV,IX , db1_n_da1);
392 deriv.set(IV,IY , db1_n_da2);
393 deriv.set(IV,IDXDZ , db1_n_da3);
394 deriv.set(IV,IDYDZ , db1_n_da4);
395 deriv.set(IV,IQP_Z , db1_n_da5);
396 deriv.set(IZ,IX , db2_n_da1);
397 deriv.set(IZ,IY , db2_n_da2);
398 deriv.set(IZ,IDXDZ , db2_n_da3);
399 deriv.set(IZ,IDYDZ , db2_n_da4);
400 deriv.set(IZ,IQP_Z , db2_n_da5);
401 deriv.set(IDVDU,IX , db3_n_da1);
402 deriv.set(IDVDU,IY , db3_n_da2);
403 deriv.set(IDVDU,IDXDZ , db3_n_da3);
404 deriv.set(IDVDU,IDYDZ , db3_n_da4);
405 deriv.set(IDVDU,IQP_Z , db3_n_da5);
406 deriv.set(IDZDU,IX , db4_n_da1);
407 deriv.set(IDZDU,IY , db4_n_da2);
408 deriv.set(IDZDU,IDXDZ , db4_n_da3);
409 deriv.set(IDZDU,IDYDZ , db4_n_da4);
410 deriv.set(IDZDU,IQP_Z , db4_n_da5);
411 deriv.set(IQP_XY,IX , db5_n_da1);
412 deriv.set(IQP_XY,IY , db5_n_da2);
413 deriv.set(IQP_XY,IDXDZ , db5_n_da3);
414 deriv.set(IQP_XY,IDYDZ , db5_n_da4);
415 deriv.set(IQP_XY,IQP_Z , db5_n_da5);
416
417 return pstat;
418 }
419
420
421 }
422