GEANT321/ghits/gicyl.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:10  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.20  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GICYL(R,X1,X2,S1,S2,IC,XINT,SINT,PZINT,IFLAG)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *     Intersection of a Track with a Cylinder or a Plane         *
  17 C.    *     --------------------------------------------------         *
  18 C.    *                                                                *
  19 C.    *   Calculates intersection of track (x1,x2)  with cylindrical   *
  20 C.    * detector of radius R.  The track  is approximated by a cubic   *
  21 C.    * in the track length.  To  improve stability,  the coordinate   *
  22 C.    * system is shifted.                                             *
  23 C.    * R         radius of cylinder in cm                             *
  24 C.    * X1        x,y,z,xp,yp,zp of 1st point                          *
  25 C.    * X2        x,y,z,xp,yp,zp of 2nd point                          *
  26 C.    * S1(2)     S at 1st (2nd) point                                 *
  27 C.    * IC        =1 straight line defined by x+xp                     *
  28 C.    * IC        =2 straight line defined by x1+x2                    *
  29 C.    * IC        =3 cubic model                                       *
  30 C.    * XINT      x,y,z,xp,yp,zp at intersection point                 *
  31 C.    * SINT      S at intersection point                              *
  32 C.    * PZINT     phi,z,dphi/dr,dz/dr                                  *
  33 C.    * IFLAG     =1 if track intersects cylinder, =0 if not           *
  34 C.    *   Calculates  intersection  of  track  (x1,x2)   with  plane   *
  35 C.    * parallel to (X-Z).   The track is approximated by a cubic in   *
  36 C.    * the  track length.   To improve  stability,  the  coordinate   *
  37 C.    * system is shifted.                                             *
  38 C.    * YC        Y coordinate of plane                                *
  39 C.    * X1,...    as for GICYL                                         *
  40 C.    * IFLAG     =1 if track intersects plane,                        *
  41 C.    *           =0 if not                                            *
  42 C.    * Warning:  the default accuracy is  10 microns.  The value of   *
  43 C.    *         EPSI  (internal variable)   must  be  changed for  a   *
  44 C.    *         better precision.                                      *
  45 C.    *                                                                *
  46 C.    *    ==>Called by : <USER>, GUDIGI                               *
  47 C.    *                                                                *
  48 C.    *       AUTHORS:R.BRUN/JJ.DUMONT FROM AN ORIGINAL ROUTINE BY     *
  49 C.    *       H. BOERNER  KEK  OCTOBER 1982                            *
  50 C.    *                                                                *
  51 C.    *                                                                *
  52 C.    ******************************************************************
  53 C.
  54       DIMENSION X1(6),X2(6),XINT(6),PZINT(4),A(4),B(4),C(4)
  55 C
  56       DATA MAXHIT/100/
  57       DATA EPSI2/0.000001/
  58 C.
  59 C.
  60 C.    ------------------------------------------------------------------
  61 C.
  62 C.
  63       IFLAG  = 1
  64       R12    = X1(1) * X1(1) + X1(2) * X1(2)
  65       R22    = X2(1) * X2(1) + X2(2) * X2(2)
  66       R2     = R * R
  67       DR2  = R22-R2
  68 C
  69 C             TRACK CROSSING THE CYLINDER FROM INSIDE OR OUTSIDE ?
  70 C
  71       IF (R22.LT.R12)                            GO TO 5
  72       IF (R2.LT.R12)                             GO TO 90
  73       IF (R2.GT.R22)                             GO TO 90
  74       DRCTN  = 1.
  75       IF (IC.EQ.3) GO TO 7
  76 C
  77       IF(IC.EQ.2) GOTO 30
  78       S=S1
  79       DXDS=X1(4)
  80       DYDS=X1(5)
  81       DZDS=X1(6)
  82       BX=X1(1)-DXDS*S
  83       BY=X1(2)-DYDS*S
  84       BZ=X1(3)-DZDS*S
  85       GO TO 40
  86 C
  87    5  IF (R2.LT.R22)                             GO TO 90
  88       IF (R2.GT.R12)                             GO TO 90
  89       DRCTN  = - 1.
  90 C
  91       IF(IC.EQ.3) GOTO 7
  92       IF(IC.EQ.2) GOTO 30
  93       S=S2
  94       DXDS=X2(4)
  95       DYDS=X2(5)
  96       DZDS=X2(6)
  97       BX=X2(1)-DXDS*S
  98       BY=X2(2)-DYDS*S
  99       BZ=X2(3)-DZDS*S
 100       GOTO 40
 101 C
 102    30 DX=X2(1)-X1(1)
 103       DY=X2(2)-X1(2)
 104       DZ=X2(3)-X1(3)
 105       DS=SQRT(DX*DX+DY*DY+DZ*DZ)
 106       S=S1
 107       DXDS=DX/DS
 108       DYDS=DY/DS
 109       DZDS=DZ/DS
 110       BX=X1(1)-DXDS*S
 111       BY=X1(2)-DYDS*S
 112       BZ=X1(3)-DZDS*S
 113 C
 114    40 AE=DYDS*DYDS+DXDS*DXDS
 115       IF(AE.EQ.0.) GO TO 30
 116       BE=DXDS*BX+DYDS*BY
 117       CE=BY*BY+BX*BX-R2
 118       SG=SIGN(1.,DR2)
 119       XX=BE*BE-AE*CE
 120       IF(XX.LE.0.) GOTO 30
 121       TRLEN=(SG*SQRT(ABS(XX))-BE)/AE
 122       XINT(1)=DXDS*TRLEN+BX
 123       XINT(2)=DYDS*TRLEN+BY
 124       XINT(3)=DZDS*TRLEN+BZ
 125       SINT=TRLEN
 126       XINT(4)=DXDS
 127       XINT(5)=DYDS
 128       XINT(6)=DZDS
 129       GO TO 200
 130 C
 131 C             SHIFT COORDINATE SYSTEM SUCH THAT CENTER OF GRAVITY=0
 132 C
 133    7  IF(R.LE.0.)                                GO TO 90
 134       SHIFTX = (X1(1) + X2(1)) * 0.5
 135       SHIFTY = (X1(2) + X2(2)) * 0.5
 136       SHIFTZ = (X1(3) + X2(3)) * 0.5
 137       SHIFTS = (S1 + S2) * 0.5
 138 C
 139 C             ONLY ONE VALUE NECESSARY SINCE X1= -X2 ETC.
 140 C
 141       XSHFT  = X1(1) - SHIFTX
 142       YSHFT  = X1(2) - SHIFTY
 143       ZSHFT  = X1(3) - SHIFTZ
 144       SSHFT  = S1 - SHIFTS
 145 C
 146       PABS1  = SQRT(X1(4)**2 + X1(5)**2 + X1(6)**2)
 147       PABS2  = SQRT(X2(4)**2 + X2(5)**2 + X2(6)**2)
 148       IF (PABS1.EQ.0..OR.PABS2.EQ.0.)            GO TO 90
 149 C
 150 C             PARAMETRIZE THE TRACK BY A CUBIC THROUGH X1,X2
 151 C
 152       CALL GCUBS(SSHFT,XSHFT,X1(4)/PABS1,X2(4)/PABS2,A)
 153       CALL GCUBS(SSHFT,YSHFT,X1(5)/PABS1,X2(5)/PABS2,B)
 154       CALL GCUBS(SSHFT,ZSHFT,X1(6)/PABS1,X2(6)/PABS2,C)
 155 C
 156 C             ITERATE TO FIND THE TRACK LENGTH CORRESPONDING TO
 157 C             THE INTERSECTION OF TRACK AND CYLINDER.
 158 C             START AT S=0. MIDDLE OF THE SHIFTED INTERVAL.
 159 C
 160       DINTER = ABS(S2 - S1) * 0.5
 161       S      = 0.
 162 C
 163       DO 10 I = 1,MAXHIT
 164          X = SHIFTX + A(1) + S * (A(2) + S * (A(3) + S * A(4)))
 165          Y = SHIFTY + B(1) + S * (B(2) + S * (B(3) + S * B(4)))
 166          RN2    = X * X + Y * Y
 167          DR2    = (R2 - RN2) * DRCTN
 168          IF (ABS(DR2).LT.EPSI2)                     GO TO 20
 169          DINTER = DINTER * 0.5
 170          IF (DR2.LT.0.)S = S - DINTER
 171          IF (DR2.GE.0.)S = S + DINTER
 172   10  CONTINUE
 173 C
 174 C             COMPUTE INTERSECTION IN ORIGINAL COORDINATES
 175 C
 176   20  CONTINUE
 177       XINT(1) = SHIFTX + A(1) + S * (A(2) + S * (A(3) + S * A(4)))
 178       XINT(2) = SHIFTY + B(1) + S * (B(2) + S * (B(3) + S * B(4)))
 179       XINT(3) = SHIFTZ + C(1) + S * (C(2) + S * (C(3) + S * C(4)))
 180       XINT(4) = A(2) + S * (2. * A(3) + 3. * S * A(4))
 181       XINT(5) = B(2) + S * (2. * B(3) + 3. * S * B(4))
 182       XINT(6) = C(2) + S * (2. * C(3) + 3. * S * C(4))
 183 C
 184 C             COMPUTE PHIHIT,ZHIT AND CORRESPONDING DERIVATIVES
 185 C
 186       SINT   = S + SHIFTS
 187   200 TERM   = 1. / (XINT(4) * XINT(1) + XINT(5) * XINT(2))
 188       PZINT(1) = ATAN2(XINT(2),XINT(1))
 189       PZINT(2) = XINT(3)
 190       PZINT(3) = (XINT(1) * XINT(5) - XINT(2) * XINT(4)) * TERM / R
 191       PZINT(4) = TERM * XINT(6) * R
 192       RETURN
 193 C
 194   90  IFLAG  = 0
 195       END