1 |
C*****This subroutine fills the common NeutronsGR: Number_N, Neut_En(3), |
2 |
C teta_n(3), fi_n(3), |
3 |
C*****where |
4 |
C*****Number_N - number of neutrons in case of Gigantic resonance |
5 |
C (Number_N=0, if no); |
6 |
C*****NeutE_n(i) - energy of each generated neutron, MeV |
7 |
C*****teta_n(i) - teta (polar) angle in reference coordinate system |
8 |
C*****fi_n(i) - fi (azimuth) angle in reference coordinate system Input |
9 |
C parameters from GEANT: Step - GEANT step, cm En-gamma - gamma |
10 |
C energy in tungsten, GeV px_gamma, py_gamma, pz_gamma - gamma |
11 |
C momentum projections, GeV |
12 |
subroutine GPGRES(Step, En_g, px_g, py_g, pz_g) |
13 |
|
14 |
integer i |
15 |
real En_exc |
16 |
real cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n |
17 |
|
18 |
C Input parameters |
19 |
* Energy of gamma, in GeV |
20 |
real En_gamma |
21 |
* Gamma momentum projections, multiplied by c, in GeV |
22 |
real px_gamma, py_gamma, pz_gamma |
23 |
* GEANT step, in cm |
24 |
real Step |
25 |
|
26 |
common/bind/En_bind(3) |
27 |
#include "gpgneut.inc" |
28 |
c |
29 |
* |
30 |
* Let's trasnform energy and momentum in MeV. Remember, GEANT is GeV |
31 |
* |
32 |
* open(77,file='testgig.dat',access='append') |
33 |
En_gamma=En_g/1000. |
34 |
px_gamma=px_g/1000. |
35 |
py_gamma=py_g/1000. |
36 |
pz_gamma=pz_g/1000. |
37 |
* |
38 |
* Exciting energy of nucleus in Gigantic resonance, MeV |
39 |
En_exc=En_gamma |
40 |
Number_N=0 |
41 |
do i=1,3 |
42 |
Neut_En(i)=0. |
43 |
teta_n(i)=0. |
44 |
fi_n(i)=0. |
45 |
enddo |
46 |
|
47 |
C Simulate number of neutrons in Gygantic resonance process |
48 |
C (Number_n=0 if no interaction) |
49 |
call GRESsimulate(En_gamma, Step, Number_N) |
50 |
|
51 |
if (Number_N.gt.0) then |
52 |
cos_teta=pz_gamma/sqrt(px_gamma**2+py_gamma**2+pz_gamma**2) |
53 |
sin_teta=sqrt(1.-cos_teta**2) |
54 |
cos_fi=px_gamma/sqrt(px_gamma**2+py_gamma**2) |
55 |
sin_fi=py_gamma/sqrt(px_gamma**2+py_gamma**2) |
56 |
|
57 |
do i=1, Number_N |
58 |
if (En_exc.gt.En_bind(i)) then |
59 |
call SpectrumSpace(Number_N, i, En_exc, |
60 |
& cos_teta, sin_teta, cos_fi, sin_fi, |
61 |
& En_n, alfa_n, beta_n) |
62 |
Neut_En(i)=En_n |
63 |
teta_n(i)=alfa_n |
64 |
fi_n(i)=beta_n |
65 |
En_exc=En_exc-En_bind(i)-En_n |
66 |
endif |
67 |
enddo |
68 |
|
69 |
endif |
70 |
* write(77,*) Number_N, step, en_gamma |
71 |
* close(77) |
72 |
RETURN |
73 |
end |
74 |
|
75 |
c************************************************************************ |
76 |
************************************************************************* |
77 |
C*****Simulate generated number of |
78 |
C*****neutrons |
79 |
C****************************** |
80 |
subroutine GRESsimulate(En_gamma, Step, Number_N) |
81 |
c Output parameter: Number_N (as default Number_N=0) |
82 |
real En_gamma, Step |
83 |
integer Number_N |
84 |
|
85 |
real Sigma |
86 |
real*4 tmp |
87 |
real Lambda1, Lambda2, Lambda3 |
88 |
|
89 |
#include "gpgsigma.inc" |
90 |
|
91 |
if ((En_gamma.gt.9.097).and.(En_gamma.lt.28.608)) then |
92 |
CCCCC Calculate Lambda1 |
93 |
Sigma=-1. |
94 |
do i=1,63 |
95 |
if ((En_gamma.ge.energy1(i)).and.(En_gamma.lt.energy1(i+1))) |
96 |
& Sigma=cross1(i)+(cross1(i+1)-cross1(i))* |
97 |
& (En_gamma-energy1(i))/ |
98 |
& (energy1(i+1)-energy1(i)) |
99 |
enddo |
100 |
if (Sigma.le.0.) then |
101 |
Lambda1=1000. |
102 |
else |
103 |
1 call GRNDM(tmp,1) |
104 |
if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 1 |
105 |
Lambda1=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22) |
106 |
endif |
107 |
CCCCCCCalculate Lambda2 |
108 |
Sigma=-1. |
109 |
do i=1,61 |
110 |
if ((En_gamma.ge.energy2(i)).and.(En_gamma.lt.energy2(i+1))) |
111 |
& Sigma=cross2(i)+(cross2(i+1)-cross2(i))* |
112 |
& (En_gamma-energy2(i))/ |
113 |
& (energy2(i+1)-energy2(i)) |
114 |
enddo |
115 |
if (Sigma.le.0.) then |
116 |
Lambda2=1000. |
117 |
else |
118 |
2 call GRNDM(tmp,1) |
119 |
if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 2 |
120 |
Lambda2=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22) |
121 |
endif |
122 |
|
123 |
CCCCCCCalculate Lambda3 |
124 |
Sigma=-1. |
125 |
do i=1,26 |
126 |
if ((En_gamma.ge.energy3(i)).and.(En_gamma.lt.energy3(i+1))) |
127 |
& Sigma=cross3(i)+(cross3(i+1)-cross3(i))* |
128 |
& (En_gamma-energy3(i))/ |
129 |
& (energy3(i+1)-energy3(i)) |
130 |
enddo |
131 |
if (Sigma.le.0.) then |
132 |
Lambda3=1000. |
133 |
else |
134 |
3 call GRNDM(tmp,1) |
135 |
if ((tmp.lt.1.E-8).or.(tmp.gt.(1.-1.E-8))) goto 3 |
136 |
Lambda3=-log(1.-tmp)/(Sigma*1.E-27*6.31*1.E22) |
137 |
endif |
138 |
|
139 |
if (Lambda1.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=1 |
140 |
if (Lambda2.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=2 |
141 |
if (Lambda3.eq.MIN(Lambda1, Lambda2, Lambda3, Step)) Number_N=3 |
142 |
endif |
143 |
RETURN |
144 |
end |
145 |
|
146 |
c************************************************************************ |
147 |
c************************************************************************ |
148 |
******Simulate neutron energy and flight direction*********************** |
149 |
|
150 |
subroutine SpectrumSpace(Number_N, Num, En_exc, |
151 |
& cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n) |
152 |
c Output parameters: En_n, alfa_n, beta_n |
153 |
common/bind/En_bind(3) |
154 |
|
155 |
real tmp(2) |
156 |
real Pi, Temperature, lev_dens, cos_alfa, sin_alfa, beta |
157 |
real cos_alfa_n, cos_beta_n, sin_beta_n |
158 |
integer Number_N, Num |
159 |
|
160 |
real En_exc |
161 |
real cos_teta, sin_teta, cos_fi, sin_fi, En_n, alfa_n, beta_n |
162 |
|
163 |
Pi=4.*atan(1.) |
164 |
if (Number_N.eq.1) then |
165 |
En_n=En_exc-En_bind(1) |
166 |
|
167 |
else if ((Number_N.eq.2).and.(Num.eq.1)) then |
168 |
C************************************************************************ |
169 |
C SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book |
170 |
C is: V.C. Barashenkov, V.D. Toneev, "The interactions of high |
171 |
C energy particles and atomic nuclei with nuclei", in Russian |
172 |
C Calculation of density levels parameter (see (6.31) from reference |
173 |
C book) |
174 |
lev_dens=0.085*184. |
175 |
C Calculation of nucleus temperature (see (6.30) from reference book) |
176 |
C (see (6.64) and comments after) |
177 |
Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) |
178 |
|
179 |
C Neutron evaporation energy spectrum simulation by exception method |
180 |
4 call GRNDM(tmp,2) |
181 |
C The average binding energy of neutron |
182 |
x_tmp=(En_exc-En_bind(1)-En_bind(2))*tmp(1) |
183 |
y_tmp=(exp(-1.)/Temperature)*tmp(2) |
184 |
if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature))) |
185 |
& goto 4 |
186 |
En_n=x_tmp |
187 |
else if ((Number_N.eq.2).and.(Num.eq.2)) then |
188 |
En_n=En_exc-En_bind(2) |
189 |
else if ((Number_N.eq.3).and.(Num.eq.1)) then |
190 |
C************************************************************************ |
191 |
C SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book is: |
192 |
C V.C. Barashenkov, V.D. Toneev, "The interactions of high energy |
193 |
C particles and atomic nuclei with nuclei", in Russian Calculation |
194 |
C of density levels parameter (see (6.31) from reference book) |
195 |
lev_dens=0.085*184. |
196 |
C Calculation of nucleus temperature (see (6.30) from reference book) |
197 |
C (see (6.64) and comments after) |
198 |
Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) |
199 |
|
200 |
C Neutron evaporation energy spectrum simulation by exception method |
201 |
5 call GRNDM(tmp,2) |
202 |
C Average binding energy of neutron |
203 |
x_tmp=(En_exc-En_bind(1)-En_bind(2)-En_bind(3))*tmp(1) |
204 |
y_tmp=(exp(-1.)/Temperature)*tmp(2) |
205 |
if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature))) |
206 |
& goto 5 |
207 |
En_n=x_tmp |
208 |
|
209 |
else if ((Number_N.eq.3).and.(Num.eq.2)) then |
210 |
C************************************************************************ |
211 |
C SIMULATION OF NEUTRON EVAPORATION SPECTRUM The reference book |
212 |
C is: V.C. Barashenkov, V.D. Toneev, "The interactions of high |
213 |
C energy particles and atomic nuclei with nuclei", in Russian |
214 |
C Calculation of density levels parameter (see (6.31) from reference |
215 |
C book) |
216 |
lev_dens=0.085*183. |
217 |
C Calculation of nucleus temperature (see (6.30) from reference book) |
218 |
C (see (6.64) and comments after) |
219 |
Temperature=sqrt((En_exc-En_bind(Num))/lev_dens) |
220 |
|
221 |
C Neutron evaporation energy spectrum simulation by exception method |
222 |
6 call GRNDM(tmp,2) |
223 |
C Average binding energy of neutron |
224 |
x_tmp=(En_exc-En_bind(2)-En_bind(3))*tmp(1) |
225 |
y_tmp=(exp(-1.)/Temperature)*tmp(2) |
226 |
if (y_tmp.gt.((x_tmp/Temperature**2)*exp(-x_tmp/Temperature))) |
227 |
& goto 6 |
228 |
En_n=x_tmp |
229 |
|
230 |
else if ((Number_N.eq.3).and.(Num.eq.3)) then |
231 |
En_n=En_exc-En_bind(3) |
232 |
|
233 |
endif |
234 |
C************************************************************************ |
235 |
C Neutron direction simulation in gamma coordinate system |
236 |
call GRNDM(tmp,2) |
237 |
cos_alfa=-1.+2.*tmp(1) |
238 |
beta=2.*Pi*tmp(2) |
239 |
sin_alfa=sqrt(1.-cos_alfa**2) |
240 |
C Return to the GEANT reference coordinate system |
241 |
cos_alfa_n=-sin_alfa*cos(beta)*(sin_teta)+cos_alfa*(cos_teta) |
242 |
cos_beta_n=sin_alfa*cos(beta)*(cos_teta*cos_fi)- |
243 |
& sin_alfa*sin(beta)*(sin_fi)+cos_alfa*(sin_teta*cos_fi) |
244 |
sin_beta_n=sin_alfa*cos(beta)*(cos_teta*sin_fi)+ |
245 |
& sin_alfa*sin(beta)*(cos_fi)+cos_alfa*(sin_teta*sin_fi) |
246 |
if (abs(cos_beta_n).lt.1.E-6) cos_beta_n=0. |
247 |
if (abs(sin_beta_n).lt.1.E-6) sin_beta_n=0. |
248 |
if ((cos_beta_n.eq.0.).and.(sin_beta_n.eq.0.)) beta_n=0. |
249 |
|
250 |
alfa_n=acos(cos_alfa_n) |
251 |
if ((cos_beta_n.ge.0.).and.(sin_beta_n.gt.0.)) |
252 |
& beta_n=atan(sin_beta_n/cos_beta_n) |
253 |
if ((cos_beta_n.ge.0.).and.(sin_beta_n.lt.0.)) |
254 |
& beta_n=2.*Pi+atan(sin_beta_n/cos_beta_n) |
255 |
if ((cos_beta_n.lt.0.).and.(sin_beta_n.le.0.)) |
256 |
& beta_n=Pi+atan(sin_beta_n/cos_beta_n) |
257 |
if ((cos_beta_n.lt.0.).and.(sin_beta_n.ge.0.)) |
258 |
& beta_n=Pi+atan(sin_beta_n/cos_beta_n) |
259 |
|
260 |
RETURN |
261 |
END |