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.trfzp.SurfZPlane;
13 import org.lcsim.recon.tracking.trfzp.PropZZ;
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 ZPlane in a constant magnetic field.
20 *<p>
21 * Propagation will fail if either the origin is not a Cylinder
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
30 public class PropCylZ extends PropDirected
31 {
32
33 // Assign track parameter indices.
34 private static final int IX = SurfZPlane.IX;
35 private static final int IY = SurfZPlane.IY;
36 private static final int IDXDZ = SurfZPlane.IDXDZ;
37 private static final int IDYDZ = SurfZPlane.IDYDZ;
38 private static final int IQP_Z = SurfZPlane.IQP;
39
40 private static final int IPHI = SurfCylinder.IPHI;
41 private static final int IZ_CYL= SurfCylinder.IZ;
42 private static final int IALF = SurfCylinder.IALF;
43 private static final int ITLM = SurfCylinder.ITLM;
44 private static final int IQPT = SurfCylinder.IQPT;
45
46
47 // attributes
48
49 private double _bfac;
50 private PropZZ _propzz;
51
52
53 // static methods
54
55 /**
56 *Return a String representation of the class' type name.
57 *Included for completeness with the C++ version.
58 *
59 * @return A String representation of the class' type name.
60 */
61 public static String typeName()
62 { return "PropCylZ";
63 }
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 *Construct an instance from a constant solenoidal magnetic field in Tesla.
77 *
78 * @param bfield The magnetic field strength in Tesla.
79 */
80 public PropCylZ(double bfield)
81 {
82 _bfac = TRFMath.BFAC*bfield;
83 _propzz = new PropZZ(bfield);
84
85 }
86
87 /**
88 *Clone an instance.
89 *
90 * @return A Clone of this instance.
91 */
92 public Propagator newPropagator( )
93 {
94 return new PropCylZ( bField() );
95
96 }
97
98 /**
99 *Propagate a track without error in the specified direction.
100 *
101 * The track parameters for a cylinder are:
102 * phi z alpha tan(lambda) curvature
103 *
104 * @param trv The VTrack to propagate.
105 * @param srf The Surface to which to propagate.
106 * @param dir The direction in which to propagate.
107 * @return The propagation status.
108 */
109 public PropStat vecDirProp( VTrack trv, Surface srf,
110 PropDir dir)
111 {
112 TrackDerivative deriv = null;
113 return vecDirProp(trv, srf, dir, deriv);
114
115 }
116
117 /**
118 *Propagate a track without error in the specified direction
119 *and return the derivative matrix in deriv.
120 *
121 * The track parameters for a cylinder are:
122 * phi z alpha tan(lambda) curvature
123 *
124 * @param trv The VTrack to propagate.
125 * @param srf The Surface to which to propagate.
126 * @param dir The direction in which to propagate.
127 * @param deriv The track derivatives to update at the surface srf.
128 * @return The propagation status.
129 */
130 public PropStat vecDirProp( VTrack trv, Surface srf,
131 PropDir dir, TrackDerivative deriv)
132 {
133 PropStat pstat = new PropStat();
134 Propagator.reduceDirection(dir);
135 // Check destination is a ZPlane.
136 Assert.assertTrue( srf.pureType().equals(SurfZPlane.staticType()) );
137 if ( !srf.pureType( ).equals(SurfZPlane.staticType()) )
138 return pstat;
139
140 VTrack trv0 = trv; // want to change this vector
141 TrackDerivative deriv1 = new TrackDerivative();
142 TrackDerivative deriv2 = new TrackDerivative();
143
144 if ( deriv != null ) pstat = vecTransformCylZ( _bfac, trv0, dir, deriv1 );
145 else pstat = vecTransformCylZ( _bfac, trv0, dir, deriv );
146
147 if ( ! pstat.success() ) return pstat;
148
149 if ( deriv != null ) pstat = _propzz.vecDirProp(trv0,srf,dir,deriv2);
150 else pstat = _propzz.vecDirProp(trv0,srf,dir, deriv);
151
152 if ( pstat.success() )
153 {
154 trv = trv0;
155 if ( deriv != null ) deriv.set(deriv2.times(deriv1));
156 }
157
158 return pstat;
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 *Return a String representation of the class' type name.
174 *Included for completeness with the C++ version.
175 *
176 * @return A String representation of the class' type name.
177 */
178 public String type()
179 { return staticType();
180 }
181
182 /**
183 *output stream
184 * @return A String representation of this instance.
185 */
186 public String toString()
187 {
188 return "Cylinder-ZPlane propagation with constant "
189 + bField() + " Tesla field";
190
191 }
192
193 //**********************************************************************
194
195 // Private function to propagate a track without error
196 // The corresponding track parameters are:
197 // On ZPlane:
198 // z (cm) is fixed
199 // 0 - x (cm)
200 // 1 - y (cm)
201 // 2 - dx/dz
202 // 3 - dy/dz
203 // 4 - q/p p is momentum of a track, q is its charge
204 // On Cylinder:
205 // r (cm) is fixed
206 // 0 - phi
207 // 1 - z (cm)
208 // 2 - alpha = phi_dir - phi; tan(alpha) = r*dphi/dr
209 // 3 - sin(lambda) = dz/ds; tan(lambda) = dz/dsT
210 // 4 - q/pT (1/GeV/c) (pT is component of p parallel to cylinder)
211 // If deriviv is nonzero, return the derivative matrix there.
212
213 private PropStat
214 vecTransformCylZ( double B, VTrack trv,
215 PropDir dir,
216 TrackDerivative deriv )
217 {
218
219 // construct return status
220 PropStat pstat = new PropStat();
221
222 // fetch the originating surface and vector
223 Surface srf1 = trv.surface();
224 // TrackVector vec1 = trv.get_vector();
225
226 // Check origin is a Cylinder.
227 Assert.assertTrue( srf1.pureType().equals(SurfCylinder.staticType()) );
228 if ( !srf1.pureType( ).equals(SurfCylinder.staticType()) )
229 return pstat;
230 SurfCylinder scy1 = ( SurfCylinder ) srf1;
231
232
233 // Fetch the R of the cylinder and the starting track vector.
234 int ir = SurfCylinder.RADIUS;
235 double Rcyl = scy1.parameter(ir);
236
237 TrackVector vec = trv.vector();
238 double c1 = vec.get(IPHI); // phi
239 double c2 = vec.get(IZ_CYL); // z
240 double c3 = vec.get(IALF); // alpha
241 double c4 = vec.get(ITLM); // tan(lambda)
242 double c5 = vec.get(IQPT); // q/pt
243
244 // check that dz != 0
245 if(c4 == 0.) return pstat;
246
247 double zpos = c2;
248 double cos_c1 = Math.cos(c1);
249 double sin_c1 = Math.sin(c1);
250 double cos_dir = Math.cos(c1+c3);
251 double sin_dir = Math.sin(c1+c3);
252 double c4_hat2 = 1+c4*c4;
253 double c4_hat = Math.sqrt(c4_hat2);
254
255 double a1 = Rcyl*cos_c1;
256 double a2 = Rcyl*sin_c1;
257 double a3 = cos_dir/c4;
258 double a4 = sin_dir/c4;
259 double a5 = c5/c4_hat;
260
261
262 // if tan(lambda)=dz/dst > 0 : dz>0; dzdst < 0 dz<0
263
264 int sign_dz = 0;
265 if(c4 > 0) sign_dz = 1;
266 if(c4 < 0) sign_dz = -1;
267
268 vec.set(IX , a1);
269 vec.set(IY , a2);
270 vec.set(IDXDZ , a3);
271 vec.set(IDYDZ , a4);
272 vec.set(IQP_Z , a5);
273
274 // Update trv
275 SurfZPlane szp = new SurfZPlane(zpos);
276 trv.setSurface(szp.newPureSurface());
277 trv.setVector(vec);
278
279 // set new direction of the track
280 if(sign_dz == 1) trv.setForward();
281 if(sign_dz == -1) trv.setBackward();
282
283 // Set the return status.
284 pstat.setSame();
285
286 // exit now if user did not ask for error matrix.
287 if ( deriv == null ) return pstat;
288
289 double a34_hat = sign_dz*Math.sqrt(a3*a3 + a4*a4);
290 double a34_hat2 = a34_hat*a34_hat;
291
292 double Rcos_phi = a3/Math.sqrt(1.+a34_hat2)*sign_dz;
293 double Rsin_phi = a4/Math.sqrt(1.+a34_hat2)*sign_dz;
294 // ddphi / dc
295
296 double ddphi_dc2 = -B*a5*Math.sqrt(a34_hat2+1.)*sign_dz;
297 double ddphi_dc2_nob = -Math.sqrt(a34_hat2+1.)*sign_dz;
298
299 // da1 / dc
300
301 double da1_dc1 = -Rcyl*sin_c1;
302
303 // da2 / dc
304
305 double da2_dc1 = Rcyl*cos_c1;
306
307 // da3 / dc
308
309 double da3_dc1 = -sin_dir/c4;
310 double da3_dc3 = -sin_dir/c4;
311 double da3_dc4 = -cos_dir/(c4*c4);
312
313 // da4 / dc
314
315 double da4_dc1 = cos_dir/c4;
316 double da4_dc3 = cos_dir/c4;
317 double da4_dc4 = -sin_dir/(c4*c4);
318
319 // da5 / dc
320
321 double da5_dc4 = -c5*c4/(c4_hat*c4_hat2);
322 double da5_dc5 = 1./c4_hat;
323
324 // da1n/ dc
325
326 double da1n_dc1 = da1_dc1;
327 double da1n_dc2 = Rcos_phi*ddphi_dc2_nob;
328 double da1n_dc3 = 0.;
329 double da1n_dc4 = 0.;
330 double da1n_dc5 = 0.;
331
332 // da2n/ dc
333
334 double da2n_dc1 = da2_dc1;
335 double da2n_dc2 = Rsin_phi*ddphi_dc2_nob;
336 double da2n_dc3 = 0.;
337 double da2n_dc4 = 0.;
338 double da2n_dc5 = 0.;
339
340 // da3n/ dc
341
342 double da3n_dc1 = da3_dc1;
343 double da3n_dc2 = - a4*ddphi_dc2;
344 double da3n_dc3 = da3_dc3;
345 double da3n_dc4 = da3_dc4;
346 double da3n_dc5 = 0.;
347
348 // da4n/ dc
349
350 double da4n_dc1 = da4_dc1;
351 double da4n_dc2 = a3*ddphi_dc2;
352 double da4n_dc3 = da4_dc3;
353 double da4n_dc4 = da4_dc4;
354 double da4n_dc5 = 0.;
355
356 // da5n
357
358 double da5n_dc1 = 0.;
359 double da5n_dc2 = 0.;
360 double da5n_dc3 = 0.;
361 double da5n_dc4 = da5_dc4;
362 double da5n_dc5 = da5_dc5;
363
364 deriv.set(IX,IPHI , da1n_dc1);
365 deriv.set(IX,IZ_CYL , da1n_dc2);
366 deriv.set(IX,IALF , da1n_dc3);
367 deriv.set(IX,ITLM , da1n_dc4);
368 deriv.set(IX,IQPT , da1n_dc5);
369 deriv.set(IY,IPHI , da2n_dc1);
370 deriv.set(IY,IZ_CYL , da2n_dc2);
371 deriv.set(IY,IALF , da2n_dc3);
372 deriv.set(IY,ITLM , da2n_dc4);
373 deriv.set(IY,IQPT , da2n_dc5);
374 deriv.set(IDXDZ,IPHI , da3n_dc1);
375 deriv.set(IDXDZ,IZ_CYL , da3n_dc2);
376 deriv.set(IDXDZ,IALF , da3n_dc3);
377 deriv.set(IDXDZ,ITLM , da3n_dc4);
378 deriv.set(IDXDZ,IQPT , da3n_dc5);
379 deriv.set(IDYDZ,IPHI , da4n_dc1);
380 deriv.set(IDYDZ,IZ_CYL , da4n_dc2);
381 deriv.set(IDYDZ,IALF , da4n_dc3);
382 deriv.set(IDYDZ,ITLM , da4n_dc4);
383 deriv.set(IDYDZ,IQPT , da4n_dc5);
384 deriv.set(IQP_Z,IPHI , da5n_dc1);
385 deriv.set(IQP_Z,IZ_CYL , da5n_dc2);
386 deriv.set(IQP_Z,IALF , da5n_dc3);
387 deriv.set(IQP_Z,ITLM , da5n_dc4);
388 deriv.set(IQP_Z,IQPT , da5n_dc5);
389
390 return pstat;
391 }
392 }