| fe4da5cc |
1 | * |
| 2 | * $Id$ |
| 3 | * |
| 4 | * $Log$ |
| 5 | * Revision 1.1.1.1 1995/10/24 10:20:53 cernlib |
| 6 | * Geant |
| 7 | * |
| 8 | * |
| 9 | #include "geant321/pilot.h" |
| 10 | *CMZ : 3.21/02 29/03/94 15.41.30 by S.Giani |
| 11 | *-- Author : |
| 12 | ********************************************************************* |
| 13 | ***** GNPGO1 ******************************************************** |
| 14 | * |
| 15 | * GNPGO1 ... 15-AUG-1991 |
| 16 | * Version 1.1 |
| 17 | * Rolf Nierhaus |
| 18 | * |
| 19 | ********************************************************************* |
| 20 | * |
| 21 | * Copyright CERN, Geneva 1991 - Copyright and any other |
| 22 | * appropriate legal protection of these computer programs and |
| 23 | * associated documentation reserved in all countries of the |
| 24 | * world. |
| 25 | * |
| 26 | ********************************************************************* |
| 27 | * |
| 28 | * Subroutine GNPGO1 is called by GNPGON for the computation |
| 29 | * of SNXT, the distance from a point P along a track T to a |
| 30 | * boundary surface of a Geant volume V of shape PGON. The point |
| 31 | * P is inside the volume V. |
| 32 | * |
| 33 | * V is generally a composite volume consisting of several |
| 34 | * sections. The sections have boundary surfaces orthogonal to |
| 35 | * the Z-axis. Each section consists generally of several |
| 36 | * sectors. Each sector is an "elementary" convex volume. This |
| 37 | * package assumes it is either a hexahedron or a pentahedron. If |
| 38 | * it is a pentahedron, it has 6 vertices, of which two are on |
| 39 | * the Z-axis. All sectors of the same section are congruent. |
| 40 | * Each section has the same number of sectors. |
| 41 | * |
| 42 | * We describe each surface by 6 parameters: the first three |
| 43 | * are the coordinates of a point on the surface |
| 44 | * XS(I),YS(I),ZS(I), the other three are the components of the |
| 45 | * normal vector of the surface XN(I),YN(I),ZN(I). I is the index |
| 46 | * of the surface. We consider only one sector at a time, and the |
| 47 | * number of boundary surfaces is never larger then 6. Each |
| 48 | * surface divides the space into two regions: the positive |
| 49 | * region and the negative region. We choose the direction of the |
| 50 | * normal vectors of the boundary surfaces such that the bounded |
| 51 | * volume is within the positive region of each surface, that is, |
| 52 | * the normal vector is pointing to the inside of the volume. |
| 53 | * |
| 54 | ***** Subroutine GNPGO1 *************************** 15-AUG-1991 ***** |
| 55 | SUBROUTINE GNPGO1(X,P,SNXT) |
| 56 | #if !defined(CERNLIB_SINGLE) |
| 57 | IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
| 58 | #endif |
| 59 | REAL X(6),P(49),SNXT |
| 60 | #if defined(CERNLIB_SINGLE) |
| 61 | PARAMETER (F=0.01745329251994330,TPI=6.283185307179586) |
| 62 | #endif |
| 63 | #if !defined(CERNLIB_SINGLE) |
| 64 | PARAMETER (F=0.01745329251994330D0,TPI=6.283185307179586D0) |
| 65 | #endif |
| 66 | PARAMETER (ONE=1,HALF=ONE/2,ZERO=0.) |
| 67 | DIMENSION XS(6), YS(6), ZS(6), XN(6), YN(6), ZN(6) |
| 68 | LOGICAL FLAG, FLG |
| 69 | XP=X(1) |
| 70 | YP=X(2) |
| 71 | ZP=X(3) |
| 72 | XD=X(4) |
| 73 | YD=X(5) |
| 74 | ZD=X(6) |
| 75 | IMAX=P(4)-.5 |
| 76 | * IMAX -> number of Z-sections |
| 77 | JMAX=P(3)+.5 |
| 78 | * JMAX -> number of Phi-sectors |
| 79 | SNXT=0. |
| 80 | 10 CONTINUE |
| 81 | * Find current elementary volume |
| 82 | IF (ZP.LE.P(5)) RETURN |
| 83 | * Current point (XP,YP,ZP) is below first section. |
| 84 | DO 20 II=1,IMAX |
| 85 | IF (ZP.LT.P(5+3*II)) THEN |
| 86 | I=II |
| 87 | GO TO 30 |
| 88 | ENDIF |
| 89 | 20 CONTINUE |
| 90 | RETURN |
| 91 | * Current point (XP,YP,ZP) is above last section. |
| 92 | 30 CONTINUE |
| 93 | IF (XP.EQ.0..AND.YP.EQ.0.) XP=1.E-20 |
| 94 | PHI=ATAN2(YP,XP) |
| 95 | IF (PHI.LT.0.) PHI=PHI+TPI |
| 96 | P1=F*P(1) |
| 97 | PHI1=PHI-P1 |
| 98 | IF (PHI1.LT.0.) PHI1=PHI1+TPI |
| 99 | IF (PHI1.GE.TPI) PHI1=PHI1-TPI |
| 100 | IF (JMAX.EQ.1) THEN |
| 101 | IF (ABS(PHI1-TPI).LT.1D-7) PHI1=0. |
| 102 | ENDIF |
| 103 | J=PHI1*P(3)/(F*P(2))+ONE |
| 104 | IF (P(2).EQ.360.) THEN |
| 105 | IF (J.LT.1) THEN |
| 106 | J=J+JMAX |
| 107 | ELSEIF (JMAX.LT.J) THEN |
| 108 | J=J-JMAX |
| 109 | END IF |
| 110 | END IF |
| 111 | IF (JMAX.LT.J.OR.J.LT.1) RETURN |
| 112 | * Current point is outside Phi-range. |
| 113 | C***** Code Expanded From Routine: GNPGO2 |
| 114 | * GNPGO2 finds the vector distance to the boundary surface |
| 115 | * of the current elementary volume. |
| 116 | * I is Z-section, J is Phi-sector. |
| 117 | * GNPGO2 calls GNPGO4 five or six times for the storage of |
| 118 | * the surface coefficients of its boundary surfaces. |
| 119 | * |
| 120 | INDEX = 2 + 3*I |
| 121 | Z1 = P(INDEX) |
| 122 | D1N = P(INDEX+1) |
| 123 | D1X = P(INDEX+2) |
| 124 | Z2 = P(INDEX+3) |
| 125 | D2N = P(INDEX+4) |
| 126 | D2X = P(INDEX+5) |
| 127 | ZM = HALF*(Z1 + Z2) |
| 128 | P11X = F*(P(1)+(J-1)*P(2)/JMAX) |
| 129 | P2 = F*(P(1)+J*P(2)/JMAX) |
| 130 | PP = HALF*(P11X + P2) |
| 131 | COSP = COS(PP) |
| 132 | SINP = SIN(PP) |
| 133 | DMX = HALF*(D1X + D2X) |
| 134 | DMN = HALF*(D1N + D2N) |
| 135 | THX = ATAN((D2X - D1X)/(Z2 - Z1)) |
| 136 | COSTHX = COS(THX) |
| 137 | SINTHX = SIN(THX) |
| 138 | XNN = -COSP*COSTHX |
| 139 | YNN = -SINP*COSTHX |
| 140 | C***** Code Expanded From Routine: GNPGO4 |
| 141 | * Store surface coefficients |
| 142 | XS(5) = DMX*COSP |
| 143 | YS(5) = DMX*SINP |
| 144 | ZS(5) = ZM |
| 145 | XN(5) = XNN |
| 146 | YN(5) = YNN |
| 147 | ZN(5) = SINTHX |
| 148 | C***** End of Code Expanded From Routine: GNPGO4 |
| 149 | C***** Code Expanded From Routine: GNPGO9 |
| 150 | * Logical function GNPGO9 returns TRUE if the point |
| 151 | * (XP,YP,ZP) is within the positive region of the surface with |
| 152 | * index I. This is the case if the scalar product of |
| 153 | * (XP-XS,YP-YS,ZP-ZS) and (XN,YN,ZN) is positive (or zero). |
| 154 | RESULT=(XP-XS(5))*XN(5)+(YP-YS(5))*YN(5)+(ZP-ZS(5))*ZN(5) |
| 155 | FLG = 0. .LE. RESULT |
| 156 | IF (.NOT.FLG) GO TO 50 |
| 157 | ISMAX = 5 |
| 158 | IF (DMN .NE. 0.) THEN |
| 159 | ISMAX = 6 |
| 160 | THN = ATAN((D2N - D1N)/(Z2 - Z1)) |
| 161 | COSTHN = COS(THN) |
| 162 | SINTHN = SIN(THN) |
| 163 | XNN = COSP*COSTHN |
| 164 | YNN = SINP*COSTHN |
| 165 | C***** Code Expanded From Routine: GNPGO4 |
| 166 | XS(6) = DMN*COSP |
| 167 | YS(6) = DMN*SINP |
| 168 | ZS(6) = ZM |
| 169 | XN(6) = XNN |
| 170 | YN(6) = YNN |
| 171 | ZN(6) = -SINTHN |
| 172 | C***** End of Code Expanded From Routine: GNPGO4 |
| 173 | C***** Code Expanded From Routine: GNPGO9 |
| 174 | RESULT=(XP-XS(6))*XN(6)+(YP-YS(6))*YN(6)+(ZP-ZS(6))*ZN(6) |
| 175 | FLG = 0. .LE. RESULT |
| 176 | IF (.NOT.FLG) GO TO 50 |
| 177 | ENDIF |
| 178 | C***** Code Expanded From Routine: GNPGO4 |
| 179 | XS(1) = ZERO |
| 180 | YS(1) = ZERO |
| 181 | ZS(1) = Z1 |
| 182 | XN(1) = ZERO |
| 183 | YN(1) = ZERO |
| 184 | ZN(1) = ONE |
| 185 | C***** End of Code Expanded From Routine: GNPGO4 |
| 186 | C***** Code Expanded From Routine: GNPGO4 |
| 187 | XS(2) = ZERO |
| 188 | YS(2) = ZERO |
| 189 | ZS(2) = Z2 |
| 190 | XN(2) = ZERO |
| 191 | YN(2) = ZERO |
| 192 | ZN(2) = -ONE |
| 193 | C***** End of Code Expanded From Routine: GNPGO4 |
| 194 | C***** Code Expanded From Routine: GNPGO4 |
| 195 | XS(3) = ZERO |
| 196 | YS(3) = ZERO |
| 197 | ZS(3) = ZM |
| 198 | XN(3) = -SIN(P11X) |
| 199 | YN(3) = COS(P11X) |
| 200 | ZN(3) = ZERO |
| 201 | C***** End of Code Expanded From Routine: GNPGO4 |
| 202 | C***** Code Expanded From Routine: GNPGO4 |
| 203 | XS(4) = ZERO |
| 204 | YS(4) = ZERO |
| 205 | ZS(4) = ZM |
| 206 | XN(4) = SIN(P2) |
| 207 | YN(4) = -COS(P2) |
| 208 | ZN(4) = ZERO |
| 209 | C***** End of Code Expanded From Routine: GNPGO4 |
| 210 | C***** Code Expanded From Routine: GNPGO5 |
| 211 | * Vector distance to volume boundary |
| 212 | SNXT1 = 1.E10 |
| 213 | DO 40 IS = 1, ISMAX |
| 214 | C***** Code Expanded From Routine: GNPGO7 |
| 215 | * To find the distance from a point (XP,YP,ZP) along a |
| 216 | * track with direction cosines (XD,YD,ZD) to a surface |
| 217 | * (XS,YS,ZS)(XN,YN,ZN), we compute first the scalar product of |
| 218 | * the vector (XS-XP,YS-YP,ZS-ZP) with the normal vector |
| 219 | * (XN,YN,ZN), then the scalar product of the vectors (XD,YD,ZD) |
| 220 | * and (XN,YN,ZN). The first scalar product is the shortest |
| 221 | * distance from the point to the plane, the second scalar |
| 222 | * product is the cosine of the angle between the track and the |
| 223 | * plane normal. The quotient is the vector distance. If this |
| 224 | * vector distance is positive (or zero) we set the logical |
| 225 | * variable FLAG TRUE. GNPGO7 is called with three parameters |
| 226 | * I,FLAG and DIST. I is the index of the surface, and DIST is |
| 227 | * the vector distance if FLAG is TRUE. |
| 228 | SPPMSN = (XP - XS(IS))*XN(IS) + (YP - YS(IS))*YN(IS) + (ZP - |
| 229 | + ZS (IS))*ZN(IS) |
| 230 | SPDN = XD*XN(IS) + YD*YN(IS) + ZD*ZN(IS) |
| 231 | IF (SPDN .EQ. 0.) THEN |
| 232 | DIST1 = 0. |
| 233 | ELSE |
| 234 | DIST1 = -(SPPMSN + .0001)/SPDN |
| 235 | ENDIF |
| 236 | FLAG = 0. .LT. DIST1 |
| 237 | C***** End of Code Expanded From Routine: GNPGO7 |
| 238 | IF (FLAG) SNXT1 = MIN(DIST1,SNXT1) |
| 239 | 40 CONTINUE |
| 240 | 50 CONTINUE |
| 241 | C***** End of Code Expanded From Routine: GNPGO2 |
| 242 | IF (FLG) THEN |
| 243 | SNXT=SNXT+SNXT1 |
| 244 | XP=XP+SNXT1*XD |
| 245 | YP=YP+SNXT1*YD |
| 246 | ZP=ZP+SNXT1*ZD |
| 247 | * The current point (XP,YP,ZP) is propagated along the track |
| 248 | * to the boundary of the current elementary volume. |
| 249 | GO TO 10 |
| 250 | ENDIF |
| 251 | END |
git://git.uio.no / u / mrichter / AliRoot.git / blame