View Javadoc

1   package org.lcsim.recon.tracking.trfxyp;
2   
3   import org.lcsim.recon.tracking.trfutil.TRFMath;
4   import org.lcsim.recon.tracking.trfutil.Assert;
5   
6   import org.lcsim.recon.tracking.trfbase.Propagator;
7   import org.lcsim.recon.tracking.trfbase.PropDirected;
8   import org.lcsim.recon.tracking.trfbase.PropDir;
9   import org.lcsim.recon.tracking.trfbase.PropStat;
10  import org.lcsim.recon.tracking.trfbase.VTrack;
11  import org.lcsim.recon.tracking.trfbase.Surface;
12  import org.lcsim.recon.tracking.trfbase.TrackVector;
13  import org.lcsim.recon.tracking.trfbase.TrackDerivative;
14  
15  /**
16   * Propagates tracks from one XYPlane to another in a constant field
17   * directed along the v axis (parallel to the surfaces of origin and
18   * destination and normal to the z axis).
19   *<p>
20   * Propagation will fail if either the origin or destination is
21   * not an XYPlane or if they are not parallel.
22   *
23   *@author Norman A. Graf
24   *@version 1.0
25   *
26   */
27  public class PropXYXYBV extends PropDirected
28  {
29      
30      // attributes
31      
32      private double _bfac;
33      
34      private static final int    IV = SurfXYPlane.IV;
35      private static final int   IZ   = SurfXYPlane.IZ;
36      private static final int   IDVDU = SurfXYPlane.IDVDU;
37      private static final int   IDZDU = SurfXYPlane.IDZDU;
38      private static final int   IQP  = SurfXYPlane.IQP;
39      
40      // static methods
41      
42      /**
43       *Return a String representation of the class' type name.
44       *Included for completeness with the C++ version.
45       *
46       * @return   A String representation of the class' type name.
47       */
48      public static String typeName()
49      { return "PropXYXYBV"; }
50      
51      /**
52       *Return a String representation of the class' type name.
53       *Included for completeness with the C++ version.
54       *
55       * @return   A String representation of the class' type name.
56       */
57      public static String staticType()
58      { return typeName(); }
59      
60      /**
61       *Construct an instance from a constant solenoidal magnetic field in Tesla.
62       *
63       * @param   bfield The magnetic field strength in Tesla.
64       */
65      public PropXYXYBV(double bfield)
66      {
67          _bfac = TRFMath.BFAC*bfield;
68      }
69      
70      /**
71       *Clone an instance.
72       *
73       * @return A Clone of this instance.
74       */
75      public Propagator newPropagator( )
76      {
77          return new PropXYXYBV( bField() );
78      }
79      
80      /**
81       *Propagate a track without error in the specified direction
82       *and return the derivative matrix in deriv.
83       *
84       * @param   trv The VTrack to propagate.
85       * @param   srf The Surface to which to propagate.
86       * @param   dir The direction in which to propagate.
87       * @param   deriv The track derivatives to update at the surface srf.
88       * @return The propagation status.
89       **/
90      public PropStat vecDirProp( VTrack trv, Surface srf,
91              PropDir dir,TrackDerivative deriv )
92      {
93          PropStat pstat = vec_propagatexyxy_bv_( _bfac, trv, srf, dir, deriv );
94          return pstat;
95      }
96      
97      /**
98       *Propagate a track without error in the specified direction.
99       *
100      * @param   trv The VTrack to propagate.
101      * @param   srf The Surface to which to propagate.
102      * @param   dir The direction in which to propagate.
103      * @return The propagation status.
104      */
105     public PropStat vecDirProp( VTrack trv, Surface srf,
106             PropDir dir)
107     {
108         TrackDerivative deriv = null;
109         return vecDirProp(trv, srf, dir, deriv);
110     }
111     
112     
113     /**
114      *Return the strength of the magnetic field in Tesla.
115      *
116      * @return The strength of the magnetic field in Tesla.
117      */
118     public double bField()
119     {
120         return _bfac/TRFMath.BFAC;
121     }
122     
123     /**
124      *Return a String representation of the class' type name.
125      *Included for completeness with the C++ version.
126      *
127      * @return   A String representation of the class' type name.
128      */
129     public String type()
130     { return staticType(); }
131     
132     
133     /**
134      *output stream
135      *
136      * @return  A String representation of this instance.
137      */
138     public String toString()
139     {
140         return "XYPlane-XYPlane propagation with constant " + bField()
141         + " Tesla field normal to z and parallel to the planes";
142     }
143     
144     
145     
146     // Private function to determine dphi of the propagation
147     double directionBV(int flag_forward, int sign_dphi,
148             double du,
149             double norm, double cat,
150             double sinphi, double cosphi)
151     {
152         
153         int sign_cat = 0;
154         if(du*flag_forward>0.) sign_cat =  1;
155         if(du*flag_forward<0.) sign_cat = -1;
156         if(du == 0.)
157         {
158             if(sinphi >=0. ) sign_cat =  1;
159             if(sinphi < 0. ) sign_cat = -1;
160         }
161         
162         double sin_dphi =  norm*sinphi + cat*cosphi*sign_cat;
163         double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat;
164         
165         int sign_sindphi = 0;
166         if(sin_dphi> 0.) sign_sindphi =  1;
167         if(sin_dphi< 0.) sign_sindphi = -1;
168         if(sin_dphi == 0.) sign_sindphi = sign_dphi;
169         
170         double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi);
171         return dphi;
172         
173     }
174     //**********************************************************************
175     
176     // Private function to propagate a track without error
177     // The corresponding track parameters are:
178     // u (cm) is fixed
179     // 0 - v (cm)
180     // 1 - z (cm)
181     // 2 - dv/du
182     // 3 - dz/du
183     // 4 - q/p   p is momentum of a track, q is its charge
184     // If pderiv is nonzero, return the derivative matrix there.
185     
186     // We use the code from PropXYXY written for the case of B^ || z^.
187     // First, we rotate the track to the c.f. where B^ || z^, propagate the
188     // track and rotate the propagated track back to the original frame.
189     // We require the two planes to be parallel to avoid the corner regions
190     // where the approximation B^ || v^ is not valid.
191     
192     PropStat
193             vec_propagatexyxy_bv_( double B, VTrack trv, Surface srf,
194             PropDir dir,
195             TrackDerivative deriv )
196     {
197         
198         // construct return status
199         PropStat pstat = new PropStat();
200         
201         // fetch the originating surface and vector
202         Surface srf1 = trv.surface();
203         // TrackVector vec1 = trv.get_vector();
204         
205         // Check origin is a XYPlane.
206         Assert.assertTrue( srf1.pureType().equals(SurfXYPlane.staticType()) );
207         if ( !srf1.pureType( ).equals(SurfXYPlane.staticType()) )
208             return pstat;
209         SurfXYPlane sxyp1 = ( SurfXYPlane ) srf1;
210         
211         // Check destination is a XYPlane.
212         Assert.assertTrue( srf.pureType().equals(SurfXYPlane.staticType()) );
213         if ( !srf.pureType( ).equals(SurfXYPlane.staticType()) )
214             return pstat;
215         SurfXYPlane sxyp2 = ( SurfXYPlane ) srf;
216         
217         // If surfaces are the same, we can return now.
218         if ( srf.pureEqual(srf1) )
219         {
220             if(deriv != null)
221             {
222                 deriv.setIdentity();
223             }
224             pstat.setSame();
225             return pstat;
226         }
227         
228         // Fetch the u's and phi's of the planes and the starting track vector.
229         int iphi  = SurfXYPlane.NORMPHI;
230         int idist = SurfXYPlane.DISTNORM;
231         double phi_0 = sxyp1.parameter(iphi);
232         double   u_0 = sxyp1.parameter(idist);
233         double phi_n = sxyp2.parameter(iphi);
234         double   u_n = sxyp2.parameter(idist);
235         double phi_u = phi_0 - phi_n;
236         
237         // check here that first planes are parallel and second that
238         //planes have phi==0||PI , otherwise this propagator doesn't work.
239         
240         if (Math.abs(phi_0)!=0. && Math.abs(phi_0)!=Math.PI ) return pstat;
241         if ( Math.abs(phi_u)!=0. && Math.abs(phi_u)!=Math.PI && Math.abs(phi_u)!=TRFMath.TWOPI ) return pstat;
242         
243         //cout<<" PropXYXYBV "<<endl;
244         //cout<<trv<<endl;
245         TrackVector vec = trv.vector();
246         
247         // b_0 is rotated by 90 deg around the u axis wrt to the input track vector
248         // (v',z') = (-z,v)
249         double b1_0 = -vec.get(IZ);
250         double b2_0 = vec.get(IV);
251         double b3_0 = -vec.get(IDZDU);
252         double b4_0 = vec.get(IDVDU);
253         double b5_0 = vec.get(IQP);
254         
255         
256         double cosphi_u = Math.cos(phi_u);
257         double sinphi_u = Math.sin(phi_u);
258         
259         // check if du == 0 ( that is track moves parallel to the destination plane )
260         double du_du_0 = cosphi_u - b3_0*sinphi_u;
261         if(du_du_0==0.) return pstat;
262         
263         double a_hat_u = 1./du_du_0;
264         double a_hat_u2 = a_hat_u*a_hat_u;
265         
266         double u  = u_0*cosphi_u - b1_0*sinphi_u;
267         double b1 = b1_0*cosphi_u + u_0*sinphi_u;
268         double b2 = b2_0;
269         double b3 = (b3_0*cosphi_u + sinphi_u)*a_hat_u;
270         double b4 = b4_0*a_hat_u;
271         double b5 = b5_0;
272         
273         int sign_du0 = 0;
274         if(trv.isForward())  sign_du0 =  1;
275         if(trv.isBackward()) sign_du0 = -1;
276         if(sign_du0 == 0)
277         {
278             System.out.println("PropXYXYBV._vec_prop: Unknown direction of a track ");
279             System.exit(1);
280         }
281         int sign_du = 0;
282         if(du_du_0*sign_du0 > 0) sign_du =  1;
283         if(du_du_0*sign_du0 < 0) sign_du = -1;
284         
285         // check that q/p != 0
286         Assert.assertTrue(b5 !=0. );
287         
288         // 1/curv of the track is r
289         double r = 1/(b5*B)*Math.sqrt(1 + b3_0*b3_0)/Math.sqrt(1 + b3_0*b3_0 + b4_0*b4_0);
290         double b3_hat = Math.sqrt(1 + b3*b3);
291         double b34_hat = Math.sqrt(1 + b3*b3 + b4*b4);
292         double b3_hat2 = b3_hat*b3_hat;
293         double b34_hat2 = b34_hat*b34_hat;
294         double cosphi =         -b3*sign_du/b3_hat;
295         double sinphi =             sign_du/b3_hat;
296         double rsinphi =  1./(b5*B)*sign_du/b34_hat;
297         double rcosphi = -b3/(b5*B)*sign_du/b34_hat;
298         
299         double du = u_n - u;
300         double norm = du/r - cosphi;
301         
302         // xy-xy propagation failed : noway to the new plane
303         if(Math.abs(norm)>1.) return pstat;
304         
305         double cat = Math.sqrt(1.-norm*norm);
306         int flag_forward = 0;
307         int sign_dphi = 0;
308         
309         if ( dir.equals(PropDir.NEAREST) )
310         {
311             {
312                 // try forward propagation
313                 flag_forward = 1;
314                 if(b5>0.) sign_dphi =  1;
315                 if(b5<0.) sign_dphi = -1;
316                 double dphi1=
317                         directionBV(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi);
318                 // try backward propagation
319                 flag_forward = -1;
320                 if(b5>0.) sign_dphi = -1;
321                 if(b5<0.) sign_dphi =  1;
322                 double dphi2=
323                         directionBV(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi);
324                 if( Math.abs(dphi2)>Math.abs(dphi1) )
325                 {
326                     flag_forward = -flag_forward;
327                     sign_dphi = - sign_dphi;
328                 }
329             }
330         }
331         else if ( dir.equals(PropDir.FORWARD) )
332         {
333             flag_forward = 1;
334             if(b5>0.) sign_dphi =  1;
335             if(b5<0.) sign_dphi = -1;
336         }
337         else if ( dir.equals(PropDir.BACKWARD) )
338         {
339             flag_forward = -1;
340             if(b5>0.) sign_dphi = -1;
341             if(b5<0.) sign_dphi =  1;
342         }
343         else
344         {
345             System.out.println( "PropXYXY._vec_prop: Unknown direction.");
346             System.exit(1);
347         }
348         
349         int sign_cat = 0;
350         if(du*sign_dphi*b5>0.) sign_cat =  1;
351         if(du*sign_dphi*b5<0.) sign_cat = -1;
352         if(du == 0.)
353         {
354             if(sinphi >=0. ) sign_cat =  1;
355             if(sinphi < 0. ) sign_cat = -1;
356         }
357         
358         double sin_dphi =  norm*sinphi + cat*cosphi*sign_cat;
359         double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat;
360         if(cos_dphi>1.) cos_dphi=1.;
361         if(cos_dphi<-1.) cos_dphi= -1.;
362         
363         int sign_sindphi = 0;
364         if(sin_dphi> 0.) sign_sindphi =  1;
365         if(sin_dphi< 0.) sign_sindphi = -1;
366         if(sin_dphi == 0.) sign_sindphi = sign_dphi;
367         
368         double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi);
369         
370         // check that I didn't make any mistakes
371         
372         Assert.assertTrue(Math.abs(Math.sin(dphi) -sin_dphi)<1.e-5);
373         
374         // check if dun == 0 ( that is track moves parallel to the destination plane)
375         double du_n_du = cos_dphi - b3*sin_dphi;
376         if(du_n_du==0.) return pstat;
377         
378         double a_hat_dphi = 1./du_n_du;
379         double a_hat_dphi2 = a_hat_dphi*a_hat_dphi;
380         double c_hat_dphi =     sin_dphi + b3*cos_dphi ;
381         
382         double b1_n = b1 + rsinphi*(1-cos_dphi) - rcosphi*sin_dphi;
383         double b2_n = b2 + b4/(b5*B)*sign_du/b34_hat*dphi;
384         double b3_n = c_hat_dphi*a_hat_dphi;
385         double b4_n =         b4*a_hat_dphi;
386         double b5_n = b5;
387         
388         // double u_n_0 = u_n*cosphi_u + b1_n*sinphi_u;
389         
390         // check if track crossed original plane during the propagation
391         // switch (dir) {
392         //  if( dir.equals(PropDir.FORWARD:
393         //   if((u_n_0 - u_0)*sign_du0<0) return pstat;
394         //  break;
395         //  if( dir.equals(PropDir.BACKWARD:
396         //   if((u_n_0 - u_0)*sign_du0>0) return pstat;
397         //  break;
398         // }
399         
400         int sign_dun = 0;
401         if(du_n_du*sign_du > 0) sign_dun =  1;
402         if(du_n_du*sign_du < 0) sign_dun = -1;
403         
404         // rotate the propagated track  back by -90 deg around the u axis:
405         // (v,z) = (z',-v').
406         vec.set(IV    , b2_n);
407         vec.set(IZ    , -b1_n);
408         vec.set(IDVDU , b4_n);
409         vec.set(IDZDU , -b3_n);
410         vec.set(IQP   , b5_n);
411         
412         // Update trv
413         trv.setSurface( srf.newPureSurface() );
414         trv.setVector(vec);
415         
416         // set new direction of the track
417         if(sign_dun ==  1) trv.setForward();
418         if(sign_dun == -1) trv.setBackward();
419         //cout<<trv<<endl;
420         // Calculate sT.
421         double crv = B*b5;
422         double dv = b1_n - b1;
423         double dxy = Math.sqrt( du*du + dv*dv );
424         double arg = 0.5*crv*dxy;
425         double dst = dxy*TRFMath.asinrat(arg);
426         
427         // Calculate s.
428         double dz = b2_n - b2;
429         Assert.assertTrue( flag_forward==1 || flag_forward==-1 );
430         double ds = ((double) flag_forward)*Math.sqrt(dst*dst+dz*dz);
431         
432         // Set the return status.
433         pstat.setPathDistance(ds);
434         //(flag_forward==1)?pstat.set_forward():pstat.set_backward();
435         
436         // exit now if user did not ask for error matrix.
437         if (  deriv == null  ) return pstat;
438         
439         // du_d0
440         
441         double du_db1_0 = - sinphi_u;
442         
443         // db1_d0
444         
445         double db1_db1_0 = cosphi_u;
446         
447         // db3_d0
448         
449         double db3_db3_0 = a_hat_u2;
450         
451         // db4_d0
452         
453         double db4_db3_0 = b4_0*sinphi_u*a_hat_u2;
454         double db4_db4_0 = a_hat_u;
455         
456         // dr_d
457         
458         double dr_db3 = r*b3*b4*b4/(b3_hat2*b34_hat2);
459         double dr_db4 = -r*b4/b34_hat2;
460         double dr_db5 = -r/b5;
461         
462         // dcosphi_d
463         
464         double dcosphi_db3 = - sign_du/b3_hat - cosphi*b3/b3_hat2;
465         
466         // dsinphi_d
467         
468         double dsinphi_db3 = - sinphi*b3/b3_hat2;
469         
470         // dcat_d
471         
472         double dcat_db3 = norm/cat*(du/(r*r)*dr_db3 + dcosphi_db3 );
473         double dcat_db4 = norm/cat* du/(r*r)*dr_db4;
474         double dcat_db5 = norm/cat* du/(r*r)*dr_db5;
475         double dcat_du  = norm/(cat*r);
476         
477         // dnorm_d
478         
479         double dnorm_db3 = - du/(r*r)*dr_db3 - dcosphi_db3;
480         double dnorm_db4 = - du/(r*r)*dr_db4;
481         double dnorm_db5 = - du/(r*r)*dr_db5;
482         double dnorm_du  = - 1./r;
483         
484         // dcos_dphi_d
485         
486         double dcos_dphi_db3 = - cosphi*dnorm_db3 - norm*dcosphi_db3 +
487                 sign_cat*(sinphi*dcat_db3 + cat*dsinphi_db3);
488         double dcos_dphi_db4 = - cosphi*dnorm_db4 + sign_cat*sinphi*dcat_db4;
489         double dcos_dphi_db5 = - cosphi*dnorm_db5 + sign_cat*sinphi*dcat_db5;
490         double dcos_dphi_du  = - cosphi*dnorm_du  + sign_cat*sinphi*dcat_du;
491         
492         // dsin_dphi_d
493         
494         double dsin_dphi_db3 = sinphi*dnorm_db3 + norm*dsinphi_db3 +
495                 sign_cat*(cosphi*dcat_db3 + cat*dcosphi_db3);
496         double dsin_dphi_db4 = sinphi*dnorm_db4 + sign_cat*cosphi*dcat_db4;
497         double dsin_dphi_db5 = sinphi*dnorm_db5 + sign_cat*cosphi*dcat_db5;
498         double dsin_dphi_du  = sinphi*dnorm_du  + sign_cat*cosphi*dcat_du;
499         
500         // ddphi_d
501         
502         double ddphi_db3;
503         double ddphi_db4;
504         double ddphi_db5;
505         double ddphi_du;
506         if(Math.abs(sin_dphi)>0.5)
507         {
508             ddphi_db3 = - dcos_dphi_db3/sin_dphi;
509             ddphi_db4 = - dcos_dphi_db4/sin_dphi;
510             ddphi_db5 = - dcos_dphi_db5/sin_dphi;
511             ddphi_du  = - dcos_dphi_du /sin_dphi;
512         }
513         else
514         {
515             ddphi_db3 = dsin_dphi_db3/cos_dphi;
516             ddphi_db4 = dsin_dphi_db4/cos_dphi;
517             ddphi_db5 = dsin_dphi_db5/cos_dphi;
518             ddphi_du  = dsin_dphi_du /cos_dphi;
519         }
520         
521         // da_hat_dphi_d
522         
523         double da_hat_dphi_db3 = - a_hat_dphi2*
524                 (dcos_dphi_db3 - sin_dphi - b3*dsin_dphi_db3);
525         double da_hat_dphi_db4 = - a_hat_dphi2*(dcos_dphi_db4 - b3*dsin_dphi_db4);
526         double da_hat_dphi_db5 = - a_hat_dphi2*(dcos_dphi_db5 - b3*dsin_dphi_db5);
527         double da_hat_dphi_du  = - a_hat_dphi2*(dcos_dphi_du  - b3*dsin_dphi_du );
528         
529         // dc_hat_dphi_d
530         
531         double dc_hat_dphi_db3 = b3*dcos_dphi_db3 + dsin_dphi_db3 + cos_dphi;
532         double dc_hat_dphi_db4 = b3*dcos_dphi_db4 + dsin_dphi_db4;
533         double dc_hat_dphi_db5 = b3*dcos_dphi_db5 + dsin_dphi_db5;
534         double dc_hat_dphi_du  = b3*dcos_dphi_du  + dsin_dphi_du ;
535         
536         // db1_n_d
537         
538         double db1_n_db1 = 1;
539         double db1_n_db3 = (dr_db3*sinphi+r*dsinphi_db3)*(1-cos_dphi)
540         - rsinphi*dcos_dphi_db3
541                 - dr_db3*cosphi*sin_dphi - r*dcosphi_db3*sin_dphi
542                 - rcosphi*dsin_dphi_db3;
543         double db1_n_db4 = dr_db4*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db4
544                 - dr_db4*cosphi*sin_dphi - rcosphi*dsin_dphi_db4;
545         double db1_n_db5 = dr_db5*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db5
546                 - dr_db5*cosphi*sin_dphi - rcosphi*dsin_dphi_db5;
547         double db1_n_du  = - rsinphi*dcos_dphi_du - rcosphi*dsin_dphi_du;
548         
549         // db2_n_d
550         
551         double db2_n_db2 = 1.;
552         double db2_n_db3 = 1./(b5*B)*b4*sign_du/b34_hat*
553                 ( - dphi*b3/b34_hat2 + ddphi_db3);
554         double db2_n_db4 = 1./(b5*B)*sign_du/b34_hat*
555                 ( - dphi*b4*b4/b34_hat2 + b4*ddphi_db4 + dphi );
556         double db2_n_db5 = 1./(b5*B)*b4*sign_du/b34_hat*( ddphi_db5 - dphi/b5);
557         double db2_n_du  = 1./(b5*B)*b4*sign_du/b34_hat*ddphi_du;
558         
559         // db3_n_d
560         
561         double db3_n_db3 = a_hat_dphi*dc_hat_dphi_db3 + da_hat_dphi_db3*c_hat_dphi;
562         double db3_n_db4 = a_hat_dphi*dc_hat_dphi_db4 + da_hat_dphi_db4*c_hat_dphi;
563         double db3_n_db5 = a_hat_dphi*dc_hat_dphi_db5 + da_hat_dphi_db5*c_hat_dphi;
564         double db3_n_du  = a_hat_dphi*dc_hat_dphi_du  + da_hat_dphi_du *c_hat_dphi;
565         
566         // db4_n_d
567         
568         double db4_n_db3 = b4*da_hat_dphi_db3;
569         double db4_n_db4 = b4*da_hat_dphi_db4 + a_hat_dphi;
570         double db4_n_db5 = b4*da_hat_dphi_db5;
571         double db4_n_du  = b4*da_hat_dphi_du;
572         
573         // db5_n_d
574         
575         // db1_n_d0
576         
577         double db1_n_db1_0 = db1_n_du * du_db1_0 + db1_n_db1*db1_db1_0;
578         double db1_n_db2_0 = 0.;
579         double db1_n_db3_0 = db1_n_db3*db3_db3_0 + db1_n_db4*db4_db3_0;
580         double db1_n_db4_0 = db1_n_db4*db4_db4_0;
581         double db1_n_db5_0 = db1_n_db5;
582         
583         // db2_n_d0
584         
585         double db2_n_db1_0 = db2_n_du *du_db1_0;
586         double db2_n_db2_0 = db2_n_db2;
587         double db2_n_db3_0 = db2_n_db3*db3_db3_0 + db2_n_db4*db4_db3_0;
588         double db2_n_db4_0 = db2_n_db4*db4_db4_0;
589         double db2_n_db5_0 = db2_n_db5;
590         
591         // db3_n_d0
592         
593         double db3_n_db1_0 = db3_n_du*du_db1_0;
594         double db3_n_db2_0 = 0.;
595         double db3_n_db3_0 = db3_n_db3*db3_db3_0 + db3_n_db4*db4_db3_0;
596         double db3_n_db4_0 = db3_n_db4*db4_db4_0;
597         double db3_n_db5_0 = db3_n_db5;
598         
599         // db4_n_d0
600         
601         double db4_n_db1_0 = db4_n_du*du_db1_0;
602         double db4_n_db2_0 = 0.;
603         double db4_n_db3_0 = db4_n_db3*db3_db3_0 + db4_n_db4*db4_db3_0;
604         double db4_n_db4_0 = db4_n_db4*db4_db4_0;
605         double db4_n_db5_0 = db4_n_db5;
606         
607         // db5_n_d0
608         
609         double db5_n_db1_0 = 0.;
610         double db5_n_db2_0 = 0.;
611         double db5_n_db3_0 = 0.;
612         double db5_n_db4_0 = 0.;
613         double db5_n_db5_0 = 1.;
614         
615         // apply the rotation to the matrix db to return to the original c.f.
616         deriv.set(IV,IV,        db2_n_db2_0);
617         deriv.set(IV,IZ,        -db2_n_db1_0);
618         deriv.set(IV,IDVDU,     db2_n_db4_0);
619         deriv.set(IV,IDZDU,     -db2_n_db3_0);
620         deriv.set(IV,IQP,       db2_n_db5_0);
621         deriv.set(IZ,IV,        -db1_n_db2_0);
622         deriv.set(IZ,IZ,        db1_n_db1_0);
623         deriv.set(IZ,IDVDU,     -db1_n_db4_0);
624         deriv.set(IZ,IDZDU,     db1_n_db3_0);
625         deriv.set(IZ,IQP,        -db1_n_db5_0);
626         deriv.set(IDVDU,IV,     db4_n_db2_0);
627         deriv.set(IDVDU,IZ,     -db4_n_db1_0);
628         deriv.set(IDVDU,IDVDU,  db4_n_db4_0);
629         deriv.set(IDVDU,IDZDU,  -db4_n_db3_0);
630         deriv.set(IDVDU,IQP,    db4_n_db5_0);
631         deriv.set(IDZDU,IV,     -db3_n_db2_0);
632         deriv.set(IDZDU,IZ,     db3_n_db1_0);
633         deriv.set(IDZDU,IDVDU,  -db3_n_db4_0);
634         deriv.set(IDZDU,IDZDU,  db3_n_db3_0);
635         deriv.set(IDZDU,IQP,    -db3_n_db5_0);
636         deriv.set(IQP,IV,       db5_n_db2_0);
637         deriv.set(IQP,IZ,       -db5_n_db1_0);
638         deriv.set(IQP,IDVDU,    db5_n_db4_0);
639         deriv.set(IQP,IDZDU,     -db5_n_db3_0);
640         deriv.set(IQP,IQP,      db5_n_db5_0);
641         
642         return pstat;
643         
644     }
645     
646 }