C======================================================================== ====== SUBROUTINE CTVMVF (PRINT, Nv, VXI) C======================================================================== ====== C===Description: C Calculates contributions from track Nt (at vertex Nv) to the derivative C vector VXI and derivitive matrix VMAT for the vertex fit section of CTVMFT. C This portion of the code is somewhat awkward to follow. See the derivation C of these formulae in CDF note 19zz (J.P. Marriner) C NB: this code does not properly handle tracks that turn through more than C 90 degrees between the distance of closest approach to the z axis and C the vertex. These tracks are considered pathological. The code will C be correct everywhere, however, if the proper branch of the arcsin C function is taken. (The initial approximation may be so bad that the C desired solution will not be found, however). C===Input Arguments: C Nv vertex C===Output Arguments: C VXI first derivative contributions for this track to the vertex fit C VMAT second derivative matrix (in common, include file CTVMFT.INC) C------------------------------------------------------------------------ ------ C===Implicit None Declaration: IMPLICIT NONE C===Global Declarations: C---- the following include should be used for all definitions of pi REAL PI C value of PI PARAMETER (PI=3.141592653) REAL TWOPI C 2pi PARAMETER (TWOPI=2.0*PI) REAL HALFPI C pi/2 PARAMETER (HALFPI=0.5*PI) REAL RADDEG C conversion from radians to degrees PARAMETER (RADDEG= 180./PI) REAL DEGRAD C conversion from degrees to radians PARAMETER (DEGRAD= PI/180.) C---- the above include should be used for all definitions of pi #include C===Local Declarations: INTEGER PRINT INTEGER NV DOUBLE PRECISION VXI(MAXDIM) INTEGER I,J,K,L INTEGER NT, NtF,NvF, TI,TJ, VI,Vj REAL C,PHI,COTH, D,Z, PT REAL TWOC, S,TS,TRAD REAL SP,CP, XCYS,XSYC,XC2YS2, SINCS,COSCS,TWOCS REAL DPXDT, DPYDT, DSDC,DSDP,DSDX,DSDY REAL EMAT(3,2), FMAT(3,2) C------------------------------------------------------------------------ ------ C===Start of Code: DO 50 Nt=1,NTRACK IF (.NOT.TrkVtx(Nt,Nv)) GO TO 50 C current half curvature C = PAR0(1,NT) + PARDIF(1,NT) C current phi PHI = PAR0(2,NT) + PARDIF(2,NT) C current coth(theta) COTH = PAR0(3,NT) + PARDIF(3,NT) TWOC = 2.0*C C Xsv*cos(phi) + Ysv*sin(phi), -Xsv*sin(phi) + Ysv*cos(phi) C (Xsv*cos(phi))**2 + (Ysv*sin(phi))**2 SP = SIN(PHI) CP = COS(PHI) XCYS = XYZVRT(1,Nv)*CP + XYZVRT(2,Nv)*SP XSYC = -XYZVRT(1,Nv)*SP + XYZVRT(2,Nv)*CP XC2YS2 = (XYZVRT(1,Nv)*CP)**2 + (XYZVRT(2,Nv)*SP)**2 c----------------------------------------------------------------------- C t = 2.0 * c * (Xsv*cos(phi) + Ysv*sin(phi)) = argument of arcsin C s = projected distance in xy plane along helix to vertex C sin ( c*s), cos ( c*s ), 2.0 * sin(c*s) * cos(c*s) C d(arcsin(t))/dt = 1 / sqrt (1 - t**2 ) c----------------------------------------------------------------------- TS = TWOC*XCYS S = ASIN(TS)/TWOC SINCS = SIN(C*S) COSCS = COS(C*S) TWOCS = 2.0*SINCS*COSCS TRAD = SQRT((1.-TS)*(1.+TS)) c----------------------------------------------------------------------- C Helix parameters d0 and z0 are dependent variables. C They are functions of the vertex and track parameters. c----------------------------------------------------------------------- D = XSYC - SINCS**2/C Z = XYZVRT(3,Nv) - COTH*S C Difference between fitted and measured d0, z0 PARDIF(4,NT) = D - PAR0(4,NT) PARDIF(5,NT) = Z - PAR0(5,NT) C Remember all flavors of momentum for this track PT = PSCALE/ABS(C) C x component of momentum TrkP4(NT,1) = PT*(CP*TRAD-SP*TS) C y component of momentum TrkP4(NT,2) = PT*(SP*TRAD+CP*TS) C z component of momentum TrkP4(NT,3) = PT*COTH C Transverse momentum TrkP4(NT,5) = PT C total momentum TrkP4(NT,6) = PT*SQRT(1.0+COTH**2) C energy TrkP4(NT,4) = SQRT(TrkP4(NT,6)**2+TMASS(NT)**2) C dPx/dt DPXDT =-PT*(SP+TS*CP/TRAD) C dPy/dt DPYDT = PT*(CP-TS*SP/TRAD) C dPx/dc DDA(NT,1) =-TrkP4(NT,1)/C + 2.0*XCYS*DPXDT C dPx/dPhi DDA(NT,2) =-PT*(TRAD*SP+TS*CP) + TWOC*XSYC*DPXDT C dPx/dXsv DDA(NT,3) = TWOC*CP*DPXDT C dPx/dYsv DDA(NT,4) = TWOC*SP*DPXDT C dPy/dc DDA(NT,5) =-TrkP4(NT,2)/C + 2.0*XCYS*DPYDT C dPy/dPhi DDA(NT,6) = PT*(TRAD*CP-TS*SP) + TWOC*XSYC*DPYDT C dPy/dXsv DDA(NT,7) = TWOC*CP*DPYDT C dPy/dYsv DDA(NT,8) = TWOC*SP*DPYDT c----------------------------------------------------------------------- C Get the derivatives of s with respect to the fit parameters C See if argument of arcsin is small C (track with a small turning angle to the vertex point) C Yes. use power series approximation c----------------------------------------------------------------------- IF (ABS(TS).LT.5.E-3) THEN DSDC =-(0.666667 - 0.3*TS**2)*(TS*TS*TS) C No. Use full functional form ELSE DSDC = TS/TRAD - TWOC*S ENDIF C Get full derivative ds/dc DSDC = DSDC/(TWOC*C) C ds/dphi DSDP = XSYC/TRAD C ds/dXsv DSDX = CP/TRAD C ds/dYsv DSDY = SP/TRAD C d(d0)/dc, d(d0)/dphi, d(d0)/dcot(theta) EMAT(1,1) = -SINCS*(TWOC*S*COSCS-SINCS)/C**2 - TWOCS*DSDC EMAT(2,1) = -XCYS - TWOCS*DSDP EMAT(3,1) = 0.0 C d(z0)/dc, d(z0)/dphi, d(z0)/dcot(theta) EMAT(1,2) = -COTH*DSDC EMAT(2,2) = -COTH*DSDP EMAT(3,2) = -S C d(d0)/dXsv, d(d0)/dYsv, d(d0)/dZsv FMAT(1,1) = -SP - TWOCS*CP/TRAD FMAT(2,1) = CP - TWOCS*SP/TRAD FMAT(3,1) = 0.0 C d(z0)/dXsv, d(z0)/dYsv, d(z0)/dZsv FMAT(1,2) = -COTH*DSDX FMAT(2,2) = -COTH*DSDY FMAT(3,2) = 1.0 C Index into matrix equations for track NT NtF = TOFF(NT) C Index into matrix equations for vertex NV NvF = VOFF(Nv) DO 20 I=1,3 TI = NtF + I VI = NvF + I C loop over (d0,z0) DO K=1,2 C Ft * G2 * Xi2 VXI(Vi) = VXI(Vi) - FMAT(I,K)* (G(K+3,4,NT)*PARDIF(4,NT) & + G(K+3,5,NT)*PARDIF(5,NT)) C (Gd + Et * G2 ) * Xi2 VXI(Ti) = VXI(Ti) - PARDIF(K+3,NT) * & (G(I,K+3,NT) + EMAT(I,1)*G(4,K+3,NT) + EMAT(I,2)*G(5,K+3,NT)) END DO DO 10 J=1,3 TJ = NtF + J VJ = NvF + J C Ft * Gd * Xi3 VXI(Vi) = VXI(Vi) & - (FMAT(I,1)*G(4,J,NT)+FMAT(I,2)*G(5,J,NT))*PARDIF(J,NT) C (G3 + Et * Gd) * Xi3 VXI(Ti) = VXI(Ti) - (G(I,J,NT) & + EMAT(I,1)*G(4,J,NT)+EMAT(I,2)*G(5,J,NT))*PARDIF(J,NT) C G3 VMAT(Ti,Tj) = G(I,J,NT) DO K=1,2 C Gd * F VMAT(Vi,Tj) = VMAT(Vi,Tj)+FMAT(I,K)*G(K+3,J,NT) C Et * Gdt + Gd * E VMAT(Ti,Tj) = VMAT(Ti,Tj) & + EMAT(I,K)*G(K+3,J,NT)+ G(I,K+3,NT)*EMAT(J,K) DO L=1,2 C Ft * G2 * F VMAT(Vi,Vj) = VMAT(Vi,Vj) & + FMAT(I,K)*G(K+3,L+3,NT)*FMAT(J,L) C Et * G2 * F VMAT(Vi,Tj) = VMAT(Vi,Tj) & + FMAT(I,K)*G(K+3,L+3,NT)*EMAT(J,L) C Et * G2 * E VMAT(Ti,Tj) = VMAT(Ti,Tj) & + EMAT(I,K)*G(K+3,L+3,NT)*EMAT(J,L) C end of loop on L END DO C end of loop on K END DO C end of loop on J 10 CONTINUE C end of loop on I 20 CONTINUE c$$$D IF (PRINT.GT.0) THEN c$$$D TS = TS/C c$$$D WRITE(PRINT,1045) NT,(PARDIF(I,NT),I=1,5), D,Z,TS c$$$D1045 FORMAT(/,' Track',I3,1X,1P,5E11.3,3X,0P,3F11.4) c$$$D WRITE(PRINT,1047) (DDA(Nt,J),J=1,8) c$$$D1047 FORMAT(10X,1P,8E11.3) c$$$D WRITE(PRINT,1046) (VXI(I),I=1,MATDIM) c$$$D DO I=1,MATDIM c$$$D WRITE(PRINT,1046) (VMAT(J,I),J=1,I) c$$$D1046 FORMAT(5X,0P,10E11.3) c$$$D END DO c$$$D END IF C end of the track loop, this vertex 50 CONTINUE C------------------------------------------------------------------------ ------ C===Return to Caller: C Symmetrize the derivative matrix DO I=1,MATDIM-1 DO J=I+1,MATDIM VMAT(J,I) = VMAT(I,J) ENDDO ENDDO RETURN END