* * $Id$ * * $Log$ * Revision 1.1.1.1 1995/10/24 10:20:04 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.44 by S.Giani *-- Author : *$ CREATE SAMCST.FOR *COPY SAMCST * *=== samcst ===========================================================* * C C******************************************************************* C SUBROUTINE SAMCST(KPROJ,EKIN,CST) #include "geant321/dblprc.inc" #include "geant321/dimpar.inc" #include "geant321/iounit.inc" C C HJM 24/10/88 C C SAMPLING OF COS(THETA) C FOR NUCLEON-PROTON ELASTIC SCATTERING C ACCORDING TO HETKFA2/BERTINI PARAMETRIZATION C C------------------------------------------------------------------- C DIMENSION DCLIN(195),DCHN(143),DCHNA(36),DCHNB(60) DIMENSION PDCI(60),PDCH(55) REAL RNDM(4) SAVE DCLIN, DCHN, DCHNA, DCHNB, PDCI, PDCH C DATA (DCLIN(I),I=1,80) / C*** DCLN ARRAY * 5.000D-01, 1.000D+00, 0.000D+00, 1.000D+00, 0.000D+00, * 4.993D-01, 9.881D-01, 5.963D-02, 9.851D-01, 5.945D-02, * 4.936D-01, 8.955D-01, 5.224D-01, 8.727D-01, 5.091D-01, * 4.889D-01, 8.228D-01, 8.859D-01, 7.871D-01, 8.518D-01, * 4.874D-01, 7.580D-01, 1.210D+00, 7.207D-01, 1.117D+00, * 4.912D-01, 6.969D-01, 1.516D+00, 6.728D-01, 1.309D+00, * 5.075D-01, 6.471D-01, 1.765D+00, 6.667D-01, 1.333D+00, * 5.383D-01, 6.054D-01, 1.973D+00, 7.059D-01, 1.176D+00, * 5.397D-01, 5.990D-01, 2.005D+00, 7.023D-01, 1.191D+00, * 5.336D-01, 6.083D-01, 1.958D+00, 6.959D-01, 1.216D+00, * 5.317D-01, 6.075D-01, 1.962D+00, 6.897D-01, 1.241D+00, * 5.300D-01, 6.016D-01, 1.992D+00, 6.786D-01, 1.286D+00, * 5.281D-01, 6.063D-01, 1.969D+00, 6.786D-01, 1.286D+00, * 5.280D-01, 5.960D-01, 2.020D+00, 6.667D-01, 1.333D+00, * 5.273D-01, 5.920D-01, 2.040D+00, 6.604D-01, 1.358D+00, * 5.273D-01, 5.862D-01, 2.069D+00, 6.538D-01, 1.385D+00/ DATA (DCLIN(I),I=81,160) / C*** DCIN ARRAY * 5.223D-01, 5.980D-01, 2.814D+00, 6.538D-01, 1.385D+00, * 5.202D-01, 5.969D-01, 2.822D+00, 6.471D-01, 1.412D+00, * 5.183D-01, 5.881D-01, 2.883D+00, 6.327D-01, 1.469D+00, * 5.159D-01, 5.866D-01, 2.894D+00, 6.250D-01, 1.500D+00, * 5.133D-01, 5.850D-01, 2.905D+00, 6.170D-01, 1.532D+00, * 5.106D-01, 5.833D-01, 2.917D+00, 6.087D-01, 1.565D+00, * 5.084D-01, 5.801D-01, 2.939D+00, 6.000D-01, 1.600D+00, * 5.063D-01, 5.763D-01, 2.966D+00, 5.909D-01, 1.636D+00, * 5.036D-01, 5.730D-01, 2.989D+00, 5.814D-01, 1.674D+00, * 5.014D-01, 5.683D-01, 3.022D+00, 5.714D-01, 1.714D+00, * 4.986D-01, 5.641D-01, 3.051D+00, 5.610D-01, 1.756D+00, * 4.964D-01, 5.580D-01, 3.094D+00, 5.500D-01, 1.800D+00, * 4.936D-01, 5.573D-01, 3.099D+00, 5.431D-01, 1.827D+00, * 4.909D-01, 5.509D-01, 3.144D+00, 5.313D-01, 1.875D+00, * 4.885D-01, 5.512D-01, 3.142D+00, 5.263D-01, 1.895D+00, * 4.857D-01, 5.437D-01, 3.194D+00, 5.135D-01, 1.946D+00/ DATA (DCLIN(I),I=161,195) / * 4.830D-01, 5.353D-01, 3.253D+00, 5.000D-01, 2.000D+00, * 4.801D-01, 5.323D-01, 3.274D+00, 4.915D-01, 2.034D+00, * 4.770D-01, 5.228D-01, 3.341D+00, 4.767D-01, 2.093D+00, * 4.738D-01, 5.156D-01, 3.391D+00, 4.643D-01, 2.143D+00, * 4.701D-01, 5.010D-01, 3.493D+00, 4.444D-01, 2.222D+00, * 4.672D-01, 4.990D-01, 3.507D+00, 4.375D-01, 2.250D+00, * 4.634D-01, 4.856D-01, 3.601D+00, 4.194D-01, 2.323D+00/ C DATA PDCI / * 4.400D+02, 1.896D-01, 1.931D-01, 1.982D-01, 1.015D-01, * 1.029D-01, 4.180D-02, 4.228D-02, 4.282D-02, 4.350D-02, * 2.204D-02, 2.236D-02, 5.900D+02, 1.433D-01, 1.555D-01, * 1.774D-01, 1.000D-01, 1.128D-01, 5.132D-02, 5.600D-02, * 6.158D-02, 6.796D-02, 3.660D-02, 3.820D-02, 6.500D+02, * 1.192D-01, 1.334D-01, 1.620D-01, 9.527D-02, 1.141D-01, * 5.283D-02, 5.952D-02, 6.765D-02, 7.878D-02, 4.796D-02, * 6.957D-02, 8.000D+02, 4.872D-02, 6.694D-02, 1.152D-01, * 9.348D-02, 1.368D-01, 6.912D-02, 7.953D-02, 9.577D-02, * 1.222D-01, 7.755D-02, 9.525D-02, 1.000D+03, 3.997D-02, * 5.456D-02, 9.804D-02, 8.084D-02, 1.208D-01, 6.520D-02, * 8.233D-02, 1.084D-01, 1.474D-01, 9.328D-02, 1.093D-01/ C DATA PDCH / * 1.000D+03, 9.453D-02, 9.804D-02, 8.084D-02, 1.208D-01, * 6.520D-02, 8.233D-02, 1.084D-01, 1.474D-01, 9.328D-02, * 1.093D-01, 1.400D+03, 1.072D-01, 7.450D-02, 6.645D-02, * 1.136D-01, 6.750D-02, 8.580D-02, 1.110D-01, 1.530D-01, * 1.010D-01, 1.350D-01, 2.170D+03, 4.004D-02, 3.013D-02, * 2.664D-02, 5.511D-02, 4.240D-02, 7.660D-02, 1.364D-01, * 2.300D-01, 1.670D-01, 2.010D-01, 2.900D+03, 1.870D-02, * 1.804D-02, 1.320D-02, 2.970D-02, 2.860D-02, 5.160D-02, * 1.020D-01, 2.400D-01, 2.250D-01, 3.370D-01, 4.400D+03, * 1.196D-03, 8.784D-03, 1.517D-02, 2.874D-02, 2.488D-02, * 4.464D-02, 8.330D-02, 2.008D-01, 2.360D-01, 3.567D-01/ C DATA (DCHN(I),I=1,90) / * 4.770D-01, 4.750D-01, 4.715D-01, 4.685D-01, 4.650D-01, * 4.610D-01, 4.570D-01, 4.550D-01, 4.500D-01, 4.450D-01, * 4.405D-01, 4.350D-01, 4.300D-01, 4.250D-01, 4.200D-01, * 4.130D-01, 4.060D-01, 4.000D-01, 3.915D-01, 3.840D-01, * 3.760D-01, 3.675D-01, 3.580D-01, 3.500D-01, 3.400D-01, * 3.300D-01, 3.200D-01, 3.100D-01, 3.000D-01, 2.900D-01, * 2.800D-01, 2.700D-01, 2.600D-01, 2.500D-01, 2.400D-01, * 2.315D-01, 2.240D-01, 2.150D-01, 2.060D-01, 2.000D-01, * 1.915D-01, 1.850D-01, 1.780D-01, 1.720D-01, 1.660D-01, * 1.600D-01, 1.550D-01, 1.500D-01, 1.450D-01, 1.400D-01, * 1.360D-01, 1.320D-01, 1.280D-01, 1.250D-01, 1.210D-01, * 1.180D-01, 1.150D-01, 1.120D-01, 1.100D-01, 1.070D-01, * 1.050D-01, 1.030D-01, 1.010D-01, 9.900D-02, 9.700D-02, * 9.550D-02, 9.480D-02, 9.400D-02, 9.200D-02, 9.150D-02, * 9.100D-02, 9.000D-02, 8.990D-02, 8.900D-02, 8.850D-02, * 8.750D-02, 8.700D-02, 8.650D-02, 8.550D-02, 8.500D-02, * 8.499D-02, 8.450D-02, 8.350D-02, 8.300D-02, 8.250D-02, * 8.150D-02, 8.100D-02, 8.030D-02, 8.000D-02, 7.990D-02/ DATA (DCHN(I),I=91,143) / * 7.980D-02, 7.950D-02, 7.900D-02, 7.860D-02, 7.800D-02, * 7.750D-02, 7.650D-02, 7.620D-02, 7.600D-02, 7.550D-02, * 7.530D-02, 7.500D-02, 7.499D-02, 7.498D-02, 7.480D-02, * 7.450D-02, 7.400D-02, 7.350D-02, 7.300D-02, 7.250D-02, * 7.230D-02, 7.200D-02, 7.100D-02, 7.050D-02, 7.020D-02, * 7.000D-02, 6.999D-02, 6.995D-02, 6.993D-02, 6.991D-02, * 6.990D-02, 6.870D-02, 6.850D-02, 6.800D-02, 6.780D-02, * 6.750D-02, 6.700D-02, 6.650D-02, 6.630D-02, 6.600D-02, * 6.550D-02, 6.525D-02, 6.510D-02, 6.500D-02, 6.499D-02, * 6.498D-02, 6.496D-02, 6.494D-02, 6.493D-02, 6.490D-02, * 6.488D-02, 6.485D-02, 6.480D-02/ C DATA DCHNA / * 6.300D+02, 7.810D-02, 1.421D-01, 1.979D-01, 2.479D-01, * 3.360D-01, 5.400D-01, 7.236D-01, 1.000D+00, 1.540D+03, * 2.225D-01, 3.950D-01, 5.279D-01, 6.298D-01, 7.718D-01, * 9.405D-01, 9.835D-01, 1.000D+00, 2.560D+03, 2.625D-01, * 4.550D-01, 5.963D-01, 7.020D-01, 8.380D-01, 9.603D-01, * 9.903D-01, 1.000D+00, 3.520D+03, 4.250D-01, 6.875D-01, * 8.363D-01, 9.163D-01, 9.828D-01, 1.000D+00, 1.000D+00, * 1.000D+00/ C DATA DCHNB / * 6.300D+02, 3.800D-02, 7.164D-02, 1.275D-01, 2.171D-01, * 3.227D-01, 4.091D-01, 5.051D-01, 6.061D-01, 7.074D-01, * 8.434D-01, 1.000D+00, 2.040D+03, 1.200D-01, 2.115D-01, * 3.395D-01, 5.295D-01, 7.251D-01, 8.511D-01, 9.487D-01, * 9.987D-01, 1.000D+00, 1.000D+00, 1.000D+00, 2.200D+03, * 1.344D-01, 2.324D-01, 3.754D-01, 5.674D-01, 7.624D-01, * 8.896D-01, 9.808D-01, 1.000D+00, 1.000D+00, 1.000D+00, * 1.000D+00, 2.850D+03, 2.330D-01, 4.130D-01, 6.610D-01, * 9.010D-01, 9.970D-01, 1.000D+00, 1.000D+00, 1.000D+00, * 1.000D+00, 1.000D+00, 1.000D+00, 3.500D+03, 3.300D-01, * 5.450D-01, 7.950D-01, 1.000D+00, 1.000D+00, 1.000D+00, * 1.000D+00, 1.000D+00, 1.000D+00, 1.000D+00, 1.000D+00/ C C--------------------------------------------------------------- CST=1.0D+0 C C* IS THE KINETIC ENERGY GREATER THAN LIMIT ? C IF (EKIN .GT. 3.5D0) RETURN C IF(KPROJ.EQ.8) GOTO 101 IF(KPROJ.EQ.1) GOTO 102 C* INVALID REACTION RETURN C-------------------------------- NP ELASTIC SCATTERING---------- 101 CONTINUE IF (EKIN.GT.0.740D+0)GOTO 1000 IF (EKIN.LT.0.300D+0)THEN C EKIN .LT. 300 MEV IDAT=1 ELSE C 300 MEV < EKIN < 740 MEV IDAT=6 END IF C ENER=EKIN IE=ABS(ENER/0.020D+0) UNIV=(ENER-IE*0.020D+0)/0.020D+0 C FORWARD/BACKWARD DECISION K=IDAT+5*IE BWFW=(DCLIN(K+5)-DCLIN(K))*UNIV + DCLIN(K) CALL GRNDM(RNDM,1) RND=RNDM(1) IF (RND.LT.BWFW)THEN VALUE2=-1.0D+0 K=K+1 ELSE VALUE2=1.0D+0 K=K+3 END IF C COEF=(DCLIN(K+5)-DCLIN(K))*UNIV + DCLIN(K) CALL GRNDM(RNDM,1) RND=RNDM(1) C IF(RND.LT.COEF)THEN CALL GRNDM(RNDM,1) CST=RNDM(1) CST=CST*VALUE2 ELSE CALL GRNDM(RNDM,4) R1=RNDM(1) R2=RNDM(2) R3=RNDM(3) R4=RNDM(4) C IF(VALUE2.GT.0.D0)THEN CST=MAX(R1,R2,R3,R4) GOTO 1500 ELSE CALL GRNDM(RNDM,1) R5=RNDM(1) C IF (IDAT.EQ.1)THEN CST=-MAX(R1,R2,R3,R4,R5) ELSE CALL GRNDM(RNDM,2) R6=RNDM(1) R7=RNDM(2) CST=-MAX(R1,R2,R3,R4,R5,R6,R7) END IF C END IF C END IF C GOTO 1500 C C******** EKIN .GT. 0.74 GEV C 1000 CONTINUE ENER=EKIN - 0.66D0 C IE=ABS(ENER/0.02D0) IE=ENER/0.02D0 EMEV=EKIN*1.D+03 C UNIV=(ENER-IE*0.020D+0)/0.020D+0 K=IE BWFW=(DCHN(K+1)-DCHN(K))*UNIV + DCHN(K) CALL GRNDM(RNDM,1) RND=RNDM(1) C FORWARD NEUTRON IF (RND.GE.BWFW)THEN DO 1200 K=10,36,9 IF (DCHNA(K).GT.EMEV) THEN UNIVE=(EMEV-DCHNA(K-9))/(DCHNA(K)-DCHNA(K-9)) CALL GRNDM(RNDM,1) UNIV=RNDM(1) DO 1100 I=1,8 II=K+I P=(DCHNA(II)-DCHNA(II-9))*UNIVE + DCHNA(II-9) C IF (P.GT.UNIV)THEN CALL GRNDM(RNDM,1) UNIV=RNDM(1) FLTI=I-UNIV GOTO(290,290,290,290,330,340,350,360) I END IF 1100 CONTINUE END IF 1200 CONTINUE C ELSE C BACKWARD NEUTRON DO 1400 K=13,60,12 IF (DCHNB(K).GT.EMEV) THEN UNIVE=(EMEV-DCHNB(K-12))/(DCHNB(K)-DCHNB(K-12)) CALL GRNDM(RNDM,1) UNIV=RNDM(1) DO 1300 I=1,11 II=K+I P=(DCHNB(II)-DCHNB(II-12))*UNIVE + DCHNB(II-12) C IF (P.GT.UNIV)THEN CALL GRNDM(RNDM,1) UNIV=RNDM(1) FLTI=I-UNIV GOTO(120,120,140,150,160,160,180,190,200,210,220) I END IF 1300 CONTINUE END IF 1400 CONTINUE END IF C 120 CST=1.0D-2*FLTI-1.0D+0 GOTO 1500 140 CST=2.0D-2*UNIV-0.98D+0 GOTO 1500 150 CST=4.0D-2*UNIV-0.96D+0 GOTO 1500 160 CST=6.0D-2*FLTI-1.16D+0 GOTO 1500 180 CST=8.0D-2*UNIV-0.80D+0 GOTO 1500 190 CST=1.0D-1*UNIV-0.72D+0 GOTO 1500 200 CST=1.2D-1*UNIV-0.62D+0 GOTO 1500 210 CST=2.0D-1*UNIV-0.50D+0 GOTO 1500 220 CST=3.0D-1*(UNIV-1.0D+0) GOTO 1500 C 290 CST=1.0D0 - 2.5D-2*FLTI GOTO 1500 330 CST=0.85D0 + 0.5D-1*UNIV GOTO 1500 340 CST=0.70D0 + 1.5D-1*UNIV GOTO 1500 350 CST=0.50D0 + 2.0D-1*UNIV GOTO 1500 360 CST=0.50D0*UNIV C 1500 RETURN C C----------------------------------- PP ELASTIC SCATTERING ------- C 102 CONTINUE EMEV=EKIN*1.D+03 C IF (EKIN.LE.0.500D0) THEN CALL GRNDM(RNDM,1) RND=RNDM(1) CST=2.0D0*RND-1.0D0 RETURN C ELSE IF (EKIN.LT.1.D0) THEN DO 2200 K=13,60,12 IF (PDCI(K).GT.EMEV) THEN UNIVE=(EMEV-PDCI(K-12))/(PDCI(K)-PDCI(K-12)) CALL GRNDM(RNDM,1) UNIV=RNDM(1) SUM=0.0D0 DO 2100 I=1,11 II=K+I SUM=SUM + (PDCI(II)-PDCI(II-12))*UNIVE + PDCI(II-12) C IF (UNIV.LT.SUM)THEN CALL GRNDM(RNDM,1) UNIV=RNDM(1) FLTI=I-UNIV GOTO(55,55,55,60,60,65,65,65,65,70,70) I END IF 2100 CONTINUE END IF 2200 CONTINUE ELSE DO 2400 K=12,55,11 IF (PDCH(K).GT.EMEV) THEN UNIVE=(EMEV-PDCH(K-11))/(PDCH(K)-PDCH(K-11)) CALL GRNDM(RNDM,1) UNIV=RNDM(1) SUM=0.D0 DO 2300 I=1,10 II=K+I SUM=SUM + (PDCH(II)-PDCH(II-11))*UNIVE + PDCH(II-11) C IF (UNIV.LT.SUM)THEN CALL GRNDM(RNDM,1) UNIV=RNDM(1) FLTI=UNIV+I GOTO(50,55,60,60,65,65,65,65,70,70) I END IF 2300 CONTINUE END IF 2400 CONTINUE END IF C 50 CST=0.4D0*UNIV GOTO 2500 55 CST=0.2D0*FLTI GOTO 2500 60 CST=0.3D0 + 0.1D0*FLTI GOTO 2500 65 CST=0.6D0 + 0.04D0*FLTI GOTO 2500 70 CST=0.78D0 + 0.02D0*FLTI C 2500 CONTINUE CALL GRNDM(RNDM,1) RND=RNDM(1) IF (RND.GT.0.5D+00)CST=-CST C RETURN END