*-- Author : S.Burke / J.V. Morris
SUBROUTINE FKTRAN(DZ,Z1,S1,S2,D)
**********************************************************************
* *
* Transform track vector (S1) through DZ to (S2) in magnetic field. *
* *
**********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*KEEP,FKNPL.
CHARACTER*5 CKDBG
PARAMETER (CKDBG='FKDBG')
PARAMETER (NPL=72)
LOGICAL LTRUE,LFIRST,LTRPL,LTRPLD
DOUBLE PRECISION TRUE,RTRUE,CHITRU,SPRO,CPRO,SFIL,CFIL
&, SSMT,CSMT,SSMTR,CSMTR,DPRO,CBPRO,QPRO,QGAIN
&, RPRO,CRPRO,RFIL,CRFIL,RSMT,CRSMT,CHIFIL,CHISMT
*
* Per-track values can go in H1WORK; note that LTRUE and LFIRST must
* be set at least per event.
*
* This is about 36k words long; the remaining common blocks are
* about 3.6k in total. Some of this could be in /H1WORK/, but the
* blocks would have to be reorganised.
*
COMMON /H1WORK/
* /FKPROJ/
& SPRO(5,NPL),CPRO(5,5,NPL)
* /FKFILT/
&, SFIL(5,NPL),CFIL(5,5,NPL)
* /FKSMTH/
&, SSMT(5,NPL),CSMT(5,5,NPL)
&, SSMTR(5,NPL),CSMTR(5,5,NPL)
* /FKINT/
&, DPRO(5,5,NPL),CBPRO(5,5,NPL),QPRO(5,5,NPL)
&, QGAIN(5,5,NPL),IAPROX,LFIRST
* /FKRSID/
&, RPRO(2,NPL),CRPRO(2,2,NPL),RFIL(2,NPL)
&, CRFIL(2,2,NPL),RSMT(2,NPL),CRSMT(2,2,NPL)
&, CHIFIL(NPL),CHISMT(NPL)
* /FKTRUE/
&, TRUE(5,NPL),RTRUE(5,NPL),CHITRU(NPL),LTRUE
* /FKDBG/
&, LTRPL(NPL),LTRPLD(NPL)
*KEEP,FKINT.
*KEND.
REAL B(3),R(3),B2(3),R2(3),DB(3)
DIMENSION S1(5),S2(5),D(5,5)
*KEEP,FKPIDP.
DOUBLE PRECISION PI,TWOPI,PIBY2
PARAMETER (PI=3.141592653589793238)
PARAMETER (TWOPI=PI*2.0D0,PIBY2=PI/2.0D0)
*KEND.
PARAMETER (BFACT=2.997886E-4)
**********************************************************************
*
* Thresholds for various approximations (*** UNOPTIMISED ***)
*
* Max z step over which field is assumed constant
PARAMETER (DZMAX=1.0D0)
* If Dphi CALL GUFLD(R,B)
IF (ABS(B(3)).LT.0.01) B(3) = SIGN(0.01,B(3))
* Work out Delta phi (don't re-order these lines!)
D(5,3) = -DZSECT*B(3)*BFACT
DPHI = D(5,3)*S1(3)
D(5,4) = DPHI*STCT
S2(5) = S2(5) + DPHI
* Choose the order of approximation
IF (DABS(DZ).GT.DZMAX) THEN
* Must consider change in field over DZ
IAPROX = 5
ELSEIF (DABS(DZ*DPHI).GT.DXYMIN) THEN
* Must take helix (x and y affected by dphi)
IAPROX = 4
ENDIF
ELSE
* Ignore DPHI
D(5,3) = 0.D0
D(5,4) = 0.D0
ENDIF
IF (IAPROX.LE.2) THEN
* Ignore change in phi for DX and DY
D(1,3) = 0.D0
D(2,3) = 0.D0
D(1,4) = DZ*CPHI1
D(2,4) = DZ*SPHI1
D(1,5) = -D(2,4)*S1(4)
D(2,5) = D(1,4)*S1(4)
S2(1) = S2(1) + D(2,5)
S2(2) = S2(2) - D(1,5)
* End of straight line approx
IF (IAPROX.EQ.2) GOTO 9100
GOTO 9200
ENDIF
IF (DABS(DPHI).LT.DPHMAX) THEN
* Small angle approximations for cos and sin (to order DPHI**2)
DPHI2 = 0.5D0*DPHI
CDPHI = 1.D0 - DPHI2*DPHI
SPHI2 = SPHI1*CDPHI + CPHI1*DPHI
CPHI2 = CPHI1*CDPHI - SPHI1*DPHI
DXOTT = DZ*(CPHI1 - SPHI1*DPHI2)
DYOTT = DZ*(SPHI1 + CPHI1*DPHI2)
ELSE
* Full expressions
CPHI2 = DCOS(S2(5))
SPHI2 = DSIN(S2(5))
FAC = DZ/DPHI
DXOTT = FAC*(SPHI2 - SPHI1)
DYOTT = FAC*(CPHI1 - CPHI2)
* Remember values
PHIR = S2(5)
CPHIR = CPHI2
SPHIR = SPHI2
ENDIF
* Change in x and y
DX = DXOTT*S1(4)
DY = DYOTT*S1(4)
S2(1) = S2(1) + DX
S2(2) = S2(2) + DY
D(1,5) = -DY
D(2,5) = DX
* Some more abbreviations ...
IF (DABS(S1(3)).GT.1.0D-10) THEN
POQ = 1.D0/S1(3)
ELSE
POQ = SIGN(1.0D10,S1(3))
ENDIF
DXDY = DZ*S1(4)
DXDYC2 = DXDY*CPHI2
DXDYS2 = DXDY*SPHI2
* ... and the remaining differential coefficients
D(1,3) = POQ*(DXDYC2 - DX)
D(1,4) = DXOTT*CTSQ + DXDYC2*STCT
D(2,3) = POQ*(DXDYS2 - DY)
D(2,4) = DYOTT*CTSQ + DXDYS2*STCT
* Are corrections for Bx and By required?
IF (IAPROX.LT.4) GOTO 9100
*
* Apply 1st order corrections assuming that Bx and By are
* small but non-zero ........ (Bx,By,Bz) is still assumed
* to be constant across DZ ......
*
DPHIOB = DZSECT*S1(3)*BFACT
CON1 = DPHIOB*DZ*0.5D0
CON2 = DPHIOB*OPTSQ
IF (S1(4).NE.0.D0) THEN
CON3 = DPHIOB/S1(4)
ELSE
* If theta = 0, phi is undefined, so CON3 can be set to any value
CON3 = 0.D0
ENDIF
S2(1) = S2(1) - CON1*B(2)
S2(2) = S2(2) + CON1*B(1)
S2(4) = S2(4) + CON2*(B(1)*SPHI1 - B(2)*CPHI1)
S2(5) = S2(5) + CON3*(B(1)*CPHI1 + B(2)*SPHI1)
* Are corrections for field variation required?
IF (IAPROX.LT.5) GOTO 9000
* Get field at end point ...
R2(1) = S2(1)
R2(2) = S2(2)
R2(3) = Z1 + DZ
CALL GUFLD(R2,B2)
IF (ABS(B2(3)).LT.0.01) B2(3) = SIGN(0.01,B2(3))
DB(1) = 0.5*(B2(1) - B(1))
DB(2) = 0.5*(B2(2) - B(2))
DB(3) = 0.5*(B2(3) - B(3))
*
* DB is the mean difference from the assumed field along the trajectory,
* assuming the field components vary linearly with Z over DZ.
* ie B and dB/dZ are non-zero but d2B/dZ2 etc are assumed zero.
*
DBZTT = DB(3)*S1(4)
S2(1) = S2(1) + CON1*(DBZTT*SPHI1 - DB(2))
S2(2) = S2(2) + CON1*(DB(1) - DBZTT*CPHI1)
S2(4) = S2(4) + CON2*(DB(1)*SPHI1 - DB(2)*CPHI1)
S2(5) = S2(5) + CON3*(DB(1)*CPHI1 + DB(2)*SPHI1) - DPHIOB*DB(3)
9000 CONTINUE
IF (S2(4).LT.0.D0) THEN
S2(4) = DABS(S2(4))
S2(5) = S2(5) + PI
ENDIF
9100 CONTINUE
IF (DABS(S2(5)).GT.TWOPI) S2(5)=DMOD(S2(5),TWOPI)
IF (S2(5).LT.0.D0) S2(5) = S2(5) + TWOPI
9200 CONTINUE
* Stop values getting too large (i.e. stop overflows)
IF (ABS(S2(1)).GT.1.D4) S2(1) = SIGN(1.D4,S2(1))
IF (ABS(S2(2)).GT.1.D4) S2(2) = SIGN(1.D4,S2(2))
IF (S2(4).GT.1.D6) S2(4) = 1.D6
RETURN
END
*