********************************************************************** * * * routine per tracciare la particella di uno STEP * SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT) C. C. ****************************************************************** C. * * C. * Runge-Kutta method for tracking a particle through a magnetic * C. * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of * C. * Standards, procedure 25.5.20) * C. * * C. * Input parameters * C. * CHARGE Particle charge * C. * STEP Step size * C. * VECT Initial co-ords,direction cosines,momentum * C. * Output parameters * C. * VOUT Output co-ords,direction cosines,momentum * C. * User routine called * C. * CALL GUFLD(X,F) * C. * * C. * ==>Called by : , GUSWIM * C. * Authors R.Brun, M.Hansroul ********* * C. * V.Perevoztchikov (CUT STEP implementation) * C. * * C. * * C. ****************************************************************** C. IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/DELTAB/DELTA0,DELTA1,DLT * REAL VVV(3),FFF(3) REAL*8 CHARGE, STEP, VECT(*), VOUT(*), F(4) REAL*8 XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3) EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)), + (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3)) * PARAMETER (MAXIT = 1992, MAXCUT = 11) cPP PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32) PARAMETER (EC=2.99792458D-4) cPP PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3) PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0) PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO) PARAMETER (PISQUA=.986960440109D+01) PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6) REAL*8 DELTAB(3) REAL*8 DLT32 DLT32=DLT/32. *. *. ------------------------------------------------------------------ *. * This constant is for units CM,GEV/C and KGAUSS * ITER = 0 NCUT = 0 DO 10 J=1,7 VOUT(J)=VECT(J) 10 CONTINUE PINV = EC * CHARGE / VECT(7) TL = 0. H = STEP * * 20 REST = STEP-TL IF (DABS(H).GT.DABS(REST)) H = REST DO I=1,3 VVV(I)=SNGL(VOUT(I)) ENDDO CALL GUFLD(VVV,FFF) * print*,'GRKUTA Bx,By,Bz: ',(FFF(i),i=1,3) DO I=1,3 F(I)=DBLE(FFF(I)) ENDDO DELTAB(2) = -F(2)*VECT(7)*CHARGE*(DELTA0+DELTA1*VVV(2)) F(2) = F(2)+DELTAB(2) cPP ----------------- * * Start of integration * X = VOUT(1) Y = VOUT(2) Z = VOUT(3) A = VOUT(4) B = VOUT(5) C = VOUT(6) * H2 = HALF * H H4 = HALF * H2 PH = PINV * H PH2 = HALF * PH SECXS(1) = (B * F(3) - C * F(2)) * PH2 SECYS(1) = (C * F(1) - A * F(3)) * PH2 SECZS(1) = (A * F(2) - B * F(1)) * PH2 ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2) IF (ANG2.GT.PISQUA) GO TO 40 DXT = H2 * A + H4 * SECXS(1) DYT = H2 * B + H4 * SECYS(1) DZT = H2 * C + H4 * SECZS(1) XT = X + DXT YT = Y + DYT ZT = Z + DZT * * Second intermediate point * EST = DABS(DXT)+DABS(DYT)+DABS(DZT) IF (EST.GT.H) GO TO 30 DO I=1,3 VVV(I)=SNGL(XYZT(I)) ENDDO CALL GUFLD(VVV,FFF) DO I=1,3 F(I)=DBLE(FFF(I)) ENDDO DELTAB(2) = -F(2)*VECT(7)*CHARGE*(DELTA0+DELTA1*VVV(2)) F(2) = F(2)+DELTAB(2) cPP ----------------- C CALL GUFLD(XYZT,F) AT = A + SECXS(1) BT = B + SECYS(1) CT = C + SECZS(1) * SECXS(2) = (BT * F(3) - CT * F(2)) * PH2 SECYS(2) = (CT * F(1) - AT * F(3)) * PH2 SECZS(2) = (AT * F(2) - BT * F(1)) * PH2 AT = A + SECXS(2) BT = B + SECYS(2) CT = C + SECZS(2) SECXS(3) = (BT * F(3) - CT * F(2)) * PH2 SECYS(3) = (CT * F(1) - AT * F(3)) * PH2 SECZS(3) = (AT * F(2) - BT * F(1)) * PH2 DXT = H * (A + SECXS(3)) DYT = H * (B + SECYS(3)) DZT = H * (C + SECZS(3)) XT = X + DXT YT = Y + DYT ZT = Z + DZT AT = A + TWO*SECXS(3) BT = B + TWO*SECYS(3) CT = C + TWO*SECZS(3) * EST = ABS(DXT)+ABS(DYT)+ABS(DZT) IF (EST.GT.2.*ABS(H)) GO TO 30 DO I=1,3 VVV(I)=SNGL(XYZT(I)) ENDDO CALL GUFLD(VVV,FFF) DO I=1,3 F(I)=DBLE(FFF(I)) ENDDO DELTAB(2) = -F(2)*VECT(7)*CHARGE*(DELTA0+DELTA1*VVV(2)) F(2) = F(2)+DELTAB(2) cPP ----------------- C CALL GUFLD(XYZT,F) * Z = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H Y = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H X = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H * SECXS(4) = (BT*F(3) - CT*F(2))* PH2 SECYS(4) = (CT*F(1) - AT*F(3))* PH2 SECZS(4) = (AT*F(2) - BT*F(1))* PH2 A = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD B = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD C = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD * EST = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3))) ++ ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3))) ++ ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3))) * IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30 ITER = ITER + 1 NCUT = 0 * If too many iterations, go to HELIX IF (ITER.GT.MAXIT) GO TO 40 * TL = TL + H IF (EST.LT.(DLT32)) THEN H = H*TWO ENDIF CBA = ONE/ SQRT(A*A + B*B + C*C) VOUT(1) = X VOUT(2) = Y VOUT(3) = Z VOUT(4) = CBA*A VOUT(5) = CBA*B VOUT(6) = CBA*C REST = STEP - TL IF (STEP.LT.0.) REST = -REST IF (REST .GT. 1.E-5*DABS(STEP)) GO TO 20 * GO TO 999 * ** CUT STEP 30 NCUT = NCUT + 1 * If too many cuts , go to HELIX IF (NCUT.GT.MAXCUT) GO TO 40 H = H*HALF GO TO 20 * ** ANGLE TOO BIG, USE HELIX 40 F1 = F(1) F2 = F(2) F3 = F(3) F4 = DSQRT(F1**2+F2**2+F3**2) RHO = -F4*PINV TET = RHO * STEP IF(TET.NE.0.) THEN HNORM = ONE/F4 F1 = F1*HNORM F2 = F2*HNORM F3 = F3*HNORM * HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY) HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ) HXP(3) = F1*VECT(IPY) - F2*VECT(IPX) HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ) * RHO1 = ONE/RHO SINT = DSIN(TET) COST = TWO*DSIN(HALF*TET)**2 * G1 = SINT*RHO1 G2 = COST*RHO1 G3 = (TET-SINT) * HP*RHO1 G4 = -COST G5 = SINT G6 = COST * HP VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1) VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2) VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3) VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1) VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2) VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3) * ELSE VOUT(IX) = VECT(IX) + STEP*VECT(IPX) VOUT(IY) = VECT(IY) + STEP*VECT(IPY) VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ) * ENDIF * 999 END * * c$$$********************************************************************** c$$$* c$$$* c$$$* routine per tracciare la particella di uno STEP c$$$* *** extended version *** c$$$* it return also the track-length c$$$* c$$$ SUBROUTINE GRKUTA2 (CHARGE,STEP,VECT,VOUT) c$$$C. c$$$C. ****************************************************************** c$$$C. * * c$$$C. * Runge-Kutta method for tracking a particle through a magnetic * c$$$C. * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of * c$$$C. * Standards, procedure 25.5.20) * c$$$C. * * c$$$C. * Input parameters * c$$$C. * CHARGE Particle charge * c$$$C. * STEP Step size * c$$$C. * VECT Initial co-ords,direction cosines,momentum * c$$$C. * Output parameters * c$$$C. * VOUT Output co-ords,direction cosines,momentum * c$$$C. * User routine called * c$$$C. * CALL GUFLD(X,F) * c$$$C. * * c$$$C. * ==>Called by : , GUSWIM * c$$$C. * Authors R.Brun, M.Hansroul ********* * c$$$C. * V.Perevoztchikov (CUT STEP implementation) * c$$$C. * * c$$$C. * * c$$$C. ****************************************************************** c$$$C. c$$$ IMPLICIT DOUBLE PRECISION(A-H,O-Z) c$$$* c$$$ REAL VVV(3),FFF(3) c$$$ REAL*8 CHARGE, STEP, VECT(*), VOUT(*), F(4) c$$$ REAL*8 XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT c$$$ DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3) c$$$ EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)), c$$$ + (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3)) c$$$* c$$$ PARAMETER (MAXIT = 1992, MAXCUT = 11) c$$$ PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32) c$$$ PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3) c$$$ PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO) c$$$ PARAMETER (PISQUA=.986960440109D+01) c$$$ PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6) c$$$ c$$$* track length c$$$ REAL*8 DL c$$$ c$$$*. c$$$*. ------------------------------------------------------------------ c$$$*. c$$$* This constant is for units CM,GEV/C and KGAUSS c$$$* c$$$ ITER = 0 c$$$ NCUT = 0 c$$$ DO 10 J=1,8 c$$$ VOUT(J)=VECT(J) c$$$ 10 CONTINUE c$$$ PINV = EC * CHARGE / VECT(7) c$$$ TL = 0. c$$$ H = STEP c$$$ c$$$c print*,'===================== START GRKUTA2' c$$$ c$$$* c$$$* c$$$ 20 REST = STEP-TL c$$$ IF (DABS(H).GT.DABS(REST)) H = REST c$$$ DO I=1,3 c$$$ VVV(I)=SNGL(VOUT(I)) c$$$ ENDDO c$$$ c$$$ CALL GUFLD(VVV,FFF) c$$$* print*,'GRKUTA Bx,By,Bz: ',(FFF(i),i=1,3) c$$$ DO I=1,3 c$$$ F(I)=DBLE(FFF(I)) c$$$ ENDDO c$$$* c$$$* Start of integration c$$$* c$$$ X = VOUT(1) c$$$ Y = VOUT(2) c$$$ Z = VOUT(3) c$$$ A = VOUT(4) c$$$ B = VOUT(5) c$$$ C = VOUT(6) c$$$ c$$$ DL = VOUT(8) c$$$ c$$$* c$$$ H2 = HALF * H c$$$ H4 = HALF * H2 c$$$ PH = PINV * H c$$$ PH2 = HALF * PH c$$$ SECXS(1) = (B * F(3) - C * F(2)) * PH2 c$$$ SECYS(1) = (C * F(1) - A * F(3)) * PH2 c$$$ SECZS(1) = (A * F(2) - B * F(1)) * PH2 c$$$ ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2) c$$$ IF (ANG2.GT.PISQUA) GO TO 40 c$$$ DXT = H2 * A + H4 * SECXS(1) c$$$ DYT = H2 * B + H4 * SECYS(1) c$$$ DZT = H2 * C + H4 * SECZS(1) c$$$ XT = X + DXT c$$$ YT = Y + DYT c$$$ ZT = Z + DZT c$$$* c$$$* Second intermediate point c$$$* c$$$ EST = DABS(DXT)+DABS(DYT)+DABS(DZT) c$$$ IF (EST.GT.H) GO TO 30 c$$$ c$$$ DO I=1,3 c$$$ VVV(I)=SNGL(XYZT(I)) c$$$ ENDDO c$$$ CALL GUFLD(VVV,FFF) c$$$ DO I=1,3 c$$$ F(I)=DBLE(FFF(I)) c$$$ ENDDO c$$$C CALL GUFLD(XYZT,F) c$$$ AT = A + SECXS(1) c$$$ BT = B + SECYS(1) c$$$ CT = C + SECZS(1) c$$$* c$$$ SECXS(2) = (BT * F(3) - CT * F(2)) * PH2 c$$$ SECYS(2) = (CT * F(1) - AT * F(3)) * PH2 c$$$ SECZS(2) = (AT * F(2) - BT * F(1)) * PH2 c$$$ AT = A + SECXS(2) c$$$ BT = B + SECYS(2) c$$$ CT = C + SECZS(2) c$$$ SECXS(3) = (BT * F(3) - CT * F(2)) * PH2 c$$$ SECYS(3) = (CT * F(1) - AT * F(3)) * PH2 c$$$ SECZS(3) = (AT * F(2) - BT * F(1)) * PH2 c$$$ DXT = H * (A + SECXS(3)) c$$$ DYT = H * (B + SECYS(3)) c$$$ DZT = H * (C + SECZS(3)) c$$$ XT = X + DXT c$$$ YT = Y + DYT c$$$ ZT = Z + DZT c$$$ AT = A + TWO*SECXS(3) c$$$ BT = B + TWO*SECYS(3) c$$$ CT = C + TWO*SECZS(3) c$$$* c$$$ EST = ABS(DXT)+ABS(DYT)+ABS(DZT) c$$$ IF (EST.GT.2.*ABS(H)) GO TO 30 c$$$ c$$$ DO I=1,3 c$$$ VVV(I)=SNGL(XYZT(I)) c$$$ ENDDO c$$$ CALL GUFLD(VVV,FFF) c$$$ DO I=1,3 c$$$ F(I)=DBLE(FFF(I)) c$$$ ENDDO c$$$C CALL GUFLD(XYZT,F) c$$$* c$$$ Z = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H c$$$ Y = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H c$$$ X = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H c$$$* c$$$ SECXS(4) = (BT*F(3) - CT*F(2))* PH2 c$$$ SECYS(4) = (CT*F(1) - AT*F(3))* PH2 c$$$ SECZS(4) = (AT*F(2) - BT*F(1))* PH2 c$$$ A = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD c$$$ B = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD c$$$ C = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD c$$$* c$$$ EST = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3))) c$$$ ++ ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3))) c$$$ ++ ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3))) c$$$* c$$$ IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30 c$$$ c$$$ ITER = ITER + 1 c$$$ NCUT = 0 c$$$* If too many iterations, go to HELIX c$$$ IF (ITER.GT.MAXIT) GO TO 40 c$$$* c$$$ DL = VOUT(8) + c$$$ $ DSQRT( 0 c$$$ $ + (X-VOUT(1))**2 c$$$ $ + (Y-VOUT(2))**2 c$$$ $ + (Z-VOUT(3))**2 c$$$ $ ) c$$$c print*,'- ',VOUT(3),z,VOUT(1),x,VOUT(2),y,DL c$$$* c$$$ TL = TL + H c$$$ IF (EST.LT.(DLT32)) THEN c$$$ H = H*TWO c$$$ ENDIF c$$$ CBA = ONE/ SQRT(A*A + B*B + C*C) c$$$ VOUT(1) = X c$$$ VOUT(2) = Y c$$$ VOUT(3) = Z c$$$ VOUT(4) = CBA*A c$$$ VOUT(5) = CBA*B c$$$ VOUT(6) = CBA*C c$$$ VOUT(8) = DL c$$$ REST = STEP - TL c$$$ IF (STEP.LT.0.) REST = -REST c$$$ IF (REST .GT. 1.E-5*DABS(STEP)) GO TO 20 c$$$* c$$$ GO TO 999 c$$$* c$$$** CUT STEP c$$$ 30 NCUT = NCUT + 1 c$$$* If too many cuts , go to HELIX c$$$ IF (NCUT.GT.MAXCUT) GO TO 40 c$$$ H = H*HALF c$$$ GO TO 20 c$$$* c$$$** ANGLE TOO BIG, USE HELIX c$$$ 40 F1 = F(1) c$$$ F2 = F(2) c$$$ F3 = F(3) c$$$ F4 = DSQRT(F1**2+F2**2+F3**2) c$$$ RHO = -F4*PINV c$$$ TET = RHO * STEP c$$$ IF(TET.NE.0.) THEN c$$$ HNORM = ONE/F4 c$$$ F1 = F1*HNORM c$$$ F2 = F2*HNORM c$$$ F3 = F3*HNORM c$$$* c$$$ HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY) c$$$ HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ) c$$$ HXP(3) = F1*VECT(IPY) - F2*VECT(IPX) c$$$ c$$$ HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ) c$$$* c$$$ RHO1 = ONE/RHO c$$$ SINT = DSIN(TET) c$$$ COST = TWO*DSIN(HALF*TET)**2 c$$$* c$$$ G1 = SINT*RHO1 c$$$ G2 = COST*RHO1 c$$$ G3 = (TET-SINT) * HP*RHO1 c$$$ G4 = -COST c$$$ G5 = SINT c$$$ G6 = COST * HP c$$$ c$$$ VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1) c$$$ VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2) c$$$ VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3) c$$$ c$$$ VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1) c$$$ VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2) c$$$ VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3) c$$$* c$$$ ELSE c$$$ VOUT(IX) = VECT(IX) + STEP*VECT(IPX) c$$$ VOUT(IY) = VECT(IY) + STEP*VECT(IPY) c$$$ VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ) c$$$* c$$$ ENDIF c$$$* TEMP !!! TEMP !!! TEMP !!! TEMP !!! TEMP !!! TEMP !!! c$$$* devo mettere la lunghezza dell'elica!!!!!!!!!!!!!! c$$$* ma non mi riesce :-( c$$$ VOUT(8) = DSQRT( 0 c$$$ $ +(VOUT(IX)-VECT(IX))**2 c$$$ $ +(VOUT(IY)-VECT(IY))**2 c$$$ $ +(VOUT(IZ)-VECT(IZ))**2 c$$$ $ ) c$$$c print*,'WARNING: GRKUTA2 --> ' c$$$c $ ,'helix :-( ... length evaluated with straight line' c$$$ c$$$* c$$$ 999 END c$$$* c$$$* ********************************************************************** * * gives the value of the magnetic field in the tracking point * ********************************************************************** subroutine gufld(v,f) !coordinates in cm, B field in kGauss real v(3),f(3) !coordinates in cm, B field in kGauss, error in kGauss real*8 vv(3),ff(3) !inter_B.f works in double precision do i=1,3 vv(i)=v(i)/100. !inter_B.f works in meters enddo c inter_B: coordinates in m, B field in Tesla c$$$ print*,'GUFLD: v ',v call inter_B(vv(1),vv(2),vv(3),ff) do i=1,3 !change back the field in kGauss f(i)=REAL(ff(i)*10.) ! EM GCC4.7 enddo c$$$ print*,'GUFLD: b ',f return end