This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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