COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cSampMol.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine csampmol (mediax, teta, cond)
 
subroutine csampmol0 (mediax, teta, cond)
 
subroutine cxaic2 (mediax, xc2)
 
subroutine cmolexpb (mediax, expb)
 
subroutine cprescatcalc
 

Function/Subroutine Documentation

◆ cmolexpb()

subroutine cmolexpb ( type(xsmedia), intent(in)  mediax,
real*8  expb 
)

Definition at line 201 of file cSampMol.f.

References d0.

Referenced by csampmol0().

201 ! compute exp(b) of Moliere's scattering theory.
202 ! epPreScatCalc must have been called beforehand for
203 ! each scattering.
204  use modxsecmedia
205  implicit none
206 #include "Zglobalc.h"
207 #include "Ztrack.h"
208 #include "Ztrackv.h"
209 
210 
211 ! type(epmedia):: mediax ! input. media
212  type(xsmedia),intent(in):: mediax
213  real*8 expb ! output. exp(b) of Moliere's theory.
214 
215  real*8 gbeta2, beta2, massratio2
216  common /zcmedia/ gbeta2, beta2, massratio2
217 
218  real(8):: sum, alfai2, temp
219  integer::i
220 
221  if( beta2 > 0.99 ) then ! beta > 0.995 --> regards 1
222  sum = mediax%MoliereExpb
223  else
224  sum =0.
225 ! temp = (cTrack.p.charge/137.0)**2/beta2
226  temp = (trackbefmove%p%charge/137.0)**2/beta2
227  do i = 1, mediax%noOfElem
228 
229  alfai2 = mediax%elem(i)%Z**2 *temp
230 
231  sum = sum +
232  * mediax%elem(i)%No *
233  * mediax%elem(i)%Z**0.3333*(mediax%elem(i)%Z + 1.d0)
234  * / (1.d0 + 3.327d0*alfai2)
235  enddo
236  sum = 6702.d0 *sum/mediax%A
237  endif
238 ! expb = sum *Move.dx * cTrack.p.charge**2 /beta2
239  expb = sum *intinfarray(processno)%thickness/10. *
240  * trackbefmove%p%charge**2 /beta2
241 
nodes i
real(4), dimension(:), allocatable, save temp
Definition: cNRLAtmos.f:29
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
Here is the caller graph for this function:

◆ cprescatcalc()

subroutine cprescatcalc ( )

Definition at line 246 of file cSampMol.f.

References d0, kelec, and masele.

Referenced by csampmol0().

246 ! compute g* beta^2, beta^2, (Me/M)^2
247  implicit none
248 #include "Zglobalc.h"
249 #include "Zmass.h"
250 #include "Zcode.h"
251 #include "Ztrack.h"
252 #include "Ztrackv.h"
253 
254  real*8 gbeta2, beta2, massratio2
255  common /zcmedia/ gbeta2, beta2, massratio2
256 
257  real(8):: g1, g2, g
258 
259 ! g1 = cTrack.p.fm.p(4)/cTrack.p.mass
260  g1 = trackbefmove%p%fm%p(4)/trackbefmove%p%mass
261 ! g2 = Move.Track.p.fm.p(4)/cTrack.p.mass
262  g2 = movedtrack%p%fm%p(4)/trackbefmove%p%mass
263  beta2 = 1.d0 - 1.d0/g1/g2 ! < beta^2>
264 ! gbeta2=(g1 - 1.d0/g1) * (g2 - 1.d0/g2) ! < (g*beta^2)^2 >
265 ! gbeta2 =( sqrt(g1*g2)*beta2)**2
266  gbeta2 =g1*g2*(beta2)**2
267 ! if(gbeta2 == 0.) then
268 ! g = (g1+ g2)/2.
269 ! gbeta2 = (g-1.d0/g)**2
270 ! endif
271 
272 ! if(cTrack.p.code == kelec) then
273  if(trackbefmove%p%code == kelec) then
274  massratio2= 1.
275  else
276 ! massratio2= (masele/cTrack.p.mass)**2 ! (me/m)^2
277  massratio2= (masele/trackbefmove%p%mass)**2 ! (me/m)^2
278  endif
max ptcl codes in the kelec
Definition: Zcode.h:2
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
const double masele
Definition: Zmass.h:2
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
Here is the caller graph for this function:

◆ csampmol()

subroutine csampmol ( type(xsmedia), intent(in)  mediax,
real(8), intent(out)  teta,
integer, intent(out)  cond 
)

Definition at line 7 of file cSampMol.f.

References csampmol0(), d0, false, rndc(), and true.

Referenced by cmulscat().

7  use modxsecmedia
8  implicit none
9 ! Moliere theory of multiple scattering angle.
10 ! with improvement by Bethe. rejection factor sqrt(sin teta/teta)
11 ! is included
12 !
13 ! The other is the next factor which need not be applied ?
14 ! exp(Bxc2/16.) rejection may be tried. We use Moliere when Bxc2 <1
15 ! So, if this rejection is to be tried,
16 ! since for Bxc2=1, exp(..)=1.064 so the rejection for
17 ! smaller Bxc2 must be done by
18 ! u < exp(Bxc2/16)/1.064.
19 ! Define next for this correction;
20 !cccc # define REJBYEXPBXC2
21 ! For Pb, 1GeV e, there is no difference with this correction
22 ! of without this correction.
23 !
24 !#include "Zglobalc.h"
25 !#include "Zmass.h"
26 !#include "Ztrack.h"
27 !#include "Ztrackv.h"
28 
29 ! type(epmedia):: mediax ! input. media
30  type(xsmedia),intent(in):: mediax
31  real(8),intent(out)::teta ! sampled spatial angle in radian.
32  integer,intent(out)::cond ! 0 ok. non-0. Moliere theory not applicable
33  integer::icon
34  real(8):: u
35  logical ok
36 
37  real*8 b, xc2, bxc2
38  common /zmoliere/ b, xc2, bxc2
39 
40 #ifdef REJBYEXPBXC2
41  real(8):: expbxc2
42 #endif
43  ok = .false.
44  cond = 0
45  do while (.not. ok)
46  call csampmol0(mediax, teta, icon)
47  if( icon == 0 ) then ! Moliere
48  call rndc(u)
49 !
50 #ifdef REJBYEXPBXC2
51  expbxc2 = exp(bxc2/16.0d0)/1.064d0 ! needed???
52 #endif
53 !
54  if(teta < 0.5d0 ) then
55 ! sin(teta) /teta = (teta - teta^3/6)/teta
56 ! = 1- teta^2/6. max error 0.054 %
57  ok = u**2 < (1.0d0- teta**2/6.0d0)
58 #ifdef REJBYEXPBXC2
59 ! needed ???
60  if(ok) then
61  call rndc(u)
62  ok = u < expbxc2
63  endif
64 #endif
65  elseif( teta < 3.1415 ) then
66  ok = u**2 < sin(teta)/teta
67 #ifdef REJBYEXPBXC2
68  if(ok) then
69  call rndc(u)
70  ok = u < expbxc2
71  endif
72 #endif
73  else
74  ok = .false.
75  endif
76  elseif( icon == 1) then ! GS
77 ! rejeciton not needed; GS is used
78  ok =.true.
79  else ! should use Gaussian
80  ok =.true.
81  cond = 1
82  endif
83  enddo
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
subroutine rndc(u)
Definition: rnd.f:91
subroutine csampmol0(mediax, teta, cond)
Definition: cSampMol.f:88
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon ! knockon is considered Obsolete *PhotoProd false
Definition: cblkElemag.h:7
real(4), save b
Definition: cNRLAtmos.f:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ csampmol0()

subroutine csampmol0 ( type(xsmedia), intent(in)  mediax,
real(8), intent(out)  teta,
integer, intent(out)  cond 
)

Definition at line 88 of file cSampMol.f.

References cmolexpb(), cprescatcalc(), cxaic2(), and d0.

Referenced by csampmol().

88  use modxsecmedia
89  use cmodsampmolreduceda
90  implicit none
91 ! Moliere theory of multiple scattering angle.
92 ! with improvement by Bethe.
93 !#include "Zglobalc.h"
94 !#include "Zmass.h"
95 !#include "Ztrack.h"
96 !#include "Ztrackv.h"
97 
98 ! type(epmedia):: mediax ! input. media
99  type(xsmedia),intent(in):: mediax
100 
101  real(8),intent(out)::teta ! sampled spatial angle in radian.
102  integer,intent(out)::cond ! 0 ok. Moliere used
103  ! 1 ok. G.S is used.
104  ! 2 n.g Gauss should be used.
105 
106  real*8 b, xc2, bxc2
107  common /zmoliere/ b, xc2, bxc2
108 
109  real*8 gbeta2, beta2, massratio2
110  common /zcmedia/ gbeta2, beta2, massratio2
111 
112 
113 ! *********************
114  real*8 expb
115  real(8):: sb, tetasq
116 
117  call cprescatcalc ! get g*beta^2, beta^2, (Me/M)^2 in Zcmdeia
118 ! get Xc^2
119  call cxaic2(mediax, xc2)
120 ! we don't get xa2 but directly get exp(b)
121  call cmolexpb(mediax, expb)
122  sb =log(expb) !
123 !///////////
124 ! write(*,'(a, 1p,6g14.5)')
125 ! * ' sb ', sb, beta2, Move.dt, Move.dl, cTrack.p.fm.p(4),
126 ! * Move.track.p.fm.p(4)
127 !//////////////
128 ! B -log(B) = b
129 !/// if(sb .lt. 3.395d0) then
130  if(sb .lt. 1.5d0) then
131 ! Moliere theory cannot be appliled; use Gaussian later
132 ! (almost no scattering)
133  cond = 2
134  else
135  call cmolblogb(sb, b) ! B
136  bxc2 = b*xc2
137  if(bxc2 > 1.0d0) then
138 ! normally this happens for low E, which stops soon
139 ! B so not serious.
140 ! use GS method ? --> cond=1
141  cond =2 ! use Gauss at present
142  else
143  call csampmolreduceda(b, tetasq, cond)
144  if( cond == 0 ) then ! now alwasy cond=0
145  teta = sqrt(tetasq *bxc2)
146  endif
147 !///////////
148 ! write(*,'(a, i2, 1p,9g14.5)')
149 ! * 'B ', cond, sb,beta2,B, Move.dt, Move.dl, cTrack.p.fm.p(4),
150 ! * Move.track.p.fm.p(4), Bxc2, teta
151 !//////////////
152  endif
153  endif
subroutine cprescatcalc
Definition: cSampMol.f:246
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
subroutine cxaic2(mediax, xc2)
Definition: cSampMol.f:157
subroutine cmolexpb(mediax, expb)
Definition: cSampMol.f:201
real(4), save b
Definition: cNRLAtmos.f:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cxaic2()

subroutine cxaic2 ( type(xsmedia), intent(in)  mediax,
real*8  xc2 
)

Definition at line 157 of file cSampMol.f.

References d0, and kelec.

Referenced by csampmol0().

157 ! compute xc2 of Moliere's theory of scattring.
158 ! epPreScatCalc must have been called beforehand for
159 ! each scattering.
160  use modxsecmedia
161  implicit none
162 #include "Zglobalc.h"
163 #include "Ztrack.h"
164 #include "Ztrackv.h"
165 #include "Zcode.h"
166 
167 
168 ! type(epmedia):: mediax ! input. media
169  type(xsmedia),intent(in):: mediax
170 
171  real*8 xc2 ! output. Xc^2 in radian^2
172 
173  real*8 gbeta2, beta2, massratio2
174  common /zcmedia/ gbeta2, beta2, massratio2
175 
176 
177 !
178 ! integeral 0 to leng of 1/(beta**4 E**2)
179 ! is approximated as leng/(beta1**2 gamma1 beta2**2 gamma2)
180 ! Note: bata**2 *gamma = gamma - 1/gamma
181 !
182 
183 
184 ! if( cTrack.p.code == kelec) then
185  if( trackbefmove%p%code == kelec) then
186  xc2 = mediax%MoliereForXc2
187  else
188 ! xc2 = 0.6011d0*cTrack.p.charge**2 * massratio2 *
189  xc2 = 0.6011d0*trackbefmove%p%charge**2 * massratio2 *
190  * mediax%Z2/mediax%A
191  endif
192 ! Move.dx is for EPICS.
193 ! xc2 = xc2*Move.dx/gbeta2 ! dx is in g/cm2. no need to consider rhoc.
194  ! real length is dx/(rho*rhoc)
195 ! thickness is kg/m2;* 1000/10000 = 1/10 g/cm2
196  xc2 = xc2*intinfarray(processno)%thickness/10./gbeta2 !
197 
max ptcl codes in the kelec
Definition: Zcode.h:2
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
Here is the caller graph for this function: