/[PAMELA software]/gpamela/gpnd/gpgres.F
ViewVC logotype

Contents of /gpamela/gpnd/gpgres.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (show annotations) (download)
Thu Nov 15 13:55:11 2007 UTC (17 years ago) by cafagna
Branch: MAIN
CVS Tags: v4r14, HEAD
Changes since 1.1: +3 -3 lines
Debug file testgig.dat commented out in ND routines

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

  ViewVC Help
Powered by ViewVC 1.1.23