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

Go to the source code of this file.

Functions/Subroutines

subroutine chookeabsorbi (info)
 
subroutine chookeabsorb (a, b, dE, info)
 This is called when Eabsorb != 0 and when a charged particle runs from a to b and deposits energy dE (GeV) to the Air. More...
 
subroutine chookeabsorbd (a, dE, info)
 
subroutine chookeabsorbb (a, info)
 
subroutine chookeabsorbc (a, n, p, info)
 
subroutine chookeabsorbw (a, n, p, Eout, diff)
 

Function/Subroutine Documentation

◆ chookeabsorb()

subroutine chookeabsorb ( type(track a,
type(track b,
real*8  dE,
integer  info 
)

This is called when Eabsorb != 0 and when a charged particle runs from a to b and deposits energy dE (GeV) to the Air.

Parameters
a
b
dE
info
Todo:
test 2016/2/16

Definition at line 60 of file chookEabsorb.f.

Referenced by ctracking().

60  implicit none
61 #include "Zmaxdef.h"
62 #include "Zobs.h"
63 #include "Ztrack.h"
64 #include "Zabsorb.h"
65 ! This is called when Eabsorb != 0 and
66 ! when a charged particle runs from a
67 ! to b and deposits energy dE (GeV) to the Air.
68 !
69  type(track):: a ! input. charged particle track info. at a.
70  type(track):: b ! input. charged particle track info. at b.
71  real*8 de ! input. energy deposit in GeV; no weight is applied yet
72  integer info ! in/out. not used now
73 
74  if(a%where .ge. 1) then
75  debydedx(a%where)= debydedx(a%where) + de*a%wgt
76  endif
Definition: Ztrack.h:44
real(4), save a
Definition: cNRLAtmos.f:20
real(4), save b
Definition: cNRLAtmos.f:21
Here is the caller graph for this function:

◆ chookeabsorbb()

subroutine chookeabsorbb ( type(track a,
integer  info 
)

Definition at line 146 of file chookEabsorb.f.

References false, and true.

Referenced by cifdead(), and cobservation().

146  implicit none
147 #include "Zmaxdef.h"
148 #include "Zcode.h"
149 #include "Ztrack.h"
150 #include "Ztrackp.h"
151 #include "Ztrackv.h"
152 #include "Zobs.h"
153 #include "Zobsp.h"
154 #include "Zobsv.h"
155 #include "Zabsorb.h"
156 ! This is called when Eabsorb(1) !=0 and
157 ! a particle crosses an observation level
158 ! info=2: normal observation level
159 ! =1: upper boundary (whcih is equal to or higher than the injection
160 ! height)
161 ! =3: lower boundary (which is equal to or lower than the last
162 ! observation depth).
163 !
164  type(track):: a ! input. a particle that croses the observation level
165  integer info ! input. see above
166 
167  integer code
168  integer lv
169  logical count
170 
171  count =.false.
172  if(info .eq. 2 .and. a%where .eq. eabsorb(2) ) then
173 ! Eabsorb(2) is the layer number where the user want to
174 ! take sum of particle energy reaching there from above.
175 ! ptcl comes to the specified level OR
176 ! upper boundary. This condition neglect
177 ! ptcles reaching the real lower boundary
178 ! (if Eabosrb(2) is not NoOfSites+1)
179  if( a%vec%coszenith .gt. 0.) then
180 ! if info=2, this is to count energy reaching
181 ! the level from above, we discard going up ptcls
182 ! (but it may come down, so some over count may happen)
183  count = .true.
184  lv = 1
185  endif
186  elseif( info .eq. 1 ) then
187  lv = 2 ! escape to space
188  count = .true.
189  endif
190 !
191 !
192  if(count) then
193  code = a%p%code
194 ! for other ptcls than g,e,mu, pi, K, N, use 7
195  if(code .gt. 7) code=7
196 ! at the last layer we see sum of the particle energy
197 ! what energy we should use here is somewhat annoying
198 ! point. We use total energy here.
199  if(lv .eq. 1) then
200 ! Ecrash(code) = Ecrash(code) + a.p.fm.p(4)*a.wgt
201  ecrash(code) = ecrash(code) + (a%p%fm%p(4)-a%p%mass)*a%wgt
202  else
203 ! Espace(code) = Espace(code) + a.p.fm.p(4)*a.wgt
204  espace(code) = espace(code) + (a%p%fm%p(4)-a%p%mass)*a%wgt
205  endif
206  endif
Definition: Ztrack.h:44
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
*Zfirst p fm *Zfirst p Zfirst p code
Definition: ZavoidUnionMap.h:1
real(4), save a
Definition: cNRLAtmos.f:20
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
Here is the caller graph for this function:

◆ chookeabsorbc()

subroutine chookeabsorbc ( type(track a,
integer  n,
type(ptcl), dimension(n p,
integer  info 
)

Definition at line 209 of file chookEabsorb.f.

References chookeabsorbw(), masn, and masp.

Referenced by cinteraction().

209  use modxsecmedia
210  implicit none
211 #include "Zmaxdef.h"
212 #include "Zcode.h"
213 #include "Zmass.h"
214 #include "Zobs.h"
215 #include "Ztrack.h"
216 #include "Ztrackv.h"
217 #include "Zabsorb.h"
218 
219  type(track):: a ! input. incident particle track
220  ! that made the collision
221  integer n ! input. number of procuded ptcls in the collision
222  type(ptcl):: p(n) ! input. produced ptpcls.
223  integer info ! input. not used now.
224 
225  real*8 eout, diff, reldiff, ein
226  real*8 diff1, diff2, ein1, ein2
227  integer i
228 ! We check conservation above this energy (GeV).
229  real*8 ebig/5.d3/
230  save
231 ! since target is not well known and nucleon there
232 ! has Fermi momentum, we simply assume rest mass.
233 ! For E>5TeV, Ar mass ~40GeV/5000 < 4/500 < 1 %
234 ! So we can see the conservation neglecting mass term
235 !
236  if( a%p%fm%p(4) .gt. ebig ) then
237  ein1 = a%p%fm%p(4) + masn*(targetnucleonno-targetprotonno) +
238  * masp*targetprotonno
239  ein2 = a%p%fm%p(4) + masp
240  eout = 0.
241  do i = 1, n
242  eout = eout + p(i)%fm%p(4)
243  enddo
244 
245  diff1 = eout - ein1
246  diff2 = eout - ein2
247  ein = ein1
248  diff = diff1
249 ! take smaller diff one for Ein
250  if( abs(diff2) .lt. abs(diff1)) then
251  diff = diff2
252  ein = ein2
253  endif
254  reldiff = (eout/ein -1.)
255 
256  if( abs(reldiff) .gt. 0.02 ) then
257  call chookeabsorbw(a, n, p, eout, diff)
258  endif
259 
260  if( abs(diff) .gt. abs(maxebreak(1)) ) then
261  maxebreak(1) = diff
262  maxebreak(2) = reldiff
263  endif
264 
265  if( abs(reldiff) .gt. abs(maxrelebreak(1)) ) then
266  maxrelebreak(1) = reldiff
267  maxrelebreak(2) = diff
268  endif
269 
270  sumediff = sumediff + diff
271  sumabsediff = sumabsediff + abs(diff)
272  endif
subroutine chookeabsorbw(a, n, p, Eout, diff)
Definition: chookEabsorb.f:275
nodes i
masn
Definition: Zmass.h:5
block data cblkIncident data *Za1ry *HeightOfInj d3
Definition: cblkIncident.h:7
Definition: Ztrack.h:44
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
real(4), save a
Definition: cNRLAtmos.f:20
masp
Definition: Zmass.h:5
Definition: Zptcl.h:75
integer n
Definition: Zcinippxc.h:1
Here is the call graph for this function:
Here is the caller graph for this function:

◆ chookeabsorbd()

subroutine chookeabsorbd ( type(track a,
real*8  dE,
integer  info 
)

Definition at line 79 of file chookEabsorb.f.

References kelec, kkaon, kneue, kneumu, knuc, kphoton, and regptcl.

Referenced by cifdead().

79  implicit none
80 #include "Zmaxdef.h"
81 #include "Zcode.h"
82 #include "Zobs.h"
83 #include "Ztrack.h"
84 #include "Zabsorb.h"
85 ! This is called when Eabsorb != 0 and
86 ! when a particle energy becomes < Emin (info=0) for
87 ! its traveling time becomes > limit or (info=2)
88 ! its angle relative to the 1ry becomes > limit.(info=4)
89 ! Whether this is called or not depends on the particle
90 ! and bit in Eabsorb. dE is energy that can be regarded
91 ! as absorbed in the Air. (GeV) eventually.
92 ! bit 1 is the LSB of Eabsorb.
93 !
94 ! bit particle
95 ! 1 photon: used to absorb shell energy at
96 ! photoelectric effect. This bit is not used in the Air.
97 ! 2 photon.
98 ! 3 e+/e-
99 ! 4 proton
100 ! 5 neutron
101 ! 6 anti-N
102 ! 7 decaying prtcl
103 ! 8 others
104 !
105 !***** Normally Eabsorb=6 (110 in bit pattarn) is enough.****
106 !
107 !
108  type(track):: a ! input. a particle that is < Emin
109  ! at birth
110  real*8 de ! input. energy which is supposed to be emitted by
111  ! the dying particle
112  integer info ! in/out. not used now.
113  if(a%where .ge. 1) then
114  if(a%p%code .eq. kneue .or. a%p%code .eq. kneumu) then
115  ! neutrino
116  debydeathneu(a%where) = debydeathneu(a%where) + de*a%wgt
117  elseif(a%p%code .eq. knuc .and. a%p%charge .eq. 0 .and.
118  * a%p%subcode .eq. regptcl ) then
119  ! low E nutron
120  debydeathnut(a%where) = debydeathnut(a%where) + de*a%wgt
121  else
122 ! next one and above are kept same as older versions for
123 ! compativilty
124  debydeath(a%where) = debydeath(a%where) + de*a%wgt
125 
126 ! we further put details
127  if(a%p%code .eq. kphoton ) then
128  debydeathg(a%where) = debydeathg(a%where) + de*a%wgt
129  elseif(a%p%code .eq. kelec ) then
130  debydeathe(a%where) = debydeathe(a%where) + de*a%wgt
131  elseif( a%p%code .le. kkaon ) then
132  debydeathmupik(a%where) = debydeathmupik(a%where)
133  * + de*a%wgt
134  elseif(a%p%code .eq. knuc .and. a%p%charge .eq. 1 ) then
135 ! p
136  debydeathp(a%where) = debydeathp(a%where) + de*a%wgt
137  else
138 ! pbar, nbar, heavy, others
139  debydeatho(a%where) = debydeatho(a%where) + de*a%wgt
140  endif
141  endif
142  endif
const int kphoton
Definition: Zcode.h:6
Definition: Ztrack.h:44
max ptcl codes in the kkaon
Definition: Zcode.h:2
max ptcl codes in the kelec
Definition: Zcode.h:2
max ptcl codes in the kneue
Definition: Zcode.h:2
max ptcl codes in the kseethru ! subcode integer regptcl
Definition: Zcode.h:2
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code knuc
Definition: cblkHeavy.h:7
max ptcl codes in the kneumu
Definition: Zcode.h:2
real(4), save a
Definition: cNRLAtmos.f:20
Here is the caller graph for this function:

◆ chookeabsorbi()

subroutine chookeabsorbi ( integer  info)

Definition at line 2 of file chookEabsorb.f.

Referenced by ceventloop().

2 ! init for each event
3  implicit none
4 #include "Zmaxdef.h"
5 #include "Zcode.h"
6 #include "Ztrack.h"
7 #include "Ztrackv.h"
8 #include "Zobs.h"
9 #include "Zobsp.h"
10 #include "Zobsv.h"
11 #include "Zabsorb.h"
12  integer info ! not used now
13  integer i
14  do i = 0, noofsites + 1
15 ! i = 0 and NoOfSite+1 will not be used now.
16  debydedx(i) = 0.
17  debydeath(i) = 0.
18  debydeathg(i) = 0.
19  debydeathe(i) = 0.
20  debydeathmupik(i) = 0.
21  debydeathneu(i) = 0.
22  debydeathp(i) = 0.
23  debydeathnut(i) = 0.
24  debydeatho(i) = 0.
25  enddo
26  do i = 1, 7
27 ! i for g,e,mu,pi,K, n, (neu and others)
28 ! Energy crashing to the given layer
29 ! default is the NoOfSites. Can be changed by Eabosrb(2)
30  ecrash(i) = 0.
31 ! Energy escaping to the upper bound( default is injection
32 ! height + 1m)
33  espace(i) = 0.
34  enddo
35 ! counter used for checking energy conservation
36 ! at the multiple producion.
37 ! negative value is generated energy is < ininial energy
38 ! positive value is generated energy is > initial energy
39  do i = 1, 2
40 ! i=1 for max Ebreak, i=2 Relative Break value for that event
41  maxebreak(i) = 0.
42 ! i=1 for max RelEgreak, i=2 E Break value for that event
43  maxrelebreak(i) = 0.
44  enddo
45 
46  sumediff = 0.
47  sumabsediff = 0.
48 
nodes i
Here is the caller graph for this function:

◆ chookeabsorbw()

subroutine chookeabsorbw ( type(track a,
integer  n,
type(ptcl), dimension(n p,
real*8  Eout,
real*8  diff 
)

Definition at line 275 of file chookEabsorb.f.

References cqeventno().

Referenced by chookeabsorbc().

275  use modxsecmedia
276  implicit none
277 #include "Zmaxdef.h"
278 #include "Zcode.h"
279 #include "Zmass.h"
280 #include "Zobs.h"
281 #include "Ztrack.h"
282 #include "Ztrackv.h"
283 #include "Zabsorb.h"
284 
285  type(track):: a ! input. incident particle track
286  ! that made the collision
287  integer n ! input. number of procuded ptcls in the collision
288  type(ptcl):: p(n) ! input. produced ptpcls.
289  real*8 eout, diff
290  integer nevent, ntevent
291 
292 
293  call cqeventno(nevent, ntevent)
294 
295  write(0,*) "================= Warning: Ebreak ============="
296  write(0,*) " event #=", nevent, " 1ry Energy=", inci%p%fm%p(4)
297  write(0,*) "incident(not 1ry) code=", a%p%code,
298  * " subcode=", a%p%subcode, " Ein~", a%p%fm%p(4)
299  write(0,*) "Target(A,Z)=", targetnucleonno, targetprotonno
300  write(0,*) "No. of generted particles =", n
301  write(0,*) "Sum of outgoing ptcl Energy:Eout=", eout
302  write(0,*) "dE=(Eout-Ein)=", diff, " dE/Ein~",
303  * diff/a%p%fm%p(4)
304  write(0,*) "(Eout-Ein)/1ryE=", diff/inci%p%fm%p(4)
305  write(0,*) "Height=", a%pos%height," depth=",a%pos%depth/10.
306  write(0,*) " where=", a%where, " weight=", a%wgt
307  write(0,*) "================================================"
308 
Definition: Ztrack.h:44
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
subroutine cqeventno(num, cumnum)
Definition: cqEventNo.f:3
real(4), save a
Definition: cNRLAtmos.f:20
Definition: Zptcl.h:75
integer n
Definition: Zcinippxc.h:1
Here is the call graph for this function:
Here is the caller graph for this function: