GEANT321/ghits/giplan.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 GIPLAN(YC,X1,X2,S1,S2,IC,XINT,SINT,PZINT,IFLAG)
  13 C.
  14 C.
  15 C.    ******************************************************************
  16 C.    *                                                                *
  17 C.    *       Calculates intersection of track (X1,X2)                 *
  18 C.    *       with plane parallel to (X-Z)                             *
  19 C.    *        The track is approximated by a cubic in the             *
  20 C.    *       track length.                                            *
  21 C.    *       To improve stability, the coordinate system              *
  22 C.    *       is shifted.                                              *
  23 C.    *       input parameters                                         *
  24 C.    *        YC    = Y COORDINATE OF PLANE                           *
  25 C.    *        X1    = X,Y,Z,XP,YP,ZP OF 1ST POINT                     *
  26 C.    *        X2    =                   2ND                           *
  27 C.    *        S1(2) = S AT 1ST(2ND) POINT                             *
  28 C.    *        IC    = 1 STRAIGHT LINE DEFINED BY X+XP                 *
  29 C.    *        IC    = 2 STRAIGHT LINE DEFINED BY X1+X2                *
  30 C.    *        IC    = 3 CUBIC MODEL                                   *
  31 C.    *                                                                *
  32 C.    *      output parameters                                         *
  33 C.    *        XINT  = X,Y,Z,XP,YP,ZP AT INTERSECTION POINT            *
  34 C.    *        SINT  = S AT INTERSECTION POINT                         *
  35 C.    *        PZINT = PHI,Z,DPHI/DR,DZ/DR                             *
  36 C.    *        IFLAG = 1 IF TRACK INTERSECTS PLANE                     *
  37 C.    *              = 0 IF NOT                                        *
  38 C.    *                                                                *
  39 C.    *      Warning : the default accuracy is 10 microns. The value   *
  40 C.    *      of EPSI must be changed for a better precision            *
  41 C.    *                                                                *
  42 C.    *    ==>Called by : <USER>, GUDIGI                               *
  43 C.    *                                                                *
  44 C.    *        Authors: R.BRUN/JJ.DUMONT from an original routine by   *
  45 C.    *       H. BOERNER  KEK  OCTOBER 1982                            *
  46 C.    *                                                                *
  47 C.    *                                                                *
  48 C.    ******************************************************************
  49 C.
  50       DIMENSION X1(6),X2(6),XINT(6),PZINT(4),A(4),B(4),C(4)
  51 C
  52       DATA MAXHIT/100/
  53       DATA EPSI/0.001/
  54 C.
  55 C.
  56 C.    ------------------------------------------------------------------
  57 C.
  58 C.
  59       IFLAG  = 1
  60       DRCTN  = 1.
  61 C
  62 C             Track crossing the plane from above or below ?
  63 C
  64       IF (X2(2).LT.X1(2))                        GO TO 5
  65       IF (YC.LT.X1(2))                           GO TO 90
  66       IF (YC.GT.X2(2))                           GO TO 90
  67       IF (IC.EQ.2) GO TO 30
  68       IF (IC.EQ.3) GO TO 7
  69 C
  70       S=S1
  71       DXDS=X1(4)
  72       DYDS=X1(5)
  73       DZDS=X1(6)
  74       BX=X1(1)-DXDS*(X1(2)-YC)/DYDS
  75       BZ=X1(3)-DZDS*(X1(2)-YC)/DYDS
  76       TRL2=(BX-X1(1))**2+(X1(2)-YC)**2+(BZ-X1(3))**2
  77       GO TO 40
  78 C
  79    5  IF (YC.LT.X2(2))                           GO TO 90
  80       IF (YC.GT.X1(2))                           GO TO 90
  81       IF(IC.EQ.2) GO TO 30
  82       DRCTN  = - 1.
  83 C
  84       IF(IC.EQ.3) GOTO 7
  85       S=S2
  86       DXDS=X2(4)
  87       DYDS=X2(5)
  88       DZDS=X2(6)
  89       BX=X2(1)-DXDS*(X2(2)-YC)/DYDS
  90       BZ=X2(3)-DZDS*(X2(2)-YC)/DYDS
  91       TRL2=(BX-X2(1))**2+(X2(2)-YC)**2+(BZ-X2(3))**2
  92       GOTO 40
  93 C
  94    30 DX=X2(1)-X1(1)
  95       DY=X2(2)-X1(2)
  96       DZ=X2(3)-X1(3)
  97       DS=SQRT(DX*DX+DY*DY+DZ*DZ)
  98       S=S1
  99       DXDS=DX/DS
 100       DYDS=DY/DS
 101       DZDS=DZ/DS
 102       BX=X1(1)-DX*(X1(2)-YC)/DY
 103       BZ=X1(3)-DZ*(X1(2)-YC)/DY
 104       TRL2=(BX-X1(1))**2+(X1(2)-YC)**2+(BZ-X1(3))**2
 105 C
 106    40 TRLEN=SQRT(TRL2)*DRCTN+S
 107       XINT(1)=BX
 108       XINT(2)=YC
 109       XINT(3)=BZ
 110       SINT=TRLEN
 111       XINT(4)=DXDS
 112       XINT(5)=DYDS
 113       XINT(6)=DZDS
 114       GO TO 200
 115 C
 116 C               Shift coordinate system such that center of gravity=0
 117 C
 118    7  IF(YC.LE.0.)                                GO TO 90
 119       SHIFTX = (X1(1) + X2(1)) * 0.5
 120       SHIFTY = (X1(2) + X2(2)) * 0.5
 121       SHIFTZ = (X1(3) + X2(3)) * 0.5
 122       SHIFTS = (S1 + S2) * 0.5
 123 C
 124 C             Only one value necessary since X1= -X2 etc...
 125 C
 126       XSHFT  = X1(1) - SHIFTX
 127       YSHFT  = X1(2) - SHIFTY
 128       ZSHFT  = X1(3) - SHIFTZ
 129       SSHFT  = S1 - SHIFTS
 130 C
 131       PABS1  = SQRT(X1(4)**2 + X1(5)**2 + X1(6)**2)
 132       PABS2  = SQRT(X2(4)**2 + X2(5)**2 + X2(6)**2)
 133       IF (PABS1.EQ.0..OR.PABS2.EQ.0.)            GO TO 90
 134 C
 135 C              Parametrize the track by a cubic through X1, X2
 136 C
 137       CALL GCUBS(SSHFT,XSHFT,X1(4)/PABS1,X2(4)/PABS2,A)
 138       CALL GCUBS(SSHFT,YSHFT,X1(5)/PABS1,X2(5)/PABS2,B)
 139       CALL GCUBS(SSHFT,ZSHFT,X1(6)/PABS1,X2(6)/PABS2,C)
 140 C
 141 C              Iterate to find the track length corresponding to
 142 C              the intersection of track and plane.
 143 C              Start at S=0. middle of the shifted interval.
 144 C
 145       DINTER = ABS(S2 - S1) * 0.5
 146       S      = 0.
 147 C
 148       DO 10 I = 1,MAXHIT
 149          Y = SHIFTY + B(1) + S * (B(2) + S * (B(3) + S * B(4)))
 150          DR=(YC-Y)*DRCTN
 151          IF (ABS(DR).LT.EPSI)                     GO TO 20
 152          DINTER = DINTER * 0.5
 153          IF (DR.LT.0.)S = S - DINTER
 154          IF (DR.GE.0.)S = S + DINTER
 155   10  CONTINUE
 156 C
 157 C              Compute intersection in original coordinates
 158 C
 159   20  CONTINUE
 160       XINT(1) = SHIFTX + A(1) + S * (A(2) + S * (A(3) + S * A(4)))
 161       XINT(2)=YC
 162       XINT(3) = SHIFTZ + C(1) + S * (C(2) + S * (C(3) + S * C(4)))
 163       XINT(4) = A(2) + S * (2. * A(3) + 3. * S * A(4))
 164       XINT(5) = B(2) + S * (2. * B(3) + 3. * S * B(4))
 165       XINT(6) = C(2) + S * (2. * C(3) + 3. * S * C(4))
 166 C
 167 C             Compute PHIHIT,ZHIT and corresponding derivatives
 168 C
 169       SINT   = S + SHIFTS
 170   200 TERM   = 1. / (XINT(4) * XINT(1) + XINT(5) * XINT(2))
 171       PZINT(1) = ATAN2(XINT(2),XINT(1))
 172       PZINT(2) = XINT(3)
 173       PZINT(3) = (XINT(1) * XINT(5) - XINT(2) * XINT(4)) * TERM / YC
 174       PZINT(4) = TERM * XINT(6) * YC
 175       RETURN
 176 C
 177   90  IFLAG  = 0
 178       END