--- DarthVader/TrackerLevel2/src/F77/grkuta.f 2006/05/19 13:15:55 1.1 +++ DarthVader/TrackerLevel2/src/F77/grkuta.f 2008/03/05 17:00:20 1.8 @@ -28,6 +28,7 @@ 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) @@ -37,17 +38,24 @@ + (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3)) * PARAMETER (MAXIT = 1992, MAXCUT = 11) - PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32) - PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3) +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) + 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 @@ -69,6 +77,9 @@ 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 * @@ -106,7 +117,10 @@ CALL GUFLD(VVV,FFF) DO I=1,3 F(I)=DBLE(FFF(I)) - ENDDO + 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) @@ -141,6 +155,9 @@ 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 @@ -237,7 +254,276 @@ * * - +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$$$* ********************************************************************** * @@ -252,15 +538,17 @@ 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 - call inter_B(vv(1),vv(2),vv(3),ff) - do i=1,3 !change back the field in kGauss - f(i)=ff(i)*10. - enddo - + 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)=ff(i)*10. + enddo +c$$$ print*,'GUFLD: b ',f + return end