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

Annotation of /gpamela/gpnd/gpgres.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide 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 pamela 1.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 cafagna 1.2 * open(77,file='testgig.dat',access='append')
33 pamela 1.1 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 cafagna 1.2 * write(77,*) Number_N, step, en_gamma
71     * close(77)
72 pamela 1.1 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