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.trfbase.PropStat;
9 import org.lcsim.recon.tracking.trfutil.TRFMath;
10 import org.lcsim.recon.tracking.trfutil.Assert;
11 import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
12 import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
13 import org.lcsim.recon.tracking.trfxyp.PropXYXY;
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 a Cylinder to a XYPlane in a constant field.
20 *<p>
21 * Propagation will fail if either the origin is not a Cylinder
22 * or destination is not a XYPlane.
23 * Propagator works incorrectly for tracks with very small curvatures.
24 *<p>
25 *
26 *@author Norman A. Graf
27 *@version 1.0
28 *
29 */
30
31 public class PropCylXY extends PropDirected
32 {
33
34 // Assign track parameter indices.
35
36 private static final int IV = SurfXYPlane.IV;
37 private static final int IZC = 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 IPHI = SurfCylinder.IPHI;
43 private static final int IZ = SurfCylinder.IZ;
44 private static final int IALF = SurfCylinder.IALF;
45 private static final int ITLM = SurfCylinder.ITLM;
46 private static final int IQPT = SurfCylinder.IQPT;
47
48
49 // attributes
50
51 private double _bfac;
52 private PropXYXY _propxyxy;
53
54 // static methods
55
56 /**
57 *Return a String representation of the class' type name.
58 *Included for completeness with the C++ version.
59 *
60 * @return A String representation of the class' type name.
61 */
62 public static String typeName()
63 { return "PropCylXY"; }
64
65 /**
66 *Return a String representation of the class' type name.
67 *Included for completeness with the C++ version.
68 *
69 * @return A String representation of the class' type name.
70 */
71 public static String staticType()
72 { return typeName(); }
73
74
75
76 /**
77 *Construct an instance from a constant solenoidal magnetic field in Tesla.
78 *
79 * @param bfield The magnetic field strength in Tesla.
80 */
81 public PropCylXY(double bfield)
82 {
83 _bfac = TRFMath.BFAC*bfield;
84 _propxyxy = new PropXYXY(bfield);
85
86 }
87
88 /**
89 *Clone an instance.
90 *
91 * @return A Clone of this instance.
92 */
93 public Propagator newPropagator( )
94 {
95 return new PropCylXY( bField() );
96
97 }
98
99 /**
100 *Propagate a track without error in the specified direction.
101 *
102 * The track parameters for a cylinder are:
103 * phi z alpha tan(lambda) curvature
104 *
105 * @param trv The VTrack to propagate.
106 * @param srf The Surface to which to propagate.
107 * @param dir The direction in which to propagate.
108 * @return The propagation status.
109 */
110 public PropStat vecDirProp( VTrack trv, Surface srf,
111 PropDir dir)
112 {
113 TrackDerivative deriv = null;
114 return vecDirProp(trv, srf, dir, deriv);
115
116 }
117
118 /**
119 *Propagate a track without error in the specified direction
120 *and return the derivative matrix in deriv.
121 *
122 * The track parameters for a cylinder are:
123 * phi z alpha tan(lambda) curvature
124 *
125 * @param trv The VTrack to propagate.
126 * @param srf The Surface to which to propagate.
127 * @param dir The direction in which to propagate.
128 * @param deriv The track derivatives to update at the surface srf.
129 * @return The propagation status.
130 */
131 public PropStat vecDirProp( VTrack trv, Surface srf,
132 PropDir dir, TrackDerivative deriv )
133 {
134 PropStat pstat = new PropStat();
135 Propagator.reduceDirection(dir);
136
137 // Check destination is a XYPlane.
138 Assert.assertTrue( srf.pureType().equals(SurfXYPlane.staticType()) );
139 if ( !srf.pureType( ).equals(SurfXYPlane.staticType()) )
140 return pstat;
141 SurfXYPlane sxyp2 = ( SurfXYPlane ) srf;
142
143 // Fetch phi of the destination plane
144
145 int iphi = SurfXYPlane.NORMPHI;
146 double phi_n = sxyp2.parameter(iphi);
147 VTrack trv0 = trv; // want to change this vector.
148 TrackDerivative deriv1 = new TrackDerivative();
149 TrackDerivative deriv2 = new TrackDerivative();
150
151 if ( deriv != null ) pstat = vecTransformCylXY( _bfac, trv0, phi_n, dir,deriv1 );
152 else pstat = vecTransformCylXY( _bfac, trv0, phi_n, dir, deriv );
153
154 if ( ! pstat.success() ) return pstat;
155
156 if ( deriv != null ) pstat = _propxyxy.vecDirProp(trv0,srf,dir,deriv2);
157 else pstat = _propxyxy.vecDirProp(trv0,srf,dir, deriv);
158
159 if ( pstat.success() )
160 {
161 trv = trv0;
162 if ( deriv != null ) deriv.set(deriv2.times(deriv1));
163 }
164
165 return pstat;
166 }
167
168 /**
169 *Return the strength of the magnetic field in Tesla.
170 *
171 * @return The strength of the magnetic field in Tesla.
172 */
173 public double bField()
174 {
175 return _bfac/TRFMath.BFAC;
176 }
177
178 /**
179 *Return a String representation of the class' type name.
180 *Included for completeness with the C++ version.
181 *
182 * @return A String representation of the class' type name.
183 */
184 public String type()
185 { return staticType(); }
186
187 /**
188 *output stream
189 * @return A String representation of this instance.
190 */
191 public String toString()
192 {
193 return "Cylinder-XYPlane propagation with constant "
194 + bField() + " Tesla field";
195
196 }
197
198 //**********************************************************************
199
200 // Private function to propagate a track without error
201 // The corresponding track parameters are:
202 // On Cylinder:
203 // r (cm) is fixed
204 // 0 - phi
205 // 1 - z (cm)
206 // 2 - alpha = phi_dir - phi; tan(alpha) = r*dphi/dr
207 // 3 - sin(lambda) = dz/ds; tan(lambda) = dz/dsT
208 // 4 - q/pT (1/GeV/c) (pT is component of p parallel to cylinder)
209 // On XYPlane:
210 // u (cm) is fixed
211 // 0 - v (cm)
212 // 1 - z (cm)
213 // 2 - dv/du
214 // 3 - dz/du
215 // 4 - q/p p is momentum of a track, q is its charge
216 // If pderiv is nonzero, return the derivative matrix there.
217
218 PropStat
219 vecTransformCylXY( double B, VTrack trv, double phi_n,
220 PropDir dir,
221 TrackDerivative deriv )
222 {
223
224 // construct return status
225 PropStat pstat = new PropStat();
226
227 // fetch the originating surface and vector
228 Surface srf1 = trv.surface();
229 // TrackVector vec1 = trv.get_vector();
230
231 // Check origin is a Cylinder.
232 Assert.assertTrue( srf1.pureType().equals(SurfCylinder.staticType()) );
233 if ( !srf1.pureType( ).equals(SurfCylinder.staticType()) )
234 return pstat;
235 SurfCylinder scy1 = ( SurfCylinder ) srf1;
236
237 // Fetch the R of the cylinder and the starting track vector.
238 int ir = SurfCylinder.RADIUS;
239 double Rcyl = scy1.parameter(ir);
240
241 TrackVector vec = trv.vector();
242 double c1 = vec.get(IPHI); // phi
243 double c2 = vec.get(IZ); // z
244 double c3 = vec.get(IALF); // alpha
245 double c4 = vec.get(ITLM); // tan(lambda)
246 double c5 = vec.get(IQPT); // q/pt
247
248 // rotate coordinate system on phi_n
249
250 c1 -= phi_n;
251 if(c1 < 0.) c1 += TRFMath.TWOPI;
252
253 double cos_c1 = Math.cos(c1);
254
255 double u = Rcyl*cos_c1;
256 if( u < 0. )
257 {
258 u = - u;
259 cos_c1 = - cos_c1;
260 c1 -= Math.PI;
261 if(c1 < 0.) c1 += TRFMath.TWOPI;
262 phi_n = phi_n + Math.PI;
263 if( phi_n >= TRFMath.TWOPI) phi_n -= TRFMath.TWOPI;
264 }
265
266 double sin_c1 = Math.sin(c1);
267 double cos_dir = Math.cos(c1+c3);
268 double sin_dir = Math.sin(c1+c3);
269 double c4_hat2 = 1+c4*c4;
270 double c4_hat = Math.sqrt(c4_hat2);
271
272 // check if du == 0 ( that is track moves parallel to the destination plane )
273 // du = pt*cos_dir
274 if(cos_dir/c5 == 0.) return pstat;
275
276 double tan_dir = sin_dir/cos_dir;
277
278 double b1 = Rcyl*sin_c1;
279 double b2 = c2;
280 double b3 = tan_dir;
281 double b4 = c4/cos_dir;
282 double b5 = c5/c4_hat;
283
284 int sign_du = 0;
285 if(cos_dir > 0) sign_du = 1;
286 if(cos_dir < 0) sign_du = -1;
287
288 vec.set(IV , b1);
289 vec.set(IZ , b2);
290 vec.set(IDVDU , b3);
291 vec.set(IDZDU , b4);
292 vec.set(IQP_XY , b5);
293
294 // Update trv
295 SurfXYPlane sxyp = new SurfXYPlane(u,phi_n);
296 trv.setSurface(sxyp.newPureSurface());
297 trv.setVector(vec);
298
299 // set new direction of the track
300 if(sign_du == 1) trv.setForward();
301 if(sign_du == -1) trv.setBackward();
302
303 // Set the return status.
304 pstat.setSame();
305
306 // exit now if user did not ask for error matrix.
307 if ( deriv == null ) return pstat;
308
309 double b34_hat = Math.sqrt(1 + b3*b3 + b4*b4);
310 double invert_rsinphi = (b5*B)*sign_du*b34_hat;
311
312
313 // du_dc
314
315 double du_dc1 = -Rcyl*sin_c1;
316
317 // db1_dc
318
319 double db1_dc1 = Rcyl*cos_c1;
320
321 // db2_dc
322
323 double db2_dc2 = 1.;
324
325 // db3_dc
326
327 double db3_dc1 = 1./(cos_dir*cos_dir);
328 double db3_dc3 = 1./(cos_dir*cos_dir);
329
330 // db4_dc
331
332 double db4_dc1 = b4*tan_dir;
333 double db4_dc3 = b4*tan_dir;
334 double db4_dc4 = 1/cos_dir;
335
336 // db5_dc
337
338 double db5_dc4 = -c4*c5/(c4_hat*c4_hat2);
339 double db5_dc5 = 1./c4_hat;
340
341 // db3_n_db
342
343 double db3_n_du = -(1. + b3*b3)*invert_rsinphi;
344
345 // db4_n_db
346
347 double db4_n_du = -b4*b3*invert_rsinphi;
348
349
350 // db1_n_dc
351
352 double db1_n_dc1 = db1_dc1 - b3 * du_dc1;
353 double db1_n_dc2 = 0.;
354 double db1_n_dc3 = 0.;
355 double db1_n_dc4 = 0.;
356 double db1_n_dc5 = 0.;
357
358 // db2_n_dc
359
360 double db2_n_dc1 = -b4 * du_dc1;
361 double db2_n_dc2 = db2_dc2;
362 double db2_n_dc3 = 0.;
363 double db2_n_dc4 = 0.;
364 double db2_n_dc5 = 0.;
365
366 // db3_n_dc
367
368 double db3_n_dc1 = db3_dc1 + db3_n_du * du_dc1;
369 double db3_n_dc2 = 0.;
370 double db3_n_dc3 = db3_dc3;
371 double db3_n_dc4 = 0.;
372 double db3_n_dc5 = 0.;
373
374 // db4_n_dc
375
376 double db4_n_dc1 = db4_dc1 + db4_n_du * du_dc1;
377 double db4_n_dc2 = 0.;
378 double db4_n_dc3 = db4_dc3;
379 double db4_n_dc4 = db4_dc4;
380 double db4_n_dc5 = 0.;
381
382 // db5_n_dc
383
384 double db5_n_dc1 = 0.;
385 double db5_n_dc2 = 0.;
386 double db5_n_dc3 = 0.;
387 double db5_n_dc4 = db5_dc4;
388 double db5_n_dc5 = db5_dc5;
389
390 deriv.set(IV,IPHI , db1_n_dc1);
391 deriv.set(IV,IZ , db1_n_dc2);
392 deriv.set(IV,IALF , db1_n_dc3);
393 deriv.set(IV,ITLM , db1_n_dc4);
394 deriv.set(IV,IQPT , db1_n_dc5);
395 deriv.set(IZ,IPHI , db2_n_dc1);
396 deriv.set(IZ,IZ , db2_n_dc2);
397 deriv.set(IZ,IALF , db2_n_dc3);
398 deriv.set(IZ,ITLM , db2_n_dc4);
399 deriv.set(IZ,IQPT , db2_n_dc5);
400 deriv.set(IDVDU,IPHI , db3_n_dc1);
401 deriv.set(IDVDU,IZ , db3_n_dc2);
402 deriv.set(IDVDU,IALF , db3_n_dc3);
403 deriv.set(IDVDU,ITLM , db3_n_dc4);
404 deriv.set(IDVDU,IQPT , db3_n_dc5);
405 deriv.set(IDZDU,IPHI , db4_n_dc1);
406 deriv.set(IDZDU,IZ , db4_n_dc2);
407 deriv.set(IDZDU,IALF , db4_n_dc3);
408 deriv.set(IDZDU,ITLM , db4_n_dc4);
409 deriv.set(IDZDU,IQPT , db4_n_dc5);
410 deriv.set(IQP_XY,IPHI , db5_n_dc1);
411 deriv.set(IQP_XY,IZ , db5_n_dc2);
412 deriv.set(IQP_XY,IALF , db5_n_dc3);
413 deriv.set(IQP_XY,ITLM , db5_n_dc4);
414 deriv.set(IQP_XY,IQPT , db5_n_dc5);
415
416 return pstat;
417 }
418
419
420
421 }