COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cinelx.f
Go to the documentation of this file.
1 ! cinelx: x A / AA' inelastic xsection
2 !
3  subroutine cinelx(pj, A, Z, xs)
4  implicit none
5 #include "Zptcl.h"
6 #include "Zcode.h"
7 #include "Zmass.h"
8 #include "Zevhnp.h"
9  type(ptcl)::pj !input projectile particle
10  real(8),intent(in):: A ! input. effective target mass number
11  real(8),intent(in):: Z ! input. effective target charge number
12  real(8),intent(out):: xs ! output. in mb. inelastic xs.
13 
14  real(8):: gamma
15  type(ptcl):: target
16  integer:: tcharge
17  real(8):: pjA, pjZ
18 
19 
20  if( a == 1.0 .and. pj%code == kgnuc ) then ! target is p/n, pj is heavy ion
21  ! we use proj. rest system ; if not, inside cinelx0
22  ! cAAXsec2 is called and it calles cinelx and recursive
23  ! call happens.
24  tcharge = z
25  call cmkptc(knuc, -1, tcharge, target) ! p or n
26  gamma = pj%fm%p(4)/pj%mass
27  pja = pj%subcode
28  pjz = pj%charge
29  if( tcharge == 0 ) then
30  target%fm%p(4) = gamma*masn
31  else
32  target%fm%p(4) = gamma*masp
33  endif
34  call cinelx0(target, pja, pjz, xs) ! target is now proj. pjA is target
35 
36  else
37  call cinelx0(pj, a, z, xs)
38  endif
39  end
40 
41  subroutine cinelx0(pj, A, Z, xs)
42  use modcosxsec ! refer only nIsDiff
43  implicit none
44 #include "Zptcl.h"
45 #include "Zcode.h"
46 #include "Zevhnp.h"
47  type(ptcl)::pj !input projectile particle
48  real(8),intent(in):: A ! input. effective target mass number
49  real(8),intent(in):: Z ! input. effective target charge number
50  real(8),intent(out):: xs ! output. in mb. total xsection.
51 
52 
53  real(8) p
54  real(8),external:: cPDGsigmaInepA
55  real(8),save:: Aprev=-100., norm
56  real(8):: xs200, xs200A
57 
58 
59  real(8)::shp, shn
60 
61  if( pj%fm%p(4) .le. pj%mass) then
62  p =0.
63  if(pj%code .eq. knuc .and. pj%subcode .eq. antip) then
64  xs = largexs
65  elseif(pj%code .eq. kelec .and. pj%charge .eq. 1) then
66  xs = largexs
67  else
68  xs = smallxs
69  endif
70  else
71  p = sqrt(pj%fm%p(4)**2 - pj%mass**2)
72  if(pj%code .eq. knuc) then
73  if(pj%charge .eq. 1) then
74  if( a == 1.0 ) then
75  if( z == 1.0 ) then
76  call cppinelaxs(p, xs)
77  elseif( z == 0.) then
78  call cnpinelaxs(p, xs)
79  else
80  write(0,*) 'A,Z=',a,z, 'for p pj'
81  write(0,*) 'in cinelx; invalid'
82  stop
83  endif
84  elseif(a > 1.) then
85  call cppinelaxs(p, shp)
86  if( nisdiff ) then
87  call cnpinelaxs(p, shn)
88  shp = (shp + shn)/2.0
89  endif
90  call cxp2xaxsec(a, shp, xs)
91  else
92  write(0,*) ' A,Z for cinelx invalid'
93  write(0,*) ' pj is p '
94  stop
95  endif
96  elseif(pj%charge .eq. -1) then
97  call cpbarpinelaxs(p, shp)
98  if( a .gt. 1.) then
99  if(nisdiff) then
100  call cnbarpinelaxs(p, shn)
101  shp = (shp + shn)/2.
102  endif
103  call cxp2xaxsec(a, shp, xs)
104  else
105  xs = shp
106  endif
107  elseif(pj%subcode .eq. antip) then
108 ! anti-neutron; assume the same one as pbar
109  call cnbarpinelaxs(p, shp)
110  if( a .gt. 1.) then
111  call cxp2xaxsec(a, shp, xs)
112  else
113  xs = shp
114  endif
115  else
116 ! neutron
117  call cnpinelaxs(p, shp)
118  if(a .gt. 1.) then
119  call cxp2xaxsec(a, shp, xs)
120  else
121  xs = shp
122  endif
123  endif
124  elseif(pj%code == kpion ) then
125  if(pj%charge .eq. 1) then
126  call cpippinelaxs(p, shp)
127  elseif(pj%charge .eq. -1) then
128  call cpimpinelaxs(p, shp)
129  else
130  ! pi0
131  call cpimpinelaxs(p, shp)
132  endif
133  if(a .ne. 1.0) then
134  call cxp2xaxsec(a, shp, xs)
135  else
136  xs = shp
137  endif
138  elseif(pj%code .eq. kkaon) then
139  if(pj%charge .eq. 1) then
140  if(a == 1.0) then
141  if( z==1.0) then
142  call ckppinelaxs(p, xs)
143  else
144  call ckpninelaxs(p, xs)
145  endif
146  elseif(a .gt. 1.) then
147  call ckppinelaxs(p, shp)
148  if( nisdiff ) then
149  call ckpninelaxs(p, shn)
150  shp = (shp + shn)/2.0
151  endif
152  call cxp2xaxsec(a, shp, xs)
153  else
154  write(0,*) 'A,Z=',a,z, ' invalid for cinel%Kaon'
155  stop
156  endif
157  elseif(pj%charge .eq. -1) then
158  call ckmpinelaxs(p, shp)
159  if(a .gt. 1.) then
160  if( nisdiff ) then
161  call ckmninelaxs(p, shn)
162  shp = (shp + shn)/2.
163  endif
164  call cxp2xaxsec(a, shp, xs)
165  else
166  xs = shp
167  endif
168  else
169 ! k0; don't worry so much
170  call ckmpinelaxs(p, shp)
171 !! call ckmnInelaXs(p, shn)
172 !! shp = (shp+ shn)/2.
173  if(a .gt. 1.) then
174  call cxp2xaxsec(a, shp, xs)
175  else
176  xs = shp
177  endif
178  endif
179  elseif(pj%code .eq. kgnuc ) then
180 ! heavy ion
181  call caaxsec2(pj, a, z, xs)
182  elseif(pj%code .eq. kdmes) then
183  call ckppinelaxs(p, shp)
184  call cxp2xaxsec(a, shp, xs)
185  elseif(pj%code .eq. kgzai .or. pj%code .eq. ksigma .or.
186  * pj%code .eq. kbomega .or. pj%code .eq. klambda .or.
187  * pj%code .eq. klambdac) then
188 ! don't worry, not used almost at all; use proton
189 ! shp = ctotpp1(p)
190  call cppinelaxs(p, shp)
191  call cxp2xaxsec(a, shp, xs)
192  elseif(pj%code == kphoton ) then
193  call cgpxsec(a, pj%fm%p(4), xs)
194  elseif(pj%code .eq. kneumu) then
195  xs= smallxs
196  elseif(pj%code .eq. kneue) then
197  xs= smallxs
198  elseif(pj%code == kneutau) then
199  xs= smallxs
200  elseif(pj%code == knnb .or. pj%code == kddb ) then
201  xs = smallxs
202  elseif(pj%code == kmuon) then
203  write(0,*) 'cinelx: code=',pj%code, ' should not come'
204  stop
205  else
206 ! eta, ds, Xic, etc
207  xs = smallxs
208 !! use pion
209  ! call cpippInelaXs(p, shp)
210  ! call cxp2xAXsec(A, shp, xs)
211  endif
212  endif
213 !!! normalize to pdg at 200 GeV. norm const is
214 !!! applied to other collisions
215  if( xs /= smallxs ) then
216  if(a > 1. .and. a /= aprev ) then
217  call cppinelaxs(200.d0, xs200)
218  call cxp2xaxsec(a, xs200, xs200a)
219 !! write(0,*) ' xs200A=',xs200A, ' pgd=', cPDGsigmaInepA(A)
220  norm = cpdgsigmainepa(a) /xs200a
221 !! write(0,*) ' norm=', norm
222  aprev = a
223  elseif( a == 1.) then
224  norm = 1.
225  aprev = a
226  endif
227  xs = xs*norm
228  endif
229  end
230  function cinelcosbypdg(A) result(ratio)
231 ! sigmaInela(pA)cosmos/sigmaInela(pA)PDG
232 ! @ 200 GeV. SigmaInela by cosmos is the one
233 ! obtained by cinelx0 before
234 ! cxp2xAXsec included cxAbyxpXsec.
235 !
236  real(8),intent(in):: A
237  real(8):: ratio
238 
239  real(8):: temp
240 
241  if( a > 48. ) then
242  if(a >150.) then
243  ratio = 1.0
244  else
245  ratio = (1.006+0.0268*a/100.)*0.957
246  endif
247  elseif(a > 20.) then
248  ratio = 0.975
249  elseif( abs(a - 16.) < 0.1 ) then
250  ratio = 0.9427
251  elseif( abs( a- 7.0) < 0.1 ) then
252  ratio = 0.902
253  elseif( a > 7. ) then
254  temp = log(a)
255  ratio = (-0.100*temp + 0.57311)*temp +0.1578
256  elseif( abs( a- 4.) < 0.1 ) then
257  ratio = 0.975
258  elseif( abs( a- 2.) < 0.1 ) then
259  ratio = 1.137
260  elseif (abs(a-1.0) < 0.1 ) then
261  ratio =1.0138
262  else
263  ratio = -0.025*(a-4) + 0.975
264  endif
265  end
266 
267 
norm
common ZdedxAir norm
Definition: ZdedxAir.h:2
kgzai
max ptcl codes in the kgzai
Definition: Zcode.h:2
cinelx
subroutine cinelx(pj, A, Z, xs)
Definition: cinelx.f:4
kdmes
max ptcl codes in the kdmes
Definition: Zcode.h:2
kgnuc
max ptcl codes in the kgnuc
Definition: Zcode.h:2
masn
masn
Definition: Zmass.h:5
cnpinelaxs
subroutine cnpinelaxs(p, xs)
Definition: cnpTotXs.f:142
kphoton
const int kphoton
Definition: Zcode.h:6
klambdac
max ptcl codes in the klambdac
Definition: Zcode.h:2
kkaon
max ptcl codes in the kkaon
Definition: Zcode.h:2
kelec
max ptcl codes in the kelec
Definition: Zcode.h:2
kneue
max ptcl codes in the kneue
Definition: Zcode.h:2
cinelx0
subroutine cinelx0(pj, A, Z, xs)
Definition: cinelx.f:42
cxp2xaxsec
subroutine cxp2xaxsec(a, xsxp, xsxa)
Definition: cxp2xAXsec.f:30
kneutau
max ptcl codes in the kneutau
Definition: Zcode.h:2
caaxsec2
subroutine caaxsec2(pj, tgA, tgZ, xs)
Definition: cAAXsec2.f:2
modcosxsec::nisdiff
logical, save nisdiff
Definition: ctotx.f:4
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
klambda
max ptcl codes in the klambda
Definition: Zcode.h:2
knuc
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code knuc
Definition: cblkHeavy.h:7
kneumu
max ptcl codes in the kneumu
Definition: Zcode.h:2
cpimpinelaxs
subroutine cpimpinelaxs(p, xs)
Definition: cpimpTotXs.f:145
ckmpinelaxs
subroutine ckmpinelaxs(p, xs)
Definition: ckmpTotXs.f:130
ckpninelaxs
subroutine ckpninelaxs(p, xs)
Definition: ckpnTotXs.f:112
knnb
max ptcl codes in the knnb
Definition: Zcode.h:2
cinelcosbypdg
real(8) function cinelcosbypdg(A)
Definition: cinelx.f:231
cmkptc
subroutine cmkptc(code, subcode, charge, p)
Definition: cmkptc.f:15
cnbarpinelaxs
subroutine cnbarpinelaxs(p, xs)
Definition: cnbarpTotXs.f:16
masp
masp
Definition: Zmass.h:5
ptcl
Definition: Zptcl.h:75
antip
max ptcl codes in the kseethru ! subcode integer antip
Definition: Zcode.h:2
ckppinelaxs
subroutine ckppinelaxs(p, xs)
Definition: ckppTotXs.f:116
kpion
max ptcl codes in the kpion
Definition: Zcode.h:2
cpbarpinelaxs
subroutine cpbarpinelaxs(p, xs)
Definition: cpbarpTotXs.f:116
ksigma
max ptcl codes in the ksigma
Definition: Zcode.h:2
kddb
max ptcl codes in the kddb
Definition: Zcode.h:2
kmuon
max ptcl codes in the kmuon
Definition: Zcode.h:2
cppinelaxs
subroutine cppinelaxs(p, xs)
Definition: cppTotXs.f:124
kbomega
max ptcl codes in the kbomega
Definition: Zcode.h:2
ckmninelaxs
subroutine ckmninelaxs(p, xs)
Definition: ckmnTotXs.f:14
cgpxsec
subroutine cgpxsec(a, energy, xs)
Definition: cgpXsec.f:11
cpippinelaxs
subroutine cpippinelaxs(p, xs)
Definition: cpippTotXs.f:180