GEANT321/giface/gnslwd.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:15  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GNSLWD(NUCFLG,INT,NFL,TEKLOW)
  13 C
  14 C *** NEUTRON TRACKING ROUTINE FOR ENERGIES BELOW THE CUT-OFF. ***
  15 C *** TAKE ONLY ELASTIC SCATTERING, NEUTRON CAPTURE            ***
  16 C *** AND NUCLEAR FISSION.                                     ***
  17 C *** NVE 11-MAY-1988 CERN GENEVA ***
  18 C
  19 C CALLED BY : GHEISH
  20 C ORIGIN : H.FESEFELDT (ROUTINE NSLDOW 20-OCT-1987)
  21 C
  22 #include "geant321/gctrak.inc"
  23 C --- GHEISHA COMMONS ---
  24 #include "geant321/mxgkgh.inc"
  25 #include "geant321/s_consts.inc"
  26 #include "geant321/s_curpar.inc"
  27 #include "geant321/s_result.inc"
  28 #include "geant321/s_blankp.inc"
  29 #include "geant321/s_prntfl.inc"
  30       DIMENSION RNDM(2)
  31 C
  32 C --- FLAGS TO INDICATE THE NUCREC ACTION ---
  33 C NUCFLG = 0 ==> NO ACTION BY NUCREC
  34 C          1 ==> ACTION BY NUCREC ==> SPECIAL TREATMENT IN GHEISH
  35       NOPT=0
  36       NUCFLG=0
  37 C
  38 C --- IN ORDER TO AVOID TROUBLES CAUSED BY ARITHMETIC INCERTAINTIES, ---
  39 C --- RECALCULATE SOME QUANTITIES. TAKE KINETIC ENERGY EK AS MOST ---
  40 C --- RELEVANT QUANTITY. ---
  41 C
  42 C --- VERY LOW KINETIC ENERGY ==> NEUTRON CAPTURE ---
  43       IF (EK .LT. 1.E-9) GO TO 22
  44 C
  45       EN=EK+ABS(AMAS)
  46       P=SQRT(ABS(EN*EN-AMAS*AMAS))
  47       PU=SQRT(PX**2+PY**2+PZ**2)
  48       IF (PU .GE. 1.E-9) GO TO 7
  49 C
  50       PX=0.0
  51       PY=0.0
  52       PZ=0.0
  53       GO TO 22
  54 C
  55  7    CONTINUE
  56       PX=PX/PU
  57       PY=PY/PU
  58       PZ=PZ/PU
  59 C
  60 C --- SELECT PROCESS ACCORDING TO "INT" ---
  61       GO TO (23,23,21,22), INT
  62 C
  63 C *** NUCLEAR FISSION ***
  64  21   CONTINUE
  65       ISTOP=1
  66       TKIN=FISSIO(EK)
  67       GO TO 9999
  68 C
  69 C *** NEUTRON CAPTURE ***
  70  22   CONTINUE
  71       ISTOP=1
  72       CALL CAPTUR(NOPT)
  73       GO TO 9999
  74 C
  75 C *** ELASTIC AND INELASTIC SCATTERING ***
  76  23   CONTINUE
  77       PV( 1,MXGKPV)=P*PX
  78       PV( 2,MXGKPV)=P*PY
  79       PV( 3,MXGKPV)=P*PZ
  80       PV( 4,MXGKPV)=EN
  81       PV( 5,MXGKPV)=AMAS
  82       PV( 6,MXGKPV)=NCH
  83       PV( 7,MXGKPV)=TOF
  84       PV( 8,MXGKPV)=IPART
  85       PV( 9,MXGKPV)=0.0
  86       PV(10,MXGKPV)=USERW
  87 C
  88 C --- SPECIAL TREATMENT FOR INELASTIC SCATTERING IN HEAVY MEDIA ---
  89       IF ((INT .EQ. 2) .AND. (ATNO2 .GE. 1.5)) GO TO 29
  90 C
  91 C *** ELASTIC SCATTERING ***
  92  30   CONTINUE
  93 C
  94       IF (NPRT(9)) PRINT 1000
  95  1000 FORMAT(' *GNSLWD* ELASTIC SCATTERING')
  96 C
  97       DO 24 J=4,9
  98       PV(J,1)=PV(J,MXGKPV)
  99  24   CONTINUE
 100       PV(10,1)=0.0
 101 C
 102 C --- VERY SIMPLE SIMULATION OF SCATTERING ANGLE AND ENERGY ---
 103 C --- NONRELATIVISTIC APPROXIMATION WITH ISOTROPIC ANGULAR ---
 104 C --- DISTRIBUTION IN THE CMS SYSTEM ---
 105   25  CALL GRNDM(RNDM,2)
 106       RAN=RNDM(1)
 107       COST1=-1.0+2.0*RAN
 108       EKA=1.0+2.0*COST1*ATNO2+ATNO2**2
 109       IF(EKA.LE.0.) GOTO 25
 110       COST=(ATNO2*COST1+1.0)/SQRT(EKA)
 111       IF (COST .LT. -1.0) COST=-1.0
 112       IF (COST .GT. 1.0) COST=1.0
 113       EKA=EKA/(1.0+ATNO2)**2
 114       EK=EK*EKA
 115       EN=EK+ABS(AMAS)
 116       P=SQRT(ABS(EN*EN-AMAS*AMAS))
 117       SINT=SQRT(ABS(1.0-COST*COST))
 118       PHI=RNDM(2)*TWPI
 119       PV(1,2)=SINT*SIN(PHI)
 120       PV(2,2)=SINT*COS(PHI)
 121       PV(3,2)=COST
 122       CALL DEFS1(2,MXGKPV,2)
 123       PU=SQRT(PV(1,2)**2+PV(2,2)**2+PV(3,2)**2)
 124       PX=PV(1,2)/PU
 125       PY=PV(2,2)/PU
 126       PZ=PV(3,2)/PU
 127       PV(1,1)=PX*P
 128       PV(2,1)=PY*P
 129       PV(3,1)=PZ*P
 130       PV(4,1)=EN
 131 C
 132 C --- STORE BACKSCATTERED PARTICLE FOR ATNO < 4.5 ---
 133       IF (ATNO2 .GT. 4.5) GO TO 27
 134 C
 135       IF (NPRT(9)) PRINT 1001,ATNO2
 136  1001 FORMAT(' *GNSLWD* BACKSCATTERED PARTICLE STORED FOR ATNO ',G12.5)
 137 C
 138       PV(1,2)=PV(1,MXGKPV)-PV(1,1)
 139       PV(2,2)=PV(2,MXGKPV)-PV(2,1)
 140       PV(3,2)=PV(3,MXGKPV)-PV(3,1)
 141       CALL LENGTX(2,PP)
 142       PV(9,2)=0.0
 143       PV(10,2)=0.0
 144       PV(7,2)=TOF
 145 C
 146       IF (ATNO2 .GT. 3.5) GO TO 274
 147       IF (ATNO2 .GT. 2.5) GO TO 273
 148       IF (ATNO2 .GT. 1.5) GO TO 272
 149 C
 150  271  CONTINUE
 151       PV(5,2)=RMASS(14)
 152       PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2))
 153       PV(6,2)=RCHARG(14)
 154       PV(8,2)=14.0
 155       GO TO 275
 156 C
 157  272  CONTINUE
 158       PV(5,2)=RMASS(30)
 159       PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2))
 160       PV(6,2)=RCHARG(30)
 161       PV(8,2)=30.0
 162       GO TO 275
 163 C
 164  273  CONTINUE
 165       PV(5,2)=RMASS(31)
 166       PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2))
 167       PV(6,2)=RCHARG(31)
 168       PV(8,2)=31.0
 169       GO TO 275
 170 C
 171  274  CONTINUE
 172       PV(5,2)=RMASS(32)
 173       PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2))
 174       PV(6,2)=RCHARG(32)
 175       PV(8,2)=32.0
 176 C
 177  275  CONTINUE
 178       INTCT=INTCT+1.0
 179       CALL SETCUR(1)
 180       NTK=NTK+1
 181       CALL SETTRK(2)
 182       GO TO 9999
 183 C
 184 C --- PUT QUANTITIES IN COMMON /RESULT/ ---
 185  27   CONTINUE
 186       IF (PV(10,1) .NE. 0.0) USERW=PV(10,1)
 187       SINL=PZ
 188       COSL=SQRT(ABS(1.0-SINL*SINL))
 189       IF (ABS(COSL) .LT. 1.E-10) GO TO 28
 190 C
 191       SINP=PY/COSL
 192       COSP=PX/COSL
 193       GO TO 9999
 194 C
 195  28   CONTINUE
 196       CALL GRNDM(RNDM,1)
 197       PHI=RNDM(1)*TWPI
 198       SINP=SIN(PHI)
 199       COSP=COS(PHI)
 200       GO TO 9999
 201 C
 202 C *** INELASTIC SCATTERING ON HEAVY NUCLEI ***
 203  29   CONTINUE
 204 C
 205       IF (NPRT(9)) PRINT 1002
 206  1002 FORMAT(' *GNSLWD* INELASTIC SCATTERING ON HEAVY NUCLEUS')
 207 C
 208 C --- DECIDE BETWEEN SPALLATION OR SIMPLE NUCLEAR REACTION ---
 209       CALL GRNDM(RNDM,1)
 210       TEST1=RNDM(1)
 211       TEST2=4.5*(EK-0.01)
 212       IF (TEST1 .GT. TEST2) GO TO 40
 213 C
 214 C *** SPALLATION ***
 215 C
 216       IF (NPRT(9)) PRINT 1003
 217  1003 FORMAT(' *GNSLWD* SPALLATION')
 218 C
 219       PV( 1,MXGKPV)=P*PX
 220       PV( 2,MXGKPV)=P*PY
 221       PV( 3,MXGKPV)=P*PZ
 222       PV( 4,MXGKPV)=EN
 223       PV( 5,MXGKPV)=AMAS
 224       PV( 6,MXGKPV)=NCH
 225       PV( 7,MXGKPV)=TOF
 226       PV( 8,MXGKPV)=IPART
 227       PV( 9,MXGKPV)=0.0
 228       PV(10,MXGKPV)=USERW
 229 C
 230 C --- FERMI-MOTION AND EVAPORATION ---
 231       TKIN=CINEMA(EK)
 232       ENP(5)=EK+TKIN
 233 C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
 234       IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
 235       ENP(6)=ENP(5)+ABS(AMAS)
 236       ENP(7)=ENP(6)*ENP(6)-AMASQ
 237       ENP(7)=SQRT(ENP(7))
 238       TKIN=FERMI(ENP(5))
 239       ENP(5)=ENP(5)+TKIN
 240 C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
 241       IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
 242       ENP(6)=ENP(5)+ABS(AMAS)
 243       ENP(7)=ENP(6)*ENP(6)-AMASQ
 244       ENP(7)=SQRT(ENP(7))
 245       TKIN=EXNU(ENP(5))
 246       ENP(5)=ENP(5)-TKIN
 247 C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
 248       IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
 249       ENP(6)=ENP(5)+ABS(AMAS)
 250       ENP(7)=ENP(6)*ENP(6)-AMASQ
 251       ENP(7)=SQRT(ENP(7))
 252 C
 253 C --- NEUTRON CASCADE ---
 254       K=2
 255       CALL VZERO(IPA(1),MXGKCU)
 256       CALL CASN(K,INT,NFL)
 257       GO TO 9999
 258 C
 259  40   CONTINUE
 260       IF (NPRT(9)) PRINT 1004
 261  1004 FORMAT(' *GNSLWD* NUCLEAR REACTION')
 262       CALL NUCREC(NOPT,1)
 263       IF (NOPT .NE. 0) NUCFLG=1
 264       IF (NOPT .EQ. 0) GO TO 30
 265 C
 266  9999 CONTINUE
 267       END