MUON/extrap.F

   1 * One step extrapolation routines in Fortran:
   2 * EXTRAP_ONESTEP_HELIX, EXTRAP_ONESTEP_HELIX3, EXTRAP_ONESTEP_RUNGEKUTTA.
   3 * Taken respectively from Geant (gtrak) routines:
   4 * GHELIX, GHELIX3, GRKUTA.
   5 * Everything in double precision,
   6 * in order that the track fit with Minuit is converging.
   7 * Modifications from Geant are indicated with "Cmodif".
   8
   9 Cmodif: SUBROUTINE GHELIX (CHARGE, STEP, VECT, VOUT) changed into:
  10       SUBROUTINE EXTRAP_ONESTEP_HELIX (CHARGE, STEP, VECT, VOUT)
  11 C.
  12 C.    ******************************************************************
  13 C.    *                                                                *
  14 C.    *  Performs the tracking of one step in a magnetic field         *
  15 C.    *  The trajectory is assumed to be a helix in a constant field   *
  16 C.    *  taken at the mid point of the step.                           *
  17 C.    *  Parameters:                                                   *
  18 C.    *   input                                                        *
  19 C.    *     STEP =arc length of the step asked                         *
  20 C.    *     VECT =input vector (position,direction cos and momentum)   *
  21 C.    *     CHARGE=  electric charge of the particle                   *
  22 C.    *   output                                                       *
  23 C.    *     VOUT = same as VECT after completion of the step           *
  24 C.    *                                                                *
  25 C.    *    ==>Called by : <USER>, GUSWIM                               *
  26 C.    *       Author    M.Hansroul  *********                          *
  27 C.    *       Modified  S.Egli, S.V.Levonian                           *
  28 C.    *       Modified  V.Perevoztchikov
  29 C.    *                                                                *
  30 C.    ******************************************************************
  31 C.
  32
  33 Cmodif: everything in double precision
  34       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  35
  36       DIMENSION      VECT(7),VOUT(7)
  37       DIMENSION      XYZ(3),H(4),HXP(3)
  38       PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
  39       PARAMETER      (SIXTH = 1./6.)
  40       PARAMETER      (EC=2.9979251E-4)
  41 C.
  42 C.    ------------------------------------------------------------------
  43 C.
  44 C       units are kgauss,centimeters,gev/c
  45 C
  46       VOUT(IPP) = VECT(IPP)
  47       IF (CHARGE.EQ.0.)         GO TO 10
  48       XYZ(1)    = VECT(IX) + 0.5 * STEP * VECT(IPX)
  49       XYZ(2)    = VECT(IY) + 0.5 * STEP * VECT(IPY)
  50       XYZ(3)    = VECT(IZ) + 0.5 * STEP * VECT(IPZ)
  51 C
  52 Cmodif: CALL GUFLD (XYZ, H) changed into:
  53       CALL GUFLD_DOUBLE (XYZ, H)
  54
  55       H2XY = H(1)**2 + H(2)**2
  56       H(4) = H(3)**2 + H2XY
  57       IF (H(4).LE.1.E-12)       GO TO 10
  58       IF (H2XY.LE.1.E-12*H(4))  THEN
  59 Cmodif: CALL GHELX3 (CHARGE*H(3), STEP, VECT, VOUT) changed into:
  60          CALL EXTRAP_ONESTEP_HELIX3 (CHARGE*H(3), STEP, VECT, VOUT)
  61          GO TO 999
  62       ENDIF
  63       H(4) = SQRT(H(4))
  64       H(1) = H(1)/H(4)
  65       H(2) = H(2)/H(4)
  66       H(3) = H(3)/H(4)
  67       H(4) = H(4)*EC
  68 *
  69       HXP(1) = H(2)*VECT(IPZ) - H(3)*VECT(IPY)
  70       HXP(2) = H(3)*VECT(IPX) - H(1)*VECT(IPZ)
  71       HXP(3) = H(1)*VECT(IPY) - H(2)*VECT(IPX)
  72
  73       HP = H(1)*VECT(IPX) + H(2)*VECT(IPY) + H(3)*VECT(IPZ)
  74 *
  75       RHO = -CHARGE*H(4)/VECT(IPP)
  76       TET = RHO * STEP
  77       IF (ABS(TET).GT.0.15)     THEN
  78          SINT = SIN(TET)
  79          SINTT = (SINT/TET)
  80          TSINT = (TET-SINT)/TET
  81          COS1T = 2.*(SIN(0.5*TET))**2/TET
  82       ELSE
  83          TSINT = SIXTH*TET**2
  84          SINTT = (1. - TSINT)
  85          SINT = TET*SINTT
  86          COS1T = 0.5*TET
  87       ENDIF
  88 *
  89       F1 = STEP * SINTT
  90       F2 = STEP * COS1T
  91       F3 = STEP * TSINT * HP
  92       F4 = -TET*COS1T
  93       F5 = SINT
  94       F6 = TET * COS1T * HP
  95
  96       VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1) + F3*H(1))
  97       VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2) + F3*H(2))
  98       VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F2*HXP(3) + F3*H(3))
  99
 100       VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1) + F6*H(1))
 101       VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2) + F6*H(2))
 102       VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F5*HXP(3) + F6*H(3))
 103
 104       GO TO 999
 105
 106    10 CONTINUE
 107       DO 20 I   = 1,3
 108          VOUT(I) = VECT(I) + STEP * VECT(I+3)
 109          VOUT(I+3) = VECT(I+3)
 110    20 CONTINUE
 111 C
 112   999 END
 113
 114 Cmodif: SUBROUTINE GHELX3 (FIELD, STEP, VECT, VOUT) changed into:
 115       SUBROUTINE EXTRAP_ONESTEP_HELIX3 (FIELD, STEP, VECT, VOUT)
 116 C.
 117 C.    ******************************************************************
 118 C.    *                                                                *
 119 C.    *       Tracking routine in a constant field oriented            *
 120 C.    *       along axis 3                                             *
 121 C.    *       Tracking is performed with a conventional                *
 122 C.    *       helix step method                                        *
 123 C.    *                                                                *
 124 C.    *    ==>Called by : <USER>, GUSWIM                               *
 125 C.    *       Authors    R.Brun, M.Hansroul  *********                 *
 126 C     *       Rewritten  V.Perevoztchikov
 127 C.    *                                                                *
 128 C.    ******************************************************************
 129 C.
 130
 131 Cmodif: everything in double precision
 132       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 133
 134       DIMENSION      VECT(7),VOUT(7),HXP(3)
 135       PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6,IPP=7)
 136       PARAMETER      (SIXTH = 1./6.)
 137       PARAMETER      (EC=2.9979251E-4)
 138 C.
 139 C.    ------------------------------------------------------------------
 140 C.
 141 C       units are kgauss,centimeters,gev/c
 142 C
 143       VOUT(IPP) = VECT(IPP)
 144       H4 = FIELD * EC
 145 *
 146       HXP(1) = - VECT(IPY)
 147       HXP(2) = + VECT(IPX)
 148
 149       HP = VECT(IPZ)
 150 *
 151       RHO = -H4/VECT(IPP)
 152       TET = RHO * STEP
 153       IF (ABS(TET).GT.0.15)     THEN
 154          SINT = SIN(TET)
 155          SINTT = (SINT/TET)
 156          TSINT = (TET-SINT)/TET
 157          COS1T = 2.*(SIN(0.5*TET))**2/TET
 158       ELSE
 159          TSINT = SIXTH*TET**2
 160          SINTT = (1. - TSINT)
 161          SINT = TET*SINTT
 162          COS1T = 0.5*TET
 163       ENDIF
 164 *
 165       F1 = STEP * SINTT
 166       F2 = STEP * COS1T
 167       F3 = STEP * TSINT * HP
 168       F4 = -TET*COS1T
 169       F5 = SINT
 170       F6 = TET * COS1T * HP
 171
 172       VOUT(IX) = VECT(IX) + (F1*VECT(IPX) + F2*HXP(1))
 173       VOUT(IY) = VECT(IY) + (F1*VECT(IPY) + F2*HXP(2))
 174       VOUT(IZ) = VECT(IZ) + (F1*VECT(IPZ) + F3)
 175
 176       VOUT(IPX) = VECT(IPX) + (F4*VECT(IPX) + F5*HXP(1))
 177       VOUT(IPY) = VECT(IPY) + (F4*VECT(IPY) + F5*HXP(2))
 178       VOUT(IPZ) = VECT(IPZ) + (F4*VECT(IPZ) + F6)
 179
 180 C
 181   999 END
 182
 183 Cmodif: SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT) changed into:
 184       SUBROUTINE EXTRAP_ONESTEP_RUNGEKUTTA (CHARGE,STEP,VECT,VOUT)
 185 C.
 186 C.    ******************************************************************
 187 C.    *                                                                *
 188 C.    *  Runge-Kutta method for tracking a particle through a magnetic *
 189 C.    *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
 190 C.    *  Standards, procedure 25.5.20)                                 *
 191 C.    *                                                                *
 192 C.    *  Input parameters                                              *
 193 C.    *       CHARGE    Particle charge                                *
 194 C.    *       STEP      Step size                                      *
 195 C.    *       VECT      Initial co-ords,direction cosines,momentum     *
 196 C.    *  Output parameters                                             *
 197 C.    *       VOUT      Output co-ords,direction cosines,momentum      *
 198 C.    *  User routine called                                           *
 199 C.    *       CALL GUFLD(X,F)                                          *
 200 C.    *                                                                *
 201 C.    *    ==>Called by : <USER>, GUSWIM                               *
 202 C.    *       Authors    R.Brun, M.Hansroul  *********                 *
 203 C.    *                  V.Perevoztchikov (CUT STEP implementation)    *
 204 C.    *                                                                *
 205 C.    *                                                                *
 206 C.    ******************************************************************
 207 C.
 208 Cmodif: no condition from CERNLIB_SINGLE for double precision
 209       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 210 Cmodif: REAL changed into DOUBLE PRECISION in following 2 lines
 211       DOUBLE PRECISION CHARGE, STEP, VECT(*), VOUT(*), F(4)
 212       DOUBLE PRECISION XYZT(3), XYZ(3), X, Y, Z, XT, YT, ZT
 213       DIMENSION SECXS(4),SECYS(4),SECZS(4),HXP(3)
 214       EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
 215      +            (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
 216 *
 217       PARAMETER (MAXIT = 1992, MAXCUT = 11)
 218       PARAMETER (EC=2.9979251D-4,DLT=1D-4,DLT32=DLT/32)
 219       PARAMETER (ZERO=0, ONE=1, TWO=2, THREE=3)
 220       PARAMETER (THIRD=ONE/THREE, HALF=ONE/TWO)
 221       PARAMETER (PISQUA=.986960440109D+01)
 222       PARAMETER      (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
 223 *.
 224 *.    ------------------------------------------------------------------
 225 *.
 226 *             This constant is for units CM,GEV/C and KGAUSS
 227 *
 228       ITER = 0
 229       NCUT = 0
 230       DO 10 J=1,7
 231          VOUT(J)=VECT(J)
 232    10 CONTINUE
 233       PINV   = EC * CHARGE / VECT(7)
 234       TL = 0.
 235       H      = STEP
 236 *
 237 *
 238    20 REST  = STEP-TL
 239       IF (ABS(H).GT.ABS(REST)) H = REST
 240 Cmodif: CALL GUFLD(VOUT,F) changed into:
 241       CALL GUFLD_DOUBLE(VOUT,F)
 242 *
 243 *             Start of integration
 244 *
 245       X      = VOUT(1)
 246       Y      = VOUT(2)
 247       Z      = VOUT(3)
 248       A      = VOUT(4)
 249       B      = VOUT(5)
 250       C      = VOUT(6)
 251 *
 252       H2     = HALF * H
 253       H4     = HALF * H2
 254       PH     = PINV * H
 255       PH2    = HALF * PH
 256       SECXS(1) = (B * F(3) - C * F(2)) * PH2
 257       SECYS(1) = (C * F(1) - A * F(3)) * PH2
 258       SECZS(1) = (A * F(2) - B * F(1)) * PH2
 259       ANG2 = (SECXS(1)**2 + SECYS(1)**2 + SECZS(1)**2)
 260       IF (ANG2.GT.PISQUA) GO TO 40
 261       DXT    = H2 * A + H4 * SECXS(1)
 262       DYT    = H2 * B + H4 * SECYS(1)
 263       DZT    = H2 * C + H4 * SECZS(1)
 264       XT     = X + DXT
 265       YT     = Y + DYT
 266       ZT     = Z + DZT
 267 *
 268 *              Second intermediate point
 269 *
 270       EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
 271       IF (EST.GT.H) GO TO 30
 272
 273 Cmodif: CALL GUFLD(XYZT,F) changed into:
 274       CALL GUFLD_DOUBLE(XYZT,F)
 275       AT     = A + SECXS(1)
 276       BT     = B + SECYS(1)
 277       CT     = C + SECZS(1)
 278 *
 279       SECXS(2) = (BT * F(3) - CT * F(2)) * PH2
 280       SECYS(2) = (CT * F(1) - AT * F(3)) * PH2
 281       SECZS(2) = (AT * F(2) - BT * F(1)) * PH2
 282       AT     = A + SECXS(2)
 283       BT     = B + SECYS(2)
 284       CT     = C + SECZS(2)
 285       SECXS(3) = (BT * F(3) - CT * F(2)) * PH2
 286       SECYS(3) = (CT * F(1) - AT * F(3)) * PH2
 287       SECZS(3) = (AT * F(2) - BT * F(1)) * PH2
 288       DXT    = H * (A + SECXS(3))
 289       DYT    = H * (B + SECYS(3))
 290       DZT    = H * (C + SECZS(3))
 291       XT     = X + DXT
 292       YT     = Y + DYT
 293       ZT     = Z + DZT
 294       AT     = A + TWO*SECXS(3)
 295       BT     = B + TWO*SECYS(3)
 296       CT     = C + TWO*SECZS(3)
 297 *
 298       EST = ABS(DXT)+ABS(DYT)+ABS(DZT)
 299       IF (EST.GT.2.*ABS(H)) GO TO 30
 300
 301 Cmodif: CALL GUFLD(XYZT,F) changed into:
 302       CALL GUFLD_DOUBLE(XYZT,F)
 303 *
 304       Z      = Z + (C + (SECZS(1) + SECZS(2) + SECZS(3)) * THIRD) * H
 305       Y      = Y + (B + (SECYS(1) + SECYS(2) + SECYS(3)) * THIRD) * H
 306       X      = X + (A + (SECXS(1) + SECXS(2) + SECXS(3)) * THIRD) * H
 307 *
 308       SECXS(4) = (BT*F(3) - CT*F(2))* PH2
 309       SECYS(4) = (CT*F(1) - AT*F(3))* PH2
 310       SECZS(4) = (AT*F(2) - BT*F(1))* PH2
 311       A      = A+(SECXS(1)+SECXS(4)+TWO * (SECXS(2)+SECXS(3))) * THIRD
 312       B      = B+(SECYS(1)+SECYS(4)+TWO * (SECYS(2)+SECYS(3))) * THIRD
 313       C      = C+(SECZS(1)+SECZS(4)+TWO * (SECZS(2)+SECZS(3))) * THIRD
 314 *
 315       EST    = ABS(SECXS(1)+SECXS(4) - (SECXS(2)+SECXS(3)))
 316      ++        ABS(SECYS(1)+SECYS(4) - (SECYS(2)+SECYS(3)))
 317      ++        ABS(SECZS(1)+SECZS(4) - (SECZS(2)+SECZS(3)))
 318 *
 319       IF (EST.GT.DLT .AND. ABS(H).GT.1.E-4) GO TO 30
 320       ITER = ITER + 1
 321       NCUT = 0
 322 *               If too many iterations, go to HELIX
 323       IF (ITER.GT.MAXIT) GO TO 40
 324 *
 325       TL = TL + H
 326       IF (EST.LT.(DLT32)) THEN
 327          H = H*TWO
 328       ENDIF
 329       CBA    = ONE/ SQRT(A*A + B*B + C*C)
 330       VOUT(1) = X
 331       VOUT(2) = Y
 332       VOUT(3) = Z
 333       VOUT(4) = CBA*A
 334       VOUT(5) = CBA*B
 335       VOUT(6) = CBA*C
 336       REST = STEP - TL
 337       IF (STEP.LT.0.) REST = -REST
 338       IF (REST .GT. 1.E-5*ABS(STEP)) GO TO 20
 339 *
 340       GO TO 999
 341 *
 342 **              CUT STEP
 343    30 NCUT = NCUT + 1
 344 *               If too many cuts , go to HELIX
 345       IF (NCUT.GT.MAXCUT)       GO TO 40
 346       H = H*HALF
 347       GO TO 20
 348 *
 349 **              ANGLE TOO BIG, USE HELIX
 350    40 F1  = F(1)
 351       F2  = F(2)
 352       F3  = F(3)
 353       F4  = SQRT(F1**2+F2**2+F3**2)
 354       RHO = -F4*PINV
 355       TET = RHO * STEP
 356       IF(TET.NE.0.) THEN
 357          HNORM = ONE/F4
 358          F1 = F1*HNORM
 359          F2 = F2*HNORM
 360          F3 = F3*HNORM
 361 *
 362          HXP(1) = F2*VECT(IPZ) - F3*VECT(IPY)
 363          HXP(2) = F3*VECT(IPX) - F1*VECT(IPZ)
 364          HXP(3) = F1*VECT(IPY) - F2*VECT(IPX)
 365
 366          HP = F1*VECT(IPX) + F2*VECT(IPY) + F3*VECT(IPZ)
 367 *
 368          RHO1 = ONE/RHO
 369          SINT = SIN(TET)
 370          COST = TWO*SIN(HALF*TET)**2
 371 *
 372          G1 = SINT*RHO1
 373          G2 = COST*RHO1
 374          G3 = (TET-SINT) * HP*RHO1
 375          G4 = -COST
 376          G5 = SINT
 377          G6 = COST * HP
 378
 379          VOUT(IX) = VECT(IX) + (G1*VECT(IPX) + G2*HXP(1) + G3*F1)
 380          VOUT(IY) = VECT(IY) + (G1*VECT(IPY) + G2*HXP(2) + G3*F2)
 381          VOUT(IZ) = VECT(IZ) + (G1*VECT(IPZ) + G2*HXP(3) + G3*F3)
 382
 383          VOUT(IPX) = VECT(IPX) + (G4*VECT(IPX) + G5*HXP(1) + G6*F1)
 384          VOUT(IPY) = VECT(IPY) + (G4*VECT(IPY) + G5*HXP(2) + G6*F2)
 385          VOUT(IPZ) = VECT(IPZ) + (G4*VECT(IPZ) + G5*HXP(3) + G6*F3)
 386 *
 387       ELSE
 388          VOUT(IX) = VECT(IX) + STEP*VECT(IPX)
 389          VOUT(IY) = VECT(IY) + STEP*VECT(IPY)
 390          VOUT(IZ) = VECT(IZ) + STEP*VECT(IPZ)
 391 *
 392       ENDIF
 393 *
 394   999 END