C=============================================================================== SUBROUTINE CTVMFA (Nv, PRINT,IERR) C=============================================================================== C C Calculates the initial approximation to the vertex position. C This subroutine tries to intersect circles that are (R,Phi) projections C of the helical trajectories of the two tracks specified in the calling C sequence. C If the circles intersect, the intersection with smaller z difference C between the 2 helices extrapolated to the intersection is chosen. C If the circles do not intersect, the vertex approximation is taken as C the point of closest appraoch of the two circles. Both cases of non- C intersection ("interiour", "exteriour") are treated. C Note that the tracks are required to travel "forwards" to the vertex, C and this requirement may override the "smallest delta Z" criterion. C Tests are made on the vertex radius and Z difference between the two C tracks, and error codes are returned for unacceptable vertices. C Note that, where possible, the calculation is completed, even if an C error condition is flagged. C Test parameter settings used to accept/reject vertex approximations; C DRMAX,DZMAX,RVMAX, TRNMAX,DSMIN C appear as data statements in the include file CTVMFT.INC. C IERR = 2 = IJKERR(1) marks failure in vertex first approximation C IJKERR(2)= 1 Concentric circles, in the first approximation C 2 (conversion) widely separated exteriour circles at the midpoint C 3 (conversion) widely separated interiour circles at the midpoint C 4 widely separated exteriour circles at the approximate vertex C 5 widely separated interiour circles at the approximate vertex C 6 Rv is too large at the chosen intersection point. C 7 delta(Z) is too large at the chosen intersection point. C 8 a track turning angle to the chosen vertex is too large C 19 no acceptable solution with an adequately positive arc length C Edit log: C ---- --- C 10/xx/94 WJA Allow NTRK = 1 or 2, for possible 1-track vertex; add C line-intersects-circle approximation for 1-track vertex; C add split-the-difference approximation for conversions; C detect negative JJ() values to trigger these special-case C first approximations. C=============================================================================== IMPLICIT NONE 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 INTEGER Nv, PRINT, IERR INTEGER I,J,K,L,M,N, JJ(2),II(2), NSOL, KPRINT INTEGER Nt, Mv, NTRK REAL HELIX(5,2), P(5),Q(5,5) REAL A,B, RA(2), RVI(2), TANG(2,2) REAL Xsvt,Ysvt, AA,BB,CC,DD,RR, Pt REAL TMP,PZV,FITA,FITB COMPLEX XYP0,PXYV,XYVR,PXYVR,XYTMP C Local variables for accumulating first-approximation momentum sums REAL SPh0,CPh0,SPhs,CPhs,SPhi,CPhi LOGICAL CONV DOUBLE PRECISION XC(2),YC(2),RC(2) DOUBLE PRECISION AB(2), XX,YY,Y0,YY2,COST,SINT,DX,DY,D,U,V INTEGER BEAM C Map the parameter order to CTVMFT form: INTEGER MAP(5) SAVE BEAM, MAP DATA BEAM /0/ C (Ctg,Crv,Z0,D0,Phi)->(Crv,Phi,Ctg,D0,Z0) DATA MAP / 3,1,5,4,2 / C Start of code C------------------------------------------------------------------------------- KPRINT = 0 IF (PRINT.GT.0 .AND. DBGPRT.GT.10) KPRINT=PRINT IERR = 0 IJKERR(1) = 0 C Find first two tracks (if they exist) JJ(1) = 0 JJ(2) = 0 DO Nt=1,NTRACK IF (TrkVtx(Nt,Nv)) THEN IF (JJ(1).EQ.0) THEN JJ(1) = Nt ELSE JJ(2) = Nt GO TO 1 ENDIF ENDIF ENDDO 1 CONTINUE CONV = .FALSE. NTRK = 2 C special case; conversion vertex IF (Cvtx(Nv).GT.0) CONV=.TRUE. C special case; "single track vertex" IF (JJ(2).EQ.0) THEN DO Mv=Nv+1,NVERTX IF (VtxPnt(Mv,1).EQ.Nv) THEN C For 1-track vertex, JJ(2) is a vertex NTRK = 1 C that points to this vertex JJ(2) = Mv END IF END DO END IF C loop over the two vertex-tracks DO I=1,NTRK C Track bank number J = IABS(LIST(JJ(I))) C dummy out MASS,RADIUS A = 0.0 CALL GETTRK (J,A,A,0, K,P,Q, IERR) C Quit, if error finding track data IF (IERR .NE. 0) THEN TKERR(JJ(I)) = IJKERR(2) GO TO 1000 ENDIF DO K=1,5 HELIX(MAP(K),I) = P(K) END DO END DO c----------------------------------------------------------------------- C For convenience, order the circles so that the first has smaller radius... c----------------------------------------------------------------------- II(1) = 1 II(2) = 2 IF (NTRK.EQ.2) THEN IF (ABS(HELIX(1,1)).LT.ABS(HELIX(1,2))) THEN II(1) = 2 II(2) = 1 END IF END IF C make the circle radii, centers DO I=1,NTRK C there are two tracks (circles) RC(I) = 0.5/HELIX(1,II(I)) U = RC(I)+HELIX(4,II(I)) XC(I) =-U*SIN(HELIX(2,II(I))) YC(I) = U*COS(HELIX(2,II(I))) RA(I) = ABS(RC(I)) END DO C there are usually two intersections DO J=1,2 XSVI(J) = 0.0 YSVI(J) = 0.0 RVI (J) = 2.0*RVMAX DZSVI(J)= 2.0*DZMAX DO I=1,2 SS(I,J) = 0.0 ZZ(I,J) = 0.0 TANG(I,J) = PI END DO END DO IF (NTRK.GT.1) GOTO 10 c----------------------------------------------------------------------- C Special handling of single-track vertex c----------------------------------------------------------------------- TMP = 1/(1+2*HELIX(1,1)*HELIX(4,1)) C Track parametrization FITA = HELIX(1,1)*TMP FITB = HELIX(4,1)*(1+HELIX(1,1)*HELIX(4,1))*TMP XYP0 = CMPLX(COS(HELIX(2,1)),SIN(HELIX(2,1))) C Daughter vertex momentum PXYV = CMPLX(VTXP4(1,JJ(2)),VTXP4(2,JJ(2))) TMP = 1/ABS(PXYV) PXYV = PXYV*TMP PZV = VTXP4(3,JJ(2))*TMP C Position & momentum, rotated to PHI0=0 coordinate system XYVR = CMPLX(XYZVRT(1,JJ(2)),XYZVRT(2,JJ(2)))*CONJG(XYP0) PXYVR = PXYV*CONJG(XYP0) C Solve quadratic equation for displacement DD from vertex to track AA = FITA BB = 0.5*IMAG(PXYVR)-FITA*REAL(XYVR*CONJG(PXYVR)) CC = FITA*ABS(XYVR)**2+FITB-IMAG(XYVR) TMP = BB*BB-AA*CC IF (TMP.GT.0) THEN NSOL = 2 IF (BB.LT.0) THEN AA = -AA BB = -BB CC = -CC ENDIF BB = BB+SQRT(TMP) DD = CC/BB DO I=1,2 XSVI(I) = XYZVRT(1,JJ(2))+DD*REAL(PXYV) YSVI(I) = XYZVRT(2,JJ(2))+DD*IMAG(PXYV) ZZ(2,I) = XYZVRT(3,JJ(2))+DD*PZV TANG(2,I) = 0 DD = BB/AA ENDDO GOTO 20 ENDIF c----------------------------------------------------------------------- C No intersection; find closest approach and use half the distance c----------------------------------------------------------------------- NSOL = 1 C Check for line and circle antiparallel at closest approach IF (REAL(PXYVR).LT.0) THEN PXYV = -PXYV PXYVR = -PXYVR PZV = -PZV ENDIF C Use SIN(PHI)**2/(1+COS(PHI)) for (1-COS(PHI)) to reduce roundoff CC = IMAG(PXYVR)*IMAG(PXYVR)/(1+REAL(PXYVR)) C Store 1/(2c) to avoid recalculation AA = 0.5/HELIX(1,1) C TMP is signed distance from circle to line, along (-SIN(PHI),COS(PHI)) TMP = (IMAG(XYVR)-HELIX(4,1))*REAL(PXYVR) & -REAL(XYVR)*IMAG(PXYVR)+CC*AA XYTMP = CMPLX(IMAG(PXYV),CC)*AA & +CMPLX(0.0,HELIX(4,1))+0.5*TMP*CMPLX(-IMAG(PXYV),REAL(PXYV)) DD = REAL((XYTMP-XYVR)*CONJG(PXYVR)) XYTMP = XYTMP*XYP0 XSVI(1) = REAL(XYTMP) YSVI(1) = IMAG(XYTMP) ZZ(2,1) = XYZVRT(3,JJ(2))+DD*PZV TANG(2,1) = 0 GOTO 20 C two track (or more) vertex 10 CONTINUE C get the circle center separation DX = XC(2) - XC(1) DY = YC(2) - YC(1) D = DX*DX+DY*DY C the circles are concentric IF (D.LE.0.0) THEN IJKERR(2) = 1 GO TO 1000 END IF D = DSQRT(D) C Separation (signed quantity); if >0, the U = D-RA(1)-RA(2) C two circles do not intersect, <0, they may FMCDIF(1) = U C Special handling of conversion vertex IF (CONV) THEN C the circles are too far appart to accept IF (ABS(U).GT.DRMAX) THEN IF (U.GT.0.0) THEN C ...exteriour IJKERR(2) = 2 ELSE C ...interiour IJKERR(2) = 3 END IF GO TO 1000 END IF NSOL = 1 C Vertex is track radius plus half of separation away from each center XSVI(1) = (XC(1)*(RA(2)+0.5*U)+XC(2)*(RA(1)+0.5*U))/D YSVI(1) = (YC(1)*(RA(2)+0.5*U)+YC(2)*(RA(1)+0.5*U))/D C Branch down to vertex acceptability checking GO TO 20 END IF c----------------------------------------------------------------------- C Rotate, translate to a system where the two circle centers lie c on the X axis. c----------------------------------------------------------------------- COST = DX/D SINT = DY/D ! Y` = Y1` = Y2` Y0 = (-XC(1)*YC(2)+XC(2)*YC(1))/D ! X1', X2' DO i=1,2 ab(i) = cost*xc(i)+sint*yc(i) ENDDO u = (xc(2)+xc(1))*(xc(2)-xc(1)) + (yc(2)+yc(1))*(yc(2)-yc(1)) v = (ra(2)+ra(1))*(ra(2)-ra(1)) ! the common circle X' (if they intersect) xx = 0.5 * (u-v)/d u = dsqrt((xx-ab(1))*(xx-ab(1))) ! y''**2 is positive if they intersect yy2 = (ra(1)+u)*(ra(1)-u) ! the circles intersect; IF (YY2.GT.0.0) THEN ! two intersection points (+/- Y) YY = DSQRT(YY2) ! invert the translation, rotation DO J=1,2 U = YY+Y0 XSVI(J) = COST*XX - SINT*U YSVI(J) = SINT*XX + COST*U YY =-YY END DO NSOL = 2 GO TO 20 END IF c----------------------------------------------------------------------- C We get here if the two circles do not intersect; C Find how close in the XY plane they approach each other, C and take the point half way between them as our vertex c approximation. c----------------------------------------------------------------------- U = D - (RA(1)+RA(2)) C A is outside of B IF (U.GT.0.0) THEN V = U J = 2 IF (AB(1).LT.AB(2)) J=1 XX = AB(J) + RA(J) C A is inside of B ELSE IF (AB(1).LT.AB(2)) THEN XX = AB(2) - RA(2) V = AB(1) - RA(1) - XX ELSE XX = AB(1) + RA(1) V = AB(2) + RA(2) - XX END IF END IF XX = XX + 0.5*V C rotate back to the original system XSVI(1) = COST*XX - SINT*Y0 YSVI(1) = SINT*XX + COST*Y0 NSOL = 1 C the circles are too far appart to accept IF (V.GT.DRMAX) THEN IF (U.GT.0.0) THEN C ...exteriour IJKERR(2) = 4 ELSE C ...interiour IJKERR(2) = 5 END IF GO TO 1000 END IF 20 CONTINUE c----------------------------------------------------------------------- C loop over solutions c----------------------------------------------------------------------- DO 30 J=1,NSOL ! loop over tracks DO I=1,NTRK ! point from circle center U = (XSVI(J)-XC(I))/RC(I) ! to the intersection vertex V =-(YSVI(J)-YC(I))/RC(I) ! turning angle from the track origin U = ATAN2(U,V) - HELIX(2,II(I)) IF (U.LT.-PI) THEN U = U + TWOPI ELSE IF (U.GT. PI) THEN U = U - TWOPI END IF TANG(I,J) = U C arc length SS(I,J) = RC(I)*U ZZ(I,J) = HELIX(5,II(I)) + SS(I,J)*HELIX(3,II(I)) END DO RVI(J) = SQRT(XSVI(J)*XSVI(J)+YSVI(J)*YSVI(J)) DZSVI(J) = ZZ(2,J) - ZZ(1,J) 30 CONTINUE C Check that there is at least one potentially acceptable solution! A = AMIN1(RVI(1),RVI(2)) C check the vertex radius is acceptable IF (A.GT.RVMAX) THEN IJKERR(2) = 6 GO TO 1000 END IF A = AMIN1(ABS(DZSVI(1)),ABS(DZSVI(2))) C check the Z difference is acceptable IF (A.GT.DZMAX) THEN IJKERR(2) = 7 GO TO 1000 END IF A = AMAX1(ABS(TANG(1,1)),ABS(TANG(2,1))) B = AMAX1(ABS(TANG(1,2)),ABS(TANG(2,2))) C check there is an acceptable TANG IF (AMIN1(A,B).GT.TRNMAX) THEN IJKERR(2) = 8 GO TO 1000 END IF C minimum track arc length, sol`n 1 A = AMIN1(SS(1,1),SS(2,1)) C minimum track arc length, sol`n 2 B = AMIN1(SS(1,2),SS(2,2)) C limit the minimum arc length IF (AMAX1(A,B).LT.DSMIN) THEN IJKERR(2) = 9 GO TO 1000 END IF C there may be a possible acceptable solution, in (R,deltaZ,S,Tang) J = 1 IF (NSOL.EQ.1) GO TO 40 IF (ABS(DZSVI(2)).LT.ABS(DZSVI(1))) J=2 FLIP = 0 A = AMAX1(ABS(TANG(1,J)),ABS(TANG(2,J))) B = AMIN1(SS(1,J),SS(2,J)) IF (RVI(J).GT.RVMAX) FLIP = 1 IF (ABS(DZSVI(J)).GT.DZMAX) FLIP = FLIP+2 IF (A.GT.TRNMAX) FLIP = FLIP+4 IF (B.LT.DSMIN) FLIP = FLIP+8 IF (FLIP.NE.0) THEN IF (J.EQ.1) THEN J = 2 ELSE J = 1 END IF END IF C final checks of the final solution... 40 CONTINUE A = AMAX1(ABS(TANG(1,J)),ABS(TANG(2,J))) B = AMIN1(SS(1,J),SS(2,J)) C check the vertex radius IF (RVI(J).GT.RVMAX) THEN IJKERR(2) = 6 GO TO 1000 END IF C check the Z difference IF (ABS(DZSVI(J)).GT.DZMAX) THEN IJKERR(2) = 7 GO TO 1000 END IF C check that TANG is acceptable IF (A.GT.TRNMAX) THEN IJKERR(2) = 8 GO TO 1000 END IF C limit the minimum arc length IF (B.LT.DSMIN) THEN IJKERR(2) = 9 GO TO 1000 END IF C "acceptable" solution. Set initial vertex approximation XYZVRT(1,Nv) = XSVI(J) XYZVRT(2,Nv) = YSVI(J) XYZVRT(3,Nv) = 0.5 * (ZZ(1,J) + ZZ(2,J)) c----------------------------------------------------------------------- C Go off and collect (GETTRK) the track data, and make the first c approximation C vertex momentum sums (which will be needed for single track vertex c first approximation calculations). C GETTRK failures or failures in transporting the tracks not used c in the vertex finding to the vertex will return IERR=IJKERR(1)=3. c----------------------------------------------------------------------- 50 CONTINUE D IF (DBGPRT.GT.0) THEN D RR = SQRT(XYZVRT(1,Nv)**2+XYZVRT(2,Nv)**2) D WRITE (PRINT,1010) Nv,(XYZVRT(I,Nv),I=1,3),RR D1010 FORMAT (/,' Vertex',I2, D & ' Approximation; ',3F10.4,5X,'RR',F7.2) D END IF C Initialize vertex 3-momentum to the sum VtxP4(1,Nv) = 0 C of daughter vertex 3-momenta VtxP4(2,Nv) = 0 C (needed for 1-track vertex CTVMFA) VtxP4(3,Nv) = 0 DO Mv=Nv+1,NVERTX IF (VtxPnt(Mv,1).EQ.Nv) THEN VtxP4(1,Nv) = VtxP4(1,Nv)+VtxP4(1,Mv) VtxP4(2,Nv) = VtxP4(2,Nv)+VtxP4(2,Mv) VtxP4(3,Nv) = VtxP4(3,Nv)+VtxP4(3,Mv) END IF END DO c----------------------------------------------------------------------- C Collect track information for this vertex c----------------------------------------------------------------------- DO Nt=1,NTRACK IF (TrkVtx(Nt,Nv)) THEN CALL CTVM01(KPRINT,Nt,BEAM,PARERR(1,Nt),IERR) IF (IERR.NE.0) THEN C GETTRK failure, or, track cannot reach vertex GO TO 1000 END IF TKERR(Nt) = 0 SPh0 = SIN(Par0(2,Nt)) CPh0 = COS(Par0(2,Nt)) Xsvt = XYZVRT(1,Nv)*CPh0+XYZVRT(2,Nv)*SPh0 Ysvt = XYZVRT(2,Nv)*CPh0-XYZVRT(1,Nv)*SPh0 C Proceed blissfully SPhs = 2*Par0(1,Nt)*Xsvt CPhs = SQRT((1-SPhs)*(1+SPhs)) SPhi = CPh0*Sphs+SPh0*CPhs CPhi = CPh0*CPhs-SPh0*SPhs C add this track's contribution to vertex momentum PT = PSCALE/ABS(Par0(1,Nt)) VtxP4(1,Nv) = VtxP4(1,Nv)+PT*CPhi VtxP4(2,Nv) = VtxP4(2,Nv)+PT*SPhi VtxP4(3,Nv) = VtxP4(3,Nv)+PT*Par0(3,Nt) END IF END DO 1000 CONTINUE IF (IERR.NE.0) THEN IJKERR(1) = IERR ELSE IF (IJKERR(2).NE.0) THEN IERR = 2 IJKERR(1) = 2 IJKERR(3) = Nv END IF RETURN END