GEANT321/gphys/gmolie.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1996/05/02 08:17:11  ravndal
   6 * correct spheric coordinate sampling
   7 *
   8 * Revision 1.1.1.1  1995/10/24 10:21:27  cernlib
   9 * Geant
  10 *
  11 *
  12 *    EX GMOL
  13 #include "geant321/pilot.h"
  14 *CMZ :  3.21/02 29/03/94  15.41.22  by  S.Giani
  15 *-- Author :
  16       SUBROUTINE GMOLIE(OMEGA,BETA2,DIN)
  17 C.
  18 C.    ******************************************************************
  19 C.    *                                                                *
  20 C     *       Computes MOLIERE multiple scattering for a particle      *
  21 C.    *       with parameters VECT in common /GCTRAK/                  *
  22 C.    *                                                                *
  23 C.    *       This subroutine must be called with the correct values   *
  24 C.    *       of the constants OMC & CHC which depend of the medium    *
  25 C.    *                                                                *
  26 C.    *       OMC AND CHC are computed at initialisation time (GMOLI)  *
  27 C.    *       No lateral displacement of the particle with respect     *
  28 C.    *       the incident direction is included.                      *
  29 C.    *       No path length correction is included                    *
  30 C.    *                                                                *
  31 C.    *  Generation of multiple scattering according to                *
  32 C.    *  Moliere theory corrected for finite angle scattering          *
  33 C.    *                                                                *
  34 C.    *  evolved from Cern library program MLR                         *
  35 C.    *                                                                *
  36 C.    *  OMEGA & CHIC are number of scatterings and critical angle     *
  37 C.    *  of the medium for a given incident particle                   *
  38 C.    *                                                                *
  39 C.    *  COSTH and SINTH are the cosine and sine of the generated      *
  40 C.    *  scattering angle between 0 and 180 degrees                    *
  41 C.    *                                                                *
  42 C.    * THRED(NA)=reduced angles of Moliere theory                     *
  43 C.    *                                                                *
  44 C.    * F0I(NA),F1I(NA),F2I(NA)= integrale of Moliere functions        *
  45 C.    *                                                                *
  46 C.    *  4 point continued fraction interpolation is used to invert    *
  47 C.    *  the total distribution function integral                      *
  48 C.    *                                                                *
  49 C.    * XINT= argument value                                           *
  50 C.    * ARG= argument values of the table (arg,val)                    *
  51 C.    * VAL= function values of the table (arg,val)                    *
  52 C.    * THRI= the resulting interpolated function value                *
  53 C.    *                                                                *
  54 C.    *    ==>Called by : GMOLIE                                       *
  55 C.    *         Author M.S. Dixit NRCC Ottawa    *********             *
  56 C.    *                                                                *
  57 C.    *                                                                *
  58 C.    ******************************************************************
  59 C.
  60 #include "geant321/gctrak.inc"
  61 #include "geant321/gcphys.inc"
  62 #include "geant321/gcmulo.inc"
  63 #include "geant321/gconsp.inc"
  64       PARAMETER (ENEPER = 2.7182818)
  65       DIMENSION DIN(3),RNDM(3)
  66       DIMENSION TINT(40),ARG(4),VAL(4),THRED(40),F0I(40),F1I(40),F2I(40)
  67       DATA THRED/
  68      +   0.00, 0.10, 0.20, 0.30
  69      +,  0.40, 0.50, 0.60, 0.70
  70      +,  0.80, 0.90, 1.00, 1.10
  71      +,  1.20, 1.30, 1.40, 1.50
  72      +,  1.60, 1.70, 1.80, 1.90
  73      +,  2.00, 2.20, 2.40, 2.60
  74      +,  2.80, 3.00, 3.20, 3.40
  75      +,  3.60, 3.80, 4.00, 5.00
  76      +,  6.00, 7.00, 8.00, 9.00
  77      +, 10.00,11.00,12.00,13.00/
  78       DATA F0I/
  79      +  0.000000E+00 ,0.995016E-02 ,0.392106E-01 ,0.860688E-01
  80      + ,0.147856E+00 ,0.221199E+00 ,0.302324E+00 ,0.387374E+00
  81      + ,0.472708E+00 ,0.555142E+00 ,0.632121E+00 ,0.701803E+00
  82      + ,0.763072E+00 ,0.815480E+00 ,0.859142E+00 ,0.894601E+00
  83      + ,0.922695E+00 ,0.944424E+00 ,0.960836E+00 ,0.972948E+00
  84      + ,0.981684E+00 ,0.992093E+00 ,0.996849E+00 ,0.998841E+00
  85      + ,0.999606E+00 ,0.999877E+00 ,0.999964E+00 ,0.999990E+00
  86      + ,0.999998E+00 ,0.999999E+00 ,0.100000E+01 ,0.100000E+01
  87      + ,0.100000E+01 ,0.100000E+01 ,0.100000E+01 ,0.100000E+01
  88      + ,1.,1.,1.,1./
  89       DATA F1I/
  90      +  0.000000E+00,0.414985E-02,0.154894E-01,0.310312E-01
  91      + ,0.464438E-01,0.569008E-01,0.580763E-01,0.468264E-01
  92      + ,0.217924E-01,-0.163419E-01,-0.651205E-01,-0.120503E+00
  93      + ,-0.178272E+00,-0.233580E+00,-0.282442E+00,-0.321901E+00
  94      + ,-0.350115E+00,-0.366534E+00,-0.371831E+00,-0.367378E+00
  95      + ,-0.354994E+00,-0.314803E+00,-0.266539E+00,-0.220551E+00
  96      + ,-0.181546E+00,-0.150427E+00,-0.126404E+00,-0.107830E+00
  97      + ,-0.933106E-01,-0.817375E-01,-0.723389E-01,-0.436650E-01
  98      + ,-0.294700E-01,-0.212940E-01,-0.161406E-01,-0.126604E-01
  99      + ,-0.102042E-01,-0.840465E-02,-0.704261E-02,-0.598886E-02/
 100       DATA F2I/
 101      +  0.0,0.121500E-01,0.454999E-01,0.913000E-01
 102      + ,0.137300E+00,0.171400E+00,0.183900E+00,0.170300E+00
 103      + ,0.132200E+00,0.763000E-01,0.126500E-01,-0.473500E-01
 104      + ,-0.936000E-01,-0.119750E+00,-0.123450E+00,-0.106300E+00
 105      + ,-0.732800E-01,-0.312400E-01,0.128450E-01,0.528800E-01
 106      + ,0.844100E-01,0.114710E+00,0.106200E+00,0.765830E-01
 107      + ,0.435800E-01,0.173950E-01,0.695001E-03,-0.809500E-02
 108      + ,-0.117355E-01,-0.125449E-01,-0.120280E-01,-0.686530E-02
 109      + ,-0.385275E-02,-0.231115E-02,-0.147056E-02,-0.982480E-03
 110      + ,-0.682440E-03,-0.489715E-03,-0.361190E-03,-0.272582E-03/
 111 *
 112 *     ------------------------------------------------------------------
 113 *
 114 *
 115 * *** Compute Theta angle from Moliere distribution
 116 *
 117       CHIC  = CHCMOL*SQRT(STMULS)/(GETOT*BETA2)
 118       COSTH=1.
 119       SINTH=0.
 120       IF(OMEGA.LE.ENEPER)GO TO 90
 121       CNST=LOG(OMEGA)
 122       B=5.
 123       DO 10 L=1,10
 124          IF(ABS(B).LT.1.E-10)THEN
 125             B=1.E-10
 126          ENDIF
 127          DB=-(B-LOG(ABS(B))-CNST)/(1.-1./B)
 128          B=B+DB
 129          IF(ABS(DB).LE.0.0001)GO TO 20
 130    10 CONTINUE
 131       GO TO 90
 132    20 CONTINUE
 133       IF(B.LE.0.)GO TO 90
 134       BINV = 1./B
 135       TINT(1) = 0.
 136       DO 30 JA=2,4
 137          TINT(JA)=F0I(JA)+(F1I(JA)+F2I(JA)*BINV)*BINV
 138    30 CONTINUE
 139       NMAX = 4
 140    40 CONTINUE
 141       CALL GRNDM(RNDM,3)
 142       XINT=RNDM(2)
 143       DO 50 NA=3,40
 144          IF(NA.GT.NMAX) THEN
 145             TINT(NA)=F0I(NA)+(F1I(NA)+F2I(NA)*BINV)*BINV
 146             NMAX=NA
 147          ENDIF
 148          IF(XINT.LE.TINT(NA-1)) GO TO 60
 149    50 CONTINUE
 150       IF(XINT.LE.TINT(40)) THEN
 151          NA=40
 152          GOTO 60
 153       ELSE
 154          TMP=1.-(1.-B*(1.-XINT))**5
 155          IF(TMP.LE.0.)GO TO 40
 156          THRI=5./TMP
 157          GO TO 80
 158       ENDIF
 159    60 CONTINUE
 160       NA = MAX(NA-1,3)
 161       NA3 = NA-3
 162       DO 70 M=1,4
 163          NA3M=NA3+M
 164          ARG(M)=TINT(NA3M)
 165          VAL(M)=THRED(NA3M)**2
 166    70 CONTINUE
 167       F=THRED(NA)*.02
 168       CALL GMOL4(THRI,XINT,VAL,ARG,F,IER)
 169    80 CONTINUE
 170       TH=CHIC*SQRT(ABS(B*THRI))
 171       IF(TH.GT.PI)GO TO 40
 172       SINTH=SIN(TH)
 173       TEST=TH*(RNDM(3))**2
 174       IF(TEST.GT.SINTH)GO TO 40
 175       COSTH=COS(TH)
 176       GOTO 100
 177    90 CONTINUE
 178       CALL GRNDM(RNDM,1)
 179 *
 180 * *** Calculate SINE and COSINE of a random angle between 0 and 360 deg
 181 *
 182   100 PHI = RNDM(1)*TWOPI
 183 *
 184       DIN(1) = SINTH*COS(PHI)
 185       DIN(2) = SINTH*SIN(PHI)
 186       DIN(3) = COSTH
 187 *
 188       END