GEANT321/gheisha/coscat.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:58  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 COSCAT
  13 C
  14 C *** MOMENTUM GENERATION FOR COHERENT ELASTIC SCATTERING ***
  15 C *** NVE 13-JUL-1988 CERN GENEVA ***
  16 C
  17 C ORIGIN : H.FESEFELDT (03-DEC-1986)
  18 C
  19 C APPROXIMATION OF BESSEL FUNCTION FOR TETA(LAB)<=20 DEG.
  20 C IS USED . THE NUCLEAR RADIUS IS TAKEN AS R=1.25*E-13*(A)**1/3FM
  21 C
  22 #include "geant321/s_defcom.inc"
  23 #include "geant321/s_coscom.inc"
  24 #include "geant321/s_kginit.inc"
  25 C
  26       EXTERNAL FCTCOS
  27       DIMENSION FF(20),ATNOX(3)
  28       DIMENSION RNDM(1)
  29 C
  30       DATA ATNOX/9.,56.,207./
  31 C
  32 C --- INITIALIZATION INDICATED BY KGINIT(14) ---
  33       IF (KGINIT(14) .NE. 0) GO TO 10
  34       KGINIT(14)=1
  35 C
  36       IF(.NOT.NPRT(10)) GOTO 10
  37       WRITE(NEWBCD,2001)
  38  2001 FORMAT(1H0,'DS/DT FOR COHERENT ELASTIC SCATTERING')
  39       DO 3 L=1,3
  40       WRITE(NEWBCD,2003) ATNOX(L),P
  41  2003 FORMAT(1H0,'CALCULATED CROSS SECTIONS FOR A=',F5.1,' AND P=',F8.2)
  42       DO 2 I=1,20
  43       TETA=(I-1)*PI/360.
  44       T=2.*P**2*(1.-COS(TETA*1.D0))
  45       IF(ATNOX(L).GT.62.) GOTO 4
  46       FF(I)=TWPI*ATNOX(L)**1.63*EXP(-14.5D0*ATNOX(L)**0.65*T)
  47      *     +TWPI*1.4*ATNOX(L)**0.33*EXP(-10.D0*T)
  48       GOTO 2
  49     4 FF(I)=TWPI*ATNOX(L)**1.33*EXP(-60.0D0*ATNOX(L)**0.33*T)
  50      *     +TWPI*0.4*ATNOX(L)**0.40*EXP(-10.D0*T)
  51     2 CONTINUE
  52       WRITE(NEWBCD,2004) FF
  53  2004 FORMAT(1H ,10E12.3)
  54     3 CONTINUE
  55    10 IF(P.LT.0.01) GO TO 9999
  56       IF(ATNO2.LT.0.5) GO TO 9999
  57       IER(46)=IER(46)+1
  58       RAN=RANRES(DUM)
  59       CALL VZERO(IPA(1),MXGKCU)
  60       IPA(1)=IPART
  61       IF(ATNO2.GT.62.) GOTO 11
  62       AA=ATNO2**1.63
  63       BB=14.5*ATNO2**0.66
  64       CC=1.4*ATNO2**0.33
  65       DD=10.
  66       AA=AA/BB
  67       CC=CC/DD
  68       RR=(AA+CC)*RAN
  69       GOTO 12
  70    11 AA=ATNO2**1.33
  71       BB=60.*ATNO2**0.33
  72       CC=0.4*ATNO2**0.40
  73       DD=10.
  74       AA=AA/BB
  75       CC=CC/DD
  76       RR=(AA+CC)*RAN
  77    12 T1=-LOG(RAN)/BB
  78       T2=-LOG(RAN)/DD
  79       EPS=0.001
  80       IND1=10
  81       CALL RTMI(T,VAL,FCTCOS,T1,T2,EPS,IND1,IER1)
  82       IF(IER1.EQ.0) GOTO 14
  83       T=0.25*(3.*T1+T2)
  84       IER(68)=IER(68)+1
  85    14 CALL GRNDM(RNDM,1)
  86       PHI=RNDM(1)*TWPI
  87       RR=0.5*T/P**2
  88       IF(RR.GT.1.) RR=0.
  89       COST=1.-RR
  90 *     SINT=SQRT(MAX((1.-COST)*(1.+COST),0.))
  91       SINT=SQRT(MAX(RR*(2.-RR),0.))
  92       IF(SINT.NE.0.) THEN
  93       PV( 1,MXGKPV-1)=P*PX
  94       PV( 2,MXGKPV-1)=P*PY
  95       PV( 3,MXGKPV-1)=P*PZ
  96       PV( 4,MXGKPV-1)=EN
  97       PV( 5,MXGKPV-1)=AMAS
  98       PV( 6,MXGKPV-1)=NCH
  99       PV( 7,MXGKPV-1)=TOF
 100       PV( 8,MXGKPV-1)=IPART
 101       PV( 9,MXGKPV-1)=0.
 102       PV(10,MXGKPV-1)=USERW
 103       PV(1,1)=P*SINT*SIN(PHI)
 104       PV(2,1)=P*SINT*COS(PHI)
 105       PV(3,1)=P*COST
 106       PV(4,1)=EN
 107       PV(5,1)=AMAS
 108       PV(6,1)=NCH
 109       PV(7,1)=TOF
 110       PV(8,1)=IPART
 111       PV(9,1)=0.
 112       PV(10,1)=0.
 113       CALL DEFS1(1,MXGKPV-1,1)
 114       SINL1=SINL
 115       COSL1=COSL
 116       SINP1=SINP
 117       COSP1=COSP
 118       CALL SETCUR(1)
 119       ELSE
 120       SINL1=SINL
 121       COSL1=COSL
 122       SINP1=SINP
 123       COSP1=COSP
 124       ENDIF
 125       IF(NPRT(4))
 126      *WRITE(NEWBCD,1004) AMAS,P,SINL1,COSL1,SINP1,COSP1,SINL,COSL,
 127      *                   SINP,COSP,T1,T,T2,IER1
 128 C
 129  1004 FORMAT(1H ,'COHERENT ELASTIC SCATTERING    MASS ',F8.3,' MOMENTUM
 130      * ',F8.3/1H ,'DIRECTION ',4F10.4,' CHANGED TO ',4F10.4/
 131      *1H ,'T1,T,T2 ',3E10.3,' IER1 ',I2)
 132 C
 133  9999 CONTINUE
 134       END