GEANT321/fluka/samcst.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:04  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.44  by  S.Giani
  11 *-- Author :
  12 *$ CREATE SAMCST.FOR
  13 *COPY SAMCST
  14 *
  15 *=== samcst ===========================================================*
  16 *
  17 C
  18 C*******************************************************************
  19 C
  20       SUBROUTINE SAMCST(KPROJ,EKIN,CST)
  21
  22 #include "geant321/dblprc.inc"
  23 #include "geant321/dimpar.inc"
  24 #include "geant321/iounit.inc"
  25 C
  26 C  HJM 24/10/88
  27 C
  28 C                 SAMPLING OF COS(THETA)
  29 C                 FOR NUCLEON-PROTON ELASTIC SCATTERING
  30 C                 ACCORDING TO HETKFA2/BERTINI PARAMETRIZATION
  31 C
  32 C-------------------------------------------------------------------
  33 C
  34       DIMENSION DCLIN(195),DCHN(143),DCHNA(36),DCHNB(60)
  35       DIMENSION PDCI(60),PDCH(55)
  36       REAL RNDM(4)
  37       SAVE DCLIN, DCHN, DCHNA, DCHNB, PDCI, PDCH
  38 C
  39       DATA (DCLIN(I),I=1,80) /
  40 C***     DCLN ARRAY
  41      *     5.000D-01,  1.000D+00,  0.000D+00,  1.000D+00,  0.000D+00,
  42      *     4.993D-01,  9.881D-01,  5.963D-02,  9.851D-01,  5.945D-02,
  43      *     4.936D-01,  8.955D-01,  5.224D-01,  8.727D-01,  5.091D-01,
  44      *     4.889D-01,  8.228D-01,  8.859D-01,  7.871D-01,  8.518D-01,
  45      *     4.874D-01,  7.580D-01,  1.210D+00,  7.207D-01,  1.117D+00,
  46      *     4.912D-01,  6.969D-01,  1.516D+00,  6.728D-01,  1.309D+00,
  47      *     5.075D-01,  6.471D-01,  1.765D+00,  6.667D-01,  1.333D+00,
  48      *     5.383D-01,  6.054D-01,  1.973D+00,  7.059D-01,  1.176D+00,
  49      *     5.397D-01,  5.990D-01,  2.005D+00,  7.023D-01,  1.191D+00,
  50      *     5.336D-01,  6.083D-01,  1.958D+00,  6.959D-01,  1.216D+00,
  51      *     5.317D-01,  6.075D-01,  1.962D+00,  6.897D-01,  1.241D+00,
  52      *     5.300D-01,  6.016D-01,  1.992D+00,  6.786D-01,  1.286D+00,
  53      *     5.281D-01,  6.063D-01,  1.969D+00,  6.786D-01,  1.286D+00,
  54      *     5.280D-01,  5.960D-01,  2.020D+00,  6.667D-01,  1.333D+00,
  55      *     5.273D-01,  5.920D-01,  2.040D+00,  6.604D-01,  1.358D+00,
  56      *     5.273D-01,  5.862D-01,  2.069D+00,  6.538D-01,  1.385D+00/
  57       DATA (DCLIN(I),I=81,160) /
  58 C***     DCIN ARRAY
  59      *     5.223D-01,  5.980D-01,  2.814D+00,  6.538D-01,  1.385D+00,
  60      *     5.202D-01,  5.969D-01,  2.822D+00,  6.471D-01,  1.412D+00,
  61      *     5.183D-01,  5.881D-01,  2.883D+00,  6.327D-01,  1.469D+00,
  62      *     5.159D-01,  5.866D-01,  2.894D+00,  6.250D-01,  1.500D+00,
  63      *     5.133D-01,  5.850D-01,  2.905D+00,  6.170D-01,  1.532D+00,
  64      *     5.106D-01,  5.833D-01,  2.917D+00,  6.087D-01,  1.565D+00,
  65      *     5.084D-01,  5.801D-01,  2.939D+00,  6.000D-01,  1.600D+00,
  66      *     5.063D-01,  5.763D-01,  2.966D+00,  5.909D-01,  1.636D+00,
  67      *     5.036D-01,  5.730D-01,  2.989D+00,  5.814D-01,  1.674D+00,
  68      *     5.014D-01,  5.683D-01,  3.022D+00,  5.714D-01,  1.714D+00,
  69      *     4.986D-01,  5.641D-01,  3.051D+00,  5.610D-01,  1.756D+00,
  70      *     4.964D-01,  5.580D-01,  3.094D+00,  5.500D-01,  1.800D+00,
  71      *     4.936D-01,  5.573D-01,  3.099D+00,  5.431D-01,  1.827D+00,
  72      *     4.909D-01,  5.509D-01,  3.144D+00,  5.313D-01,  1.875D+00,
  73      *     4.885D-01,  5.512D-01,  3.142D+00,  5.263D-01,  1.895D+00,
  74      *     4.857D-01,  5.437D-01,  3.194D+00,  5.135D-01,  1.946D+00/
  75       DATA (DCLIN(I),I=161,195) /
  76      *     4.830D-01,  5.353D-01,  3.253D+00,  5.000D-01,  2.000D+00,
  77      *     4.801D-01,  5.323D-01,  3.274D+00,  4.915D-01,  2.034D+00,
  78      *     4.770D-01,  5.228D-01,  3.341D+00,  4.767D-01,  2.093D+00,
  79      *     4.738D-01,  5.156D-01,  3.391D+00,  4.643D-01,  2.143D+00,
  80      *     4.701D-01,  5.010D-01,  3.493D+00,  4.444D-01,  2.222D+00,
  81      *     4.672D-01,  4.990D-01,  3.507D+00,  4.375D-01,  2.250D+00,
  82      *     4.634D-01,  4.856D-01,  3.601D+00,  4.194D-01,  2.323D+00/
  83 C
  84       DATA PDCI /
  85      *     4.400D+02,  1.896D-01,  1.931D-01,  1.982D-01,  1.015D-01,
  86      *     1.029D-01,  4.180D-02,  4.228D-02,  4.282D-02,  4.350D-02,
  87      *     2.204D-02,  2.236D-02,  5.900D+02,  1.433D-01,  1.555D-01,
  88      *     1.774D-01,  1.000D-01,  1.128D-01,  5.132D-02,  5.600D-02,
  89      *     6.158D-02,  6.796D-02,  3.660D-02,  3.820D-02,  6.500D+02,
  90      *     1.192D-01,  1.334D-01,  1.620D-01,  9.527D-02,  1.141D-01,
  91      *     5.283D-02,  5.952D-02,  6.765D-02,  7.878D-02,  4.796D-02,
  92      *     6.957D-02,  8.000D+02,  4.872D-02,  6.694D-02,  1.152D-01,
  93      *     9.348D-02,  1.368D-01,  6.912D-02,  7.953D-02,  9.577D-02,
  94      *     1.222D-01,  7.755D-02,  9.525D-02,  1.000D+03,  3.997D-02,
  95      *     5.456D-02,  9.804D-02,  8.084D-02,  1.208D-01,  6.520D-02,
  96      *     8.233D-02,  1.084D-01,  1.474D-01,  9.328D-02,  1.093D-01/
  97 C
  98       DATA PDCH /
  99      *     1.000D+03,  9.453D-02,  9.804D-02,  8.084D-02,  1.208D-01,
 100      *     6.520D-02,  8.233D-02,  1.084D-01,  1.474D-01,  9.328D-02,
 101      *     1.093D-01,  1.400D+03,  1.072D-01,  7.450D-02,  6.645D-02,
 102      *     1.136D-01,  6.750D-02,  8.580D-02,  1.110D-01,  1.530D-01,
 103      *     1.010D-01,  1.350D-01,  2.170D+03,  4.004D-02,  3.013D-02,
 104      *     2.664D-02,  5.511D-02,  4.240D-02,  7.660D-02,  1.364D-01,
 105      *     2.300D-01,  1.670D-01,  2.010D-01,  2.900D+03,  1.870D-02,
 106      *     1.804D-02,  1.320D-02,  2.970D-02,  2.860D-02,  5.160D-02,
 107      *     1.020D-01,  2.400D-01,  2.250D-01,  3.370D-01,  4.400D+03,
 108      *     1.196D-03,  8.784D-03,  1.517D-02,  2.874D-02,  2.488D-02,
 109      *     4.464D-02,  8.330D-02,  2.008D-01,  2.360D-01,  3.567D-01/
 110 C
 111       DATA (DCHN(I),I=1,90) /
 112      *     4.770D-01,  4.750D-01,  4.715D-01,  4.685D-01,  4.650D-01,
 113      *     4.610D-01,  4.570D-01,  4.550D-01,  4.500D-01,  4.450D-01,
 114      *     4.405D-01,  4.350D-01,  4.300D-01,  4.250D-01,  4.200D-01,
 115      *     4.130D-01,  4.060D-01,  4.000D-01,  3.915D-01,  3.840D-01,
 116      *     3.760D-01,  3.675D-01,  3.580D-01,  3.500D-01,  3.400D-01,
 117      *     3.300D-01,  3.200D-01,  3.100D-01,  3.000D-01,  2.900D-01,
 118      *     2.800D-01,  2.700D-01,  2.600D-01,  2.500D-01,  2.400D-01,
 119      *     2.315D-01,  2.240D-01,  2.150D-01,  2.060D-01,  2.000D-01,
 120      *     1.915D-01,  1.850D-01,  1.780D-01,  1.720D-01,  1.660D-01,
 121      *     1.600D-01,  1.550D-01,  1.500D-01,  1.450D-01,  1.400D-01,
 122      *     1.360D-01,  1.320D-01,  1.280D-01,  1.250D-01,  1.210D-01,
 123      *     1.180D-01,  1.150D-01,  1.120D-01,  1.100D-01,  1.070D-01,
 124      *     1.050D-01,  1.030D-01,  1.010D-01,  9.900D-02,  9.700D-02,
 125      *     9.550D-02,  9.480D-02,  9.400D-02,  9.200D-02,  9.150D-02,
 126      *     9.100D-02,  9.000D-02,  8.990D-02,  8.900D-02,  8.850D-02,
 127      *     8.750D-02,  8.700D-02,  8.650D-02,  8.550D-02,  8.500D-02,
 128      *     8.499D-02,  8.450D-02,  8.350D-02,  8.300D-02,  8.250D-02,
 129      *     8.150D-02,  8.100D-02,  8.030D-02,  8.000D-02,  7.990D-02/
 130       DATA (DCHN(I),I=91,143) /
 131      *     7.980D-02,  7.950D-02,  7.900D-02,  7.860D-02,  7.800D-02,
 132      *     7.750D-02,  7.650D-02,  7.620D-02,  7.600D-02,  7.550D-02,
 133      *     7.530D-02,  7.500D-02,  7.499D-02,  7.498D-02,  7.480D-02,
 134      *     7.450D-02,  7.400D-02,  7.350D-02,  7.300D-02,  7.250D-02,
 135      *     7.230D-02,  7.200D-02,  7.100D-02,  7.050D-02,  7.020D-02,
 136      *     7.000D-02,  6.999D-02,  6.995D-02,  6.993D-02,  6.991D-02,
 137      *     6.990D-02,  6.870D-02,  6.850D-02,  6.800D-02,  6.780D-02,
 138      *     6.750D-02,  6.700D-02,  6.650D-02,  6.630D-02,  6.600D-02,
 139      *     6.550D-02,  6.525D-02,  6.510D-02,  6.500D-02,  6.499D-02,
 140      *     6.498D-02,  6.496D-02,  6.494D-02,  6.493D-02,  6.490D-02,
 141      *     6.488D-02,  6.485D-02,  6.480D-02/
 142 C
 143       DATA DCHNA /
 144      *     6.300D+02,  7.810D-02,  1.421D-01,  1.979D-01,  2.479D-01,
 145      *     3.360D-01,  5.400D-01,  7.236D-01,  1.000D+00,  1.540D+03,
 146      *     2.225D-01,  3.950D-01,  5.279D-01,  6.298D-01,  7.718D-01,
 147      *     9.405D-01,  9.835D-01,  1.000D+00,  2.560D+03,  2.625D-01,
 148      *     4.550D-01,  5.963D-01,  7.020D-01,  8.380D-01,  9.603D-01,
 149      *     9.903D-01,  1.000D+00,  3.520D+03,  4.250D-01,  6.875D-01,
 150      *     8.363D-01,  9.163D-01,  9.828D-01,  1.000D+00,  1.000D+00,
 151      *     1.000D+00/
 152 C
 153       DATA DCHNB /
 154      *     6.300D+02,  3.800D-02,  7.164D-02,  1.275D-01,  2.171D-01,
 155      *     3.227D-01,  4.091D-01,  5.051D-01,  6.061D-01,  7.074D-01,
 156      *     8.434D-01,  1.000D+00,  2.040D+03,  1.200D-01,  2.115D-01,
 157      *     3.395D-01,  5.295D-01,  7.251D-01,  8.511D-01,  9.487D-01,
 158      *     9.987D-01,  1.000D+00,  1.000D+00,  1.000D+00,  2.200D+03,
 159      *     1.344D-01,  2.324D-01,  3.754D-01,  5.674D-01,  7.624D-01,
 160      *     8.896D-01,  9.808D-01,  1.000D+00,  1.000D+00,  1.000D+00,
 161      *     1.000D+00,  2.850D+03,  2.330D-01,  4.130D-01,  6.610D-01,
 162      *     9.010D-01,  9.970D-01,  1.000D+00,  1.000D+00,  1.000D+00,
 163      *     1.000D+00,  1.000D+00,  1.000D+00,  3.500D+03,  3.300D-01,
 164      *     5.450D-01,  7.950D-01,  1.000D+00,  1.000D+00,  1.000D+00,
 165      *     1.000D+00,  1.000D+00,  1.000D+00,  1.000D+00,  1.000D+00/
 166 C
 167 C---------------------------------------------------------------
 168         CST=1.0D+0
 169 C
 170 C*      IS THE KINETIC ENERGY GREATER THAN LIMIT ?
 171 C
 172         IF (EKIN .GT. 3.5D0) RETURN
 173 C
 174         IF(KPROJ.EQ.8) GOTO 101
 175         IF(KPROJ.EQ.1) GOTO 102
 176 C*                                INVALID REACTION
 177         RETURN
 178 C-------------------------------- NP ELASTIC SCATTERING----------
 179 101   CONTINUE
 180       IF (EKIN.GT.0.740D+0)GOTO 1000
 181       IF (EKIN.LT.0.300D+0)THEN
 182 C                                 EKIN .LT. 300 MEV
 183          IDAT=1
 184       ELSE
 185 C                                 300 MEV < EKIN < 740 MEV
 186          IDAT=6
 187       END IF
 188 C
 189       ENER=EKIN
 190       IE=ABS(ENER/0.020D+0)
 191       UNIV=(ENER-IE*0.020D+0)/0.020D+0
 192 C                                            FORWARD/BACKWARD DECISION
 193       K=IDAT+5*IE
 194       BWFW=(DCLIN(K+5)-DCLIN(K))*UNIV + DCLIN(K)
 195       CALL GRNDM(RNDM,1)
 196       RND=RNDM(1)
 197       IF (RND.LT.BWFW)THEN
 198          VALUE2=-1.0D+0
 199          K=K+1
 200       ELSE
 201          VALUE2=1.0D+0
 202          K=K+3
 203       END IF
 204 C
 205       COEF=(DCLIN(K+5)-DCLIN(K))*UNIV + DCLIN(K)
 206       CALL GRNDM(RNDM,1)
 207       RND=RNDM(1)
 208 C
 209       IF(RND.LT.COEF)THEN
 210          CALL GRNDM(RNDM,1)
 211          CST=RNDM(1)
 212          CST=CST*VALUE2
 213       ELSE
 214          CALL GRNDM(RNDM,4)
 215          R1=RNDM(1)
 216          R2=RNDM(2)
 217          R3=RNDM(3)
 218          R4=RNDM(4)
 219 C
 220          IF(VALUE2.GT.0.D0)THEN
 221             CST=MAX(R1,R2,R3,R4)
 222             GOTO 1500
 223          ELSE
 224             CALL GRNDM(RNDM,1)
 225             R5=RNDM(1)
 226 C
 227             IF (IDAT.EQ.1)THEN
 228                CST=-MAX(R1,R2,R3,R4,R5)
 229             ELSE
 230                CALL GRNDM(RNDM,2)
 231                R6=RNDM(1)
 232                R7=RNDM(2)
 233                CST=-MAX(R1,R2,R3,R4,R5,R6,R7)
 234             END IF
 235 C
 236          END IF
 237 C
 238       END IF
 239 C
 240       GOTO 1500
 241 C
 242 C********                                EKIN  .GT.  0.74 GEV
 243 C
 244 1000  CONTINUE
 245       ENER=EKIN - 0.66D0
 246 C     IE=ABS(ENER/0.02D0)
 247       IE=ENER/0.02D0
 248       EMEV=EKIN*1.D+03
 249 C
 250       UNIV=(ENER-IE*0.020D+0)/0.020D+0
 251       K=IE
 252       BWFW=(DCHN(K+1)-DCHN(K))*UNIV + DCHN(K)
 253       CALL GRNDM(RNDM,1)
 254       RND=RNDM(1)
 255 C                                        FORWARD NEUTRON
 256       IF (RND.GE.BWFW)THEN
 257          DO 1200 K=10,36,9
 258            IF (DCHNA(K).GT.EMEV) THEN
 259               UNIVE=(EMEV-DCHNA(K-9))/(DCHNA(K)-DCHNA(K-9))
 260               CALL GRNDM(RNDM,1)
 261               UNIV=RNDM(1)
 262               DO 1100 I=1,8
 263                  II=K+I
 264                  P=(DCHNA(II)-DCHNA(II-9))*UNIVE + DCHNA(II-9)
 265 C
 266                  IF (P.GT.UNIV)THEN
 267                     CALL GRNDM(RNDM,1)
 268                     UNIV=RNDM(1)
 269                     FLTI=I-UNIV
 270                     GOTO(290,290,290,290,330,340,350,360) I
 271                  END IF
 272  1100         CONTINUE
 273            END IF
 274  1200    CONTINUE
 275 C
 276       ELSE
 277 C                                        BACKWARD NEUTRON
 278          DO 1400 K=13,60,12
 279             IF (DCHNB(K).GT.EMEV) THEN
 280                UNIVE=(EMEV-DCHNB(K-12))/(DCHNB(K)-DCHNB(K-12))
 281                CALL GRNDM(RNDM,1)
 282                UNIV=RNDM(1)
 283                DO 1300 I=1,11
 284                  II=K+I
 285                  P=(DCHNB(II)-DCHNB(II-12))*UNIVE + DCHNB(II-12)
 286 C
 287                  IF (P.GT.UNIV)THEN
 288                    CALL GRNDM(RNDM,1)
 289                    UNIV=RNDM(1)
 290                    FLTI=I-UNIV
 291                    GOTO(120,120,140,150,160,160,180,190,200,210,220) I
 292                  END IF
 293  1300          CONTINUE
 294             END IF
 295  1400    CONTINUE
 296       END IF
 297 C
 298 120   CST=1.0D-2*FLTI-1.0D+0
 299       GOTO 1500
 300 140   CST=2.0D-2*UNIV-0.98D+0
 301       GOTO 1500
 302 150   CST=4.0D-2*UNIV-0.96D+0
 303       GOTO 1500
 304 160   CST=6.0D-2*FLTI-1.16D+0
 305       GOTO 1500
 306 180   CST=8.0D-2*UNIV-0.80D+0
 307       GOTO 1500
 308 190   CST=1.0D-1*UNIV-0.72D+0
 309       GOTO 1500
 310 200   CST=1.2D-1*UNIV-0.62D+0
 311       GOTO 1500
 312 210   CST=2.0D-1*UNIV-0.50D+0
 313       GOTO 1500
 314 220   CST=3.0D-1*(UNIV-1.0D+0)
 315       GOTO 1500
 316 C
 317 290   CST=1.0D0 - 2.5D-2*FLTI
 318       GOTO 1500
 319 330   CST=0.85D0 + 0.5D-1*UNIV
 320       GOTO 1500
 321 340   CST=0.70D0 + 1.5D-1*UNIV
 322       GOTO 1500
 323 350   CST=0.50D0 + 2.0D-1*UNIV
 324       GOTO 1500
 325 360   CST=0.50D0*UNIV
 326 C
 327 1500  RETURN
 328 C
 329 C-----------------------------------  PP ELASTIC SCATTERING -------
 330 C
 331  102  CONTINUE
 332       EMEV=EKIN*1.D+03
 333 C
 334       IF (EKIN.LE.0.500D0) THEN
 335          CALL GRNDM(RNDM,1)
 336          RND=RNDM(1)
 337          CST=2.0D0*RND-1.0D0
 338          RETURN
 339 C
 340       ELSE IF (EKIN.LT.1.D0) THEN
 341          DO 2200 K=13,60,12
 342             IF (PDCI(K).GT.EMEV) THEN
 343                UNIVE=(EMEV-PDCI(K-12))/(PDCI(K)-PDCI(K-12))
 344                CALL GRNDM(RNDM,1)
 345                UNIV=RNDM(1)
 346                SUM=0.0D0
 347                DO 2100 I=1,11
 348                  II=K+I
 349                  SUM=SUM + (PDCI(II)-PDCI(II-12))*UNIVE + PDCI(II-12)
 350 C
 351                  IF (UNIV.LT.SUM)THEN
 352                    CALL GRNDM(RNDM,1)
 353                    UNIV=RNDM(1)
 354                    FLTI=I-UNIV
 355                    GOTO(55,55,55,60,60,65,65,65,65,70,70) I
 356                  END IF
 357  2100          CONTINUE
 358             END IF
 359  2200    CONTINUE
 360       ELSE
 361          DO 2400 K=12,55,11
 362             IF (PDCH(K).GT.EMEV) THEN
 363               UNIVE=(EMEV-PDCH(K-11))/(PDCH(K)-PDCH(K-11))
 364               CALL GRNDM(RNDM,1)
 365               UNIV=RNDM(1)
 366               SUM=0.D0
 367               DO 2300 I=1,10
 368                 II=K+I
 369                 SUM=SUM + (PDCH(II)-PDCH(II-11))*UNIVE + PDCH(II-11)
 370 C
 371                 IF (UNIV.LT.SUM)THEN
 372                   CALL GRNDM(RNDM,1)
 373                   UNIV=RNDM(1)
 374                   FLTI=UNIV+I
 375                   GOTO(50,55,60,60,65,65,65,65,70,70) I
 376                 END IF
 377  2300         CONTINUE
 378             END IF
 379  2400    CONTINUE
 380       END IF
 381 C
 382 50    CST=0.4D0*UNIV
 383       GOTO 2500
 384 55    CST=0.2D0*FLTI
 385       GOTO 2500
 386 60    CST=0.3D0 + 0.1D0*FLTI
 387       GOTO 2500
 388 65    CST=0.6D0 + 0.04D0*FLTI
 389       GOTO 2500
 390 70    CST=0.78D0 + 0.02D0*FLTI
 391 C
 392 2500  CONTINUE
 393       CALL GRNDM(RNDM,1)
 394       RND=RNDM(1)
 395       IF (RND.GT.0.5D+00)CST=-CST
 396 C
 397       RETURN
 398       END