COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cintePhoton.f
Go to the documentation of this file.
1  subroutine cintephoton
2  implicit none
3 #include "Ztrack.h"
4 #include "Ztrackv.h"
5  integer i
6  character*70 msg
7 
8 !///////////////////
9 ! to see possibility of photon projectile
10 ! with Primary is 'gamma' ; see last
11 ! part of this file.
12 ! call cpredpmjet( MovedTrack.p, Pwork, Nproduced)
13 !//////////
14 
15  if(intinfarray(processno)%process .eq. 'pair') then
16  call cpair
17  elseif(intinfarray(processno)%process .eq. 'compt') then
18  call ccompt
19  elseif(intinfarray(processno)%process .eq. 'photoe') then
20  call cphotoee
21  elseif(intinfarray(processno)%process .eq. 'photop') then
22  call cphotop
23  elseif(intinfarray(processno)%process .eq. 'cohs') then
24  call ccohs
25  elseif(intinfarray(processno)%process .eq. 'mpair') then
26  call cmpair
27  else
28  write(msg, *) ' process for photon; Proc#=',
29  * processno,
30  * intinfarray(processno)%process, ' undef'
31  call cerrormsg(msg, 1)
32  write(0, '("energy =",g12.4, "where=",i4, " w=",g12.3)')
33  * movedtrack%p%fm%p(4), movedtrack%where, movedtrack%wgt
34  write(0, * ) " coszenith=", movedtrack%vec%coszenith
35  write(0, *) ' MoveStat=',movestat, 'No Of inte=',
36  * numberofinte
37  do i = 1, numberofinte
38  write(0,*) intinfarray(i)%process, ' dt=',
39  * intinfarray(i)%thickness,intinfarray(i)%length
40  enddo
41  write(0,*)
42  write(0, * ) " (dep,h)B//A==", trackbefmove%pos%depth,
43  * trackbefmove%pos%height,
44  * movedtrack%pos%depth, movedtrack%pos%height
45  stop
46  endif
47  end
48 ! *******************
49  subroutine cpair
50 ! *******************
51  implicit none
52 
53 #include "Zcode.h"
54 #include "Zmass.h"
55 #include "Ztrackp.h"
56 #include "Ztrack.h"
57 #include "Ztrackv.h"
58 
59 !
60  real*8 e1, e2, u, eg, cs, sn
61  real*8 teta, teta1, teta2, cos1, sin1, cos2, sin2
62  integer ica
63  real*8 den, cvh2den
64  type(track)::aTrack
65  type(coord)::dc
66  type(coord)::dce
67  real*8 temp
68 !
69  eg = movedtrack%p%fm%p(4)
70  if(lpmeffect .and. eg .gt. lpmpairemin) then
71 ! den = cthick2den(TrackBefMove.pos.depth)
72  den = cvh2den(trackbefmove%pos%height) ! better
73  call cpairerglpm(eg, den, e1)
74  else
75  call cpairenergy(eg, e1)
76  endif
77 
78  e2=eg - e1
79  if( e1 .gt. e2) then
80 ! store higher energy ptcl later
81  temp = e1
82  e1= e2
83  e2 = temp
84  endif
85 ! now e1 < e2
86 
87 ! assign charge for e1
88  call rndc(u)
89  if(u .gt. .5) then
90  ica=1
91  else
92  ica=-1
93  endif
94 !
95  atrack = movedtrack
96 ! if(eg .lt. 100.e-3) then
97 ! take pair angle if < 100 MeV
98  if(eg .lt. 10) then
99 ! take pair angle if < 10GeV
100 ! call cPairAng(1, eg, e1, teta1) ! teta1 < pi/2
101  call cpairang(e1, masele, teta1)
102  if(teta1 .lt. 0.03d0) then
103  cos1 = 1. - teta1**2/2
104  sin1 = teta1
105  else
106  cos1 = cos(teta1)
107  sin1 = sin(teta1)
108  endif
109 !
110  sin2 = sin1 * sqrt( (e1**2-masele**2)/(e2**2-masele**2) )
111  if(sin2 .lt. 0.03d0) then
112  cos2 = 1.- sin2**2/2
113  else
114  cos2 = sqrt(1.d0 - sin2**2)
115  endif
116  call kcossn(cs, sn)
117 
118 ! teta = masele/eg
119 ! teta1=teta* e2/eg
120 ! teta2=teta* e1/eg
121 ! cos1=1. - teta1**2/2
122 ! cos2=1. - teta2**2/2
123 ! sample direction cos. of 1st
124 ! sin1=teta1
125 ! sin2=teta2
126 !
127  dc%r(1) = cs * sin1
128  dc%r(2) = sn * sin1
129  dc%r(3)= cos1
130  call ctransvectz(movedtrack%vec%w, dc, dce)
131  call cmkptc(kelec, 0, ica, atrack%p)
132  atrack%p%fm%p(4) = e1
133  call csetdircos(dce, atrack)
134  call ce2p(atrack)
135  nproduced = nproduced + 1
136  pwork(nproduced) = atrack%p
137 ! another electron (higher E)
138 
139  dc%r(1) = -cs*sin2
140  dc%r(2) = -sn*sin2
141  dc%r(3) = cos2
142  call ctransvectz(movedtrack%vec%w, dc, dce)
143  atrack%p%fm%p(4) = e2
144  call cmkptc(kelec, 0, -ica, atrack%p)
145  call csetdircos(dce, atrack)
146  call ce2p(atrack)
147  nproduced = nproduced + 1
148  pwork(nproduced) = atrack%p
149  else
150 ! neglect pair angle
151  atrack%p%fm%p(4) = e1
152  call ce2p(atrack)
153  call cmkptc(kelec, 0, ica, atrack%p)
154  nproduced = nproduced + 1
155  pwork(nproduced) = atrack%p
156 !
157  atrack%p%fm%p(4) = e2
158  call ce2p(atrack)
159  call cmkptc(kelec, 0, -ica, atrack%p)
160  nproduced = nproduced + 1
161  pwork(nproduced) = atrack%p
162  endif
163  end
164 ! ***********
165  subroutine ccompt
166 ! ***********
167  implicit none
168 !---- include '../Particle/Zcode.h'
169 #include "Zcode.h"
170 !---- include 'Ztrack.h'
171 #include "Ztrack.h"
172 !---- include 'Ztrackv.h'
173 #include "Ztrackv.h"
174 !
175  type(track)::aTrack
176  real*8 eg, e1, cs, sn, cosg, cose
177  real*8 sine, tmp, sing
178  type(coord)::dc
179  type(coord)::dce
180  type(coord)::dcg
181 !
182  call ccomptea(movedtrack%p%fm%p(4), eg, e1, cosg, cose)
183 
184  call kcossn(cs,sn)
185 ! electron direction
186  tmp=max(1.d0-cose*cose, 0.d0)
187  sine=sqrt(tmp)
188  dc%r(1)=cs*sine
189  dc%r(2)=sn*sine
190  dc%r(3)=cose
191  call ctransvectz(movedtrack%vec%w, dc, dce)
192  atrack = movedtrack
193  call cmkptc(kelec, 0, -1, atrack%p)
194  atrack%p%fm%p(4) = e1
195  call csetdircos(dce, atrack)
196  call ce2p(atrack)
197  nproduced = nproduced + 1
198  pwork(nproduced) = atrack%p
199 ! gamma dicrection
200  tmp=max(1.d0-cosg*cosg, 0.d0)
201  sing=sqrt(tmp)
202  dc%r(1) = -cs*sing
203  dc%r(2) = -sn*sing
204  dc%r(3) = cosg
205  call ctransvectz(movedtrack%vec%w, dc, dcg)
206  atrack%p%fm%p(4) = eg
207  call cmkptc(kphoton, kcasg, 0, atrack%p)
208  call csetdircos(dcg, atrack)
209  call ce2p(atrack)
210  nproduced = nproduced + 1
211  pwork(nproduced) = atrack%p
212 
213 
214  end
215 ! ***********
216  subroutine cphotop
217 ! ***********
218  use modxsecmedia
219  implicit none
220 #include "Zcode.h"
221 #include "Ztrack.h"
222 #include "Ztrackv.h"
223 #include "Ztrackp.h"
224 !
225 !
226  integer ngen
227  character(8):: whichcode
228  nproduced = 0
229  if( howphotop == 0 ) then ! will not happen
230  ngen = 1
231  pwork(nproduced+1) = movedtrack%p
232  elseif( howphotop == 1 ) then
233  whichcode = "sofia"
234  elseif( howphotop == 2 ) then
235  if( movedtrack%p%fm%p(4) < 2.5 ) then
236  whichcode ="current"
237  else
238  whichcode ="sofia"
239  endif
240  elseif( howphotop == 3 ) then
241  if( movedtrack%p%fm%p(4) < 2.5 ) then
242  whichcode ="sofia"
243  else
244  whichcode ="current"
245  endif
246  elseif( howphotop == 4 ) then
247  whichcode = "current"
248  else
249  write(0,*) 'HowPhotoP =', howphotop, ' invalid '
250  stop
251  endif
252 
253  if( whichcode == "sofia" ) then
254  call csofia( targetnucleonno, targetprotonno,
255  * movedtrack%p, pwork(nproduced+1), ngen)
256  elseif( whichcode == "current" ) then
257  call cgphad(targetnucleonno, targetprotonno,
258  * movedtrack%p, pwork(nproduced+1), ngen)
259  else
260  write(0,*) ' setting mistake of whichcode=',
261  * whichcode, ' in cphotop'
262  stop
263  endif
264  nproduced = nproduced + ngen
265 
266  end
267 ! ****************
268  subroutine ccohs
269 ! ****************
270 ! coherent scattering
271  implicit none
272 #include "Zcode.h"
273 #include "Ztrack.h"
274 #include "Ztrackv.h"
275 ! ******************
276 ! since coherent scattering is effective at
277 ! low energies where angular distribution can be
278 ! approximated by (1+cos^2) dcos, we simply use
279 
280  type(track)::aTrack
281  type(coord)::w
282  real*8 cosg, tmp, cs, sn, sing
283 ! sample scattering angle from (1+cos^2)dcos
284  call ksamprsa(cosg)
285  tmp=1.d0-cosg*cosg
286  sing = sqrt(tmp)
287  call kcossn(cs,sn)
288 
289  w%r(1) = cs*sing
290  w%r(2) = sn*sing
291  w%r(2) = cosg
292  call ctransvectz(movedtrack%vec%w, w, w)
293  atrack = movedtrack
294 ! energy unchaged;
295  call csetdircos(w, atrack)
296  call ce2p(atrack)
297  nproduced = nproduced + 1
298  pwork(nproduced) = atrack%p
299  end
300 ! ****************
301  subroutine cphotoee
302 ! ****************
303 ! photo electric effect
304  implicit none
305 #include "Zcode.h"
306 #include "Zmass.h"
307 #include "Ztrack.h"
308 #include "Ztrackv.h"
309 
310 ! ******************
311  type(track)::aTrack
312  type(coord)::w
313  real*8 cs, sn, sing
314 
315  real*8 eout, eg, rEg, rEe, cost, a, tmp
316 ! essentially no shell down to 1keV.
317  eg = movedtrack%p%fm%p(4)
318  eout = masele + eg
319 
320  atrack = movedtrack
321  reg=eg/masele
322  ree=(eout-masele) /masele
323  if(ree .le. 0.) then
324  cost=1.
325  else
326  a = ( targetatomicn/137.0/137.0 + 2.0*ree + reg**2 )/
327  * (2.0*reg*sqrt(2.0*ree) )
328  call ksamppeang(a, cost)
329  endif
330 
331  tmp=1.d0-cost*cost
332  sing = sqrt(tmp)
333  call kcossn(cs,sn)
334 
335  w%r(1) = cs*sing
336  w%r(2) = sn*sing
337  w%r(2) = cost
338  call ctransvectz(movedtrack%vec%w, w, w)
339  atrack%p%fm%p(4) = eout
340  call cmkptc(kelec, 0, -1, atrack%p)
341  call csetdircos(w, atrack)
342  call ce2p(atrack)
343  nproduced = nproduced + 1
344  pwork(nproduced) = atrack%p
345  end
346 ! *******************
347  subroutine cmpair
348 ! magnetic pair creation
349 ! *******************
350  implicit none
351 #include "Zcode.h"
352 #include "Ztrack.h"
353 #include "Ztrackv.h"
354 !
355  real*8 e1, e2, u
356  integer ica, nc
357  type(track)::aTrack
358 !
359  call cmpaire(xai, e2, nc)
360 ! e2 is higher energy fraction; change to real energy
361 ! store later in working array, then higher one is
362 ! stored first in the stack to save the memory.
363  e2 = movedtrack%p%fm%p(4) * e2
364  e1=movedtrack%p%fm%p(4) - e2
365 ! assign charge for e1
366  call rndc(u)
367  if(u .gt. .5) then
368  ica=1
369  else
370  ica=-1
371  endif
372 !
373  atrack = movedtrack
374  atrack%p%fm%p(4) = e1
375  call ce2p(atrack)
376  call cmkptc(kelec, 0, ica, atrack%p)
377  nproduced = nproduced + 1
378  pwork(nproduced) = atrack%p
379 
380  atrack%p%fm%p(4) = e2
381  call ce2p(atrack)
382  call cmkptc(kelec, 0, -ica, atrack%p)
383  nproduced = nproduced + 1
384  pwork(nproduced) = atrack%p
385  end
386 
387 !/////////////
388  subroutine cpredpmjet(pj, a, ntp)
389 #include "Zair.h"
390 #include "Zptcl.h"
391  type(ptcl)::pj
392  integer ntp
393  type(ptcl)::a(*)
394  targetnucleonno=14
395  targetprotonno= 7
396  call cdpmjet( pj, targetnucleonno, targetprotonno,
397  * a, ntp)
398  write(0,*) ' ntp =',ntp
399  stop
400  end
401 !/////////////
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
subroutine cphotop
Definition: cintePhoton.f:217
subroutine cgphad(massN, atomicN, pj, a, ntp)
Definition: cgpHad.f:75
subroutine cpairang(e, m, teta)
Definition: cPBA.f:34
! for length to thickness conversion or v v ! integer maxnodes real Hinf ! This is used when making table for dim simulation ! The slant length for vertical height to km is ! divided by LenStep m steps ! It can cover the slant length of about km for cos
Definition: Zatmos.h:8
subroutine cpairerglpm(eg, rho, ee)
Definition: cpairErgLPM.f:20
subroutine cpairenergy(Eg, Ee)
Definition: cpairPath.f:37
subroutine ce2p(aTrack)
Definition: ce2p.f:5
const int kphoton
Definition: Zcode.h:6
Definition: Ztrack.h:44
subroutine cmpair
Definition: cintePhoton.f:348
max ptcl codes in the kelec
Definition: Zcode.h:2
subroutine ccohs
Definition: cintePhoton.f:269
subroutine cmpaire(xai, ee, nc)
Definition: cmPairE.f:14
subroutine rndc(u)
Definition: rnd.f:91
subroutine cphotoee
Definition: cintePhoton.f:302
max ptcl codes in the kseethru ! subcode integer kcasg
Definition: Zcode.h:2
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
subroutine ctransvectz(zax, dir1, dir2)
Definition: ctransVectZ.f:21
subroutine csetdircos(dc, aTrack)
Definition: cgetZenith.f:4
subroutine ksamprsa(costheta)
Definition: ksampRSA.f:8
subroutine kcossn(cs, sn)
Definition: kcossn.f:13
const double masele
Definition: Zmass.h:2
subroutine cmkptc(code, subcode, charge, p)
Definition: cmkptc.f:15
Definition: Zptcl.h:75
Definition: Zcoord.h:43
subroutine cintephoton
Definition: cintePhoton.f:2
subroutine cpredpmjet(pj, a, ntp)
Definition: cintePhoton.f:389
subroutine ccomptea(Eg, Egout, Ee, cosg, cose)
Definition: ccomptPath.f:50
subroutine cpair
Definition: cintePhoton.f:50
subroutine ksamppeang(ain, cost)
Definition: ksampPEang.f:16
subroutine ccompt
Definition: cintePhoton.f:166