* * $Id$ * * $Log$ * Revision 1.1.1.1 1995/10/24 10:21:15 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.38 by S.Giani *-- Author : SUBROUTINE GNSLWD(NUCFLG,INT,NFL,TEKLOW) C C *** NEUTRON TRACKING ROUTINE FOR ENERGIES BELOW THE CUT-OFF. *** C *** TAKE ONLY ELASTIC SCATTERING, NEUTRON CAPTURE *** C *** AND NUCLEAR FISSION. *** C *** NVE 11-MAY-1988 CERN GENEVA *** C C CALLED BY : GHEISH C ORIGIN : H.FESEFELDT (ROUTINE NSLDOW 20-OCT-1987) C #include "geant321/gctrak.inc" C --- GHEISHA COMMONS --- #include "geant321/mxgkgh.inc" #include "geant321/s_consts.inc" #include "geant321/s_curpar.inc" #include "geant321/s_result.inc" #include "geant321/s_blankp.inc" #include "geant321/s_prntfl.inc" DIMENSION RNDM(2) C C --- FLAGS TO INDICATE THE NUCREC ACTION --- C NUCFLG = 0 ==> NO ACTION BY NUCREC C 1 ==> ACTION BY NUCREC ==> SPECIAL TREATMENT IN GHEISH NOPT=0 NUCFLG=0 C C --- IN ORDER TO AVOID TROUBLES CAUSED BY ARITHMETIC INCERTAINTIES, --- C --- RECALCULATE SOME QUANTITIES. TAKE KINETIC ENERGY EK AS MOST --- C --- RELEVANT QUANTITY. --- C C --- VERY LOW KINETIC ENERGY ==> NEUTRON CAPTURE --- IF (EK .LT. 1.E-9) GO TO 22 C EN=EK+ABS(AMAS) P=SQRT(ABS(EN*EN-AMAS*AMAS)) PU=SQRT(PX**2+PY**2+PZ**2) IF (PU .GE. 1.E-9) GO TO 7 C PX=0.0 PY=0.0 PZ=0.0 GO TO 22 C 7 CONTINUE PX=PX/PU PY=PY/PU PZ=PZ/PU C C --- SELECT PROCESS ACCORDING TO "INT" --- GO TO (23,23,21,22), INT C C *** NUCLEAR FISSION *** 21 CONTINUE ISTOP=1 TKIN=FISSIO(EK) GO TO 9999 C C *** NEUTRON CAPTURE *** 22 CONTINUE ISTOP=1 CALL CAPTUR(NOPT) GO TO 9999 C C *** ELASTIC AND INELASTIC SCATTERING *** 23 CONTINUE PV( 1,MXGKPV)=P*PX PV( 2,MXGKPV)=P*PY PV( 3,MXGKPV)=P*PZ PV( 4,MXGKPV)=EN PV( 5,MXGKPV)=AMAS PV( 6,MXGKPV)=NCH PV( 7,MXGKPV)=TOF PV( 8,MXGKPV)=IPART PV( 9,MXGKPV)=0.0 PV(10,MXGKPV)=USERW C C --- SPECIAL TREATMENT FOR INELASTIC SCATTERING IN HEAVY MEDIA --- IF ((INT .EQ. 2) .AND. (ATNO2 .GE. 1.5)) GO TO 29 C C *** ELASTIC SCATTERING *** 30 CONTINUE C IF (NPRT(9)) PRINT 1000 1000 FORMAT(' *GNSLWD* ELASTIC SCATTERING') C DO 24 J=4,9 PV(J,1)=PV(J,MXGKPV) 24 CONTINUE PV(10,1)=0.0 C C --- VERY SIMPLE SIMULATION OF SCATTERING ANGLE AND ENERGY --- C --- NONRELATIVISTIC APPROXIMATION WITH ISOTROPIC ANGULAR --- C --- DISTRIBUTION IN THE CMS SYSTEM --- 25 CALL GRNDM(RNDM,2) RAN=RNDM(1) COST1=-1.0+2.0*RAN EKA=1.0+2.0*COST1*ATNO2+ATNO2**2 IF(EKA.LE.0.) GOTO 25 COST=(ATNO2*COST1+1.0)/SQRT(EKA) IF (COST .LT. -1.0) COST=-1.0 IF (COST .GT. 1.0) COST=1.0 EKA=EKA/(1.0+ATNO2)**2 EK=EK*EKA EN=EK+ABS(AMAS) P=SQRT(ABS(EN*EN-AMAS*AMAS)) SINT=SQRT(ABS(1.0-COST*COST)) PHI=RNDM(2)*TWPI PV(1,2)=SINT*SIN(PHI) PV(2,2)=SINT*COS(PHI) PV(3,2)=COST CALL DEFS1(2,MXGKPV,2) PU=SQRT(PV(1,2)**2+PV(2,2)**2+PV(3,2)**2) PX=PV(1,2)/PU PY=PV(2,2)/PU PZ=PV(3,2)/PU PV(1,1)=PX*P PV(2,1)=PY*P PV(3,1)=PZ*P PV(4,1)=EN C C --- STORE BACKSCATTERED PARTICLE FOR ATNO < 4.5 --- IF (ATNO2 .GT. 4.5) GO TO 27 C IF (NPRT(9)) PRINT 1001,ATNO2 1001 FORMAT(' *GNSLWD* BACKSCATTERED PARTICLE STORED FOR ATNO ',G12.5) C PV(1,2)=PV(1,MXGKPV)-PV(1,1) PV(2,2)=PV(2,MXGKPV)-PV(2,1) PV(3,2)=PV(3,MXGKPV)-PV(3,1) CALL LENGTX(2,PP) PV(9,2)=0.0 PV(10,2)=0.0 PV(7,2)=TOF C IF (ATNO2 .GT. 3.5) GO TO 274 IF (ATNO2 .GT. 2.5) GO TO 273 IF (ATNO2 .GT. 1.5) GO TO 272 C 271 CONTINUE PV(5,2)=RMASS(14) PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2)) PV(6,2)=RCHARG(14) PV(8,2)=14.0 GO TO 275 C 272 CONTINUE PV(5,2)=RMASS(30) PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2)) PV(6,2)=RCHARG(30) PV(8,2)=30.0 GO TO 275 C 273 CONTINUE PV(5,2)=RMASS(31) PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2)) PV(6,2)=RCHARG(31) PV(8,2)=31.0 GO TO 275 C 274 CONTINUE PV(5,2)=RMASS(32) PV(4,2)=SQRT(PP*PP+PV(5,2)*PV(5,2)) PV(6,2)=RCHARG(32) PV(8,2)=32.0 C 275 CONTINUE INTCT=INTCT+1.0 CALL SETCUR(1) NTK=NTK+1 CALL SETTRK(2) GO TO 9999 C C --- PUT QUANTITIES IN COMMON /RESULT/ --- 27 CONTINUE IF (PV(10,1) .NE. 0.0) USERW=PV(10,1) SINL=PZ COSL=SQRT(ABS(1.0-SINL*SINL)) IF (ABS(COSL) .LT. 1.E-10) GO TO 28 C SINP=PY/COSL COSP=PX/COSL GO TO 9999 C 28 CONTINUE CALL GRNDM(RNDM,1) PHI=RNDM(1)*TWPI SINP=SIN(PHI) COSP=COS(PHI) GO TO 9999 C C *** INELASTIC SCATTERING ON HEAVY NUCLEI *** 29 CONTINUE C IF (NPRT(9)) PRINT 1002 1002 FORMAT(' *GNSLWD* INELASTIC SCATTERING ON HEAVY NUCLEUS') C C --- DECIDE BETWEEN SPALLATION OR SIMPLE NUCLEAR REACTION --- CALL GRNDM(RNDM,1) TEST1=RNDM(1) TEST2=4.5*(EK-0.01) IF (TEST1 .GT. TEST2) GO TO 40 C C *** SPALLATION *** C IF (NPRT(9)) PRINT 1003 1003 FORMAT(' *GNSLWD* SPALLATION') C PV( 1,MXGKPV)=P*PX PV( 2,MXGKPV)=P*PY PV( 3,MXGKPV)=P*PZ PV( 4,MXGKPV)=EN PV( 5,MXGKPV)=AMAS PV( 6,MXGKPV)=NCH PV( 7,MXGKPV)=TOF PV( 8,MXGKPV)=IPART PV( 9,MXGKPV)=0.0 PV(10,MXGKPV)=USERW C C --- FERMI-MOTION AND EVAPORATION --- TKIN=CINEMA(EK) ENP(5)=EK+TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=ENP(6)*ENP(6)-AMASQ ENP(7)=SQRT(ENP(7)) TKIN=FERMI(ENP(5)) ENP(5)=ENP(5)+TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=ENP(6)*ENP(6)-AMASQ ENP(7)=SQRT(ENP(7)) TKIN=EXNU(ENP(5)) ENP(5)=ENP(5)-TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=ENP(6)*ENP(6)-AMASQ ENP(7)=SQRT(ENP(7)) C C --- NEUTRON CASCADE --- K=2 CALL VZERO(IPA(1),MXGKCU) CALL CASN(K,INT,NFL) GO TO 9999 C 40 CONTINUE IF (NPRT(9)) PRINT 1004 1004 FORMAT(' *GNSLWD* NUCLEAR REACTION') CALL NUCREC(NOPT,1) IF (NOPT .NE. 0) NUCFLG=1 IF (NOPT .EQ. 0) GO TO 30 C 9999 CONTINUE END