COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cinipipx.f
Go to the documentation of this file.
1  subroutine cinipipx
2  implicit none
3 ! initialize leading particle x sampling for pi+p-->pi leading
4 !
5 #include "Zcinippxc.h"
6 
7  external cintdndx2
8 
9  integer i, icon
10  logical first/.true./
11  save first
12 
13  real*8 s
14  real*8 x, x1, eps, ans
15  real*8 xa(npitbl)
16  real*8 dndx(npitbl) ! unnormalized dn/dx of leading particle at
17  data xa/0,
18  * .01, .02, .03, .04, .05, .06, .07, .08, .09, .1,
19  * .11, .15, .20, .25, .30, .35, .40, .45, .50, .55,
20  * .60, .65, .70, .725, .75, .77, .78, .79, .80,.82,
21  * .84, .85, .87, .88, .90, .91, .92, .93, .94, .95,
22  * .96, .97, .98, .99, 1.00/
23  data dndx/0,
24  * 31., 38., 43., 48., 52., 53.5, 55.5, 58., 59., 60.,
25  * 58.5, 58., 52, 44, 36, 30, 25, 21, 17, 14.5,
26  * 12.5, 11, 9.6, 9.0, 8.8, 8.6, 8.5, 8.4, 8.4, 8.4,
27  * 8.6, 8.8, 9.2, 9.5, 10.1, 11.5, 12.2, 14.2, 16., 21.,
28  * 26.0, 33., 50, 80, 130./
29 
30  eps = 1.d-4
31  if(first) then
32 ! move xa into xval;
33  do i = 1, npitbl
34  xval(i) = xa(i)
35  enddo
36  first = .false.
37  endif
38 
39 ! integral of dn/dx from 0 to 1.
40  call ktrpzintt2(dndx, 1, npitbl, xa, 1, 0.d0, 1.d0, s)
41 !
42 ! make normalized table
43 !
44  do i = 1, npitbl
45  ndndx(i) = dndx(i)/s
46  enddo
47 ! make table of intgral(0:x) of ndndx. for x = 0 to 1.0 step
48 ! 0.01
49  x = 0.
50  dx = 0.01d0
51  do i=2, n-1
52  x = x + dx
53  call ktrpzintt2(ndndx, 1, npitbl, xa, 1, 0.d0,
54  * x, intendndx2(i) )
55  enddo
56  intendndx2(1) = 0.
57  intendndx2(n) = 1.
58 ! solve inte(0:x) of ndndx = u for u = 0 to 1.0 step 0.01
59  x = dx
60  uconst = 0.
61  x1 = 0.
62  do i = 2, n-1
63  uconst = uconst + dx
64  call kbinchop(cintdndx2, x1, 1.d0, x, eps, ans, icon)
65  if(icon .ne. 0) then
66  call cerrormsg(
67  * 'failed in making pi-leading particle sampling table', 0)
68  endif
69  pipsx(i) = ans
70  x = ans
71  x1 = x
72  enddo
73  pipsx(1) = 0.
74  pipsx(n) = 1.
75  end
76  real*8 function cintdndx2(x)
77  implicit none
78  real*8 x
79 
80 #include "Zcinippxc.h"
81 
82  real*8 ans
83 
84  call ktrpzintt2(ndndx, 1, npitbl, xval, 1, 0.d0, x, ans)
85  cintdndx2 = ans - uconst
86  end
87 
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
cintdndx2
real *8 function cintdndx2(x)
Definition: cinipipx.f:77
dx
integer npitbl real *nx dx real dx
Definition: Zcinippxc.h:10
ndndx
integer npitbl real *nx dx real ndndx
Definition: Zcinippxc.h:10
xval
integer npitbl real *nx dx real xval
Definition: Zcinippxc.h:10
kbinchop
subroutine kbinchop(f, x1, x2, x, eps, ans, icon)
Definition: kbinChop.f:20
uconst
integer npitbl real *nx dx real uconst
Definition: Zcinippxc.h:10
ktrpzintt2
subroutine ktrpzintt2(t, intv, n, xt, intvx, a, b, ans)
Definition: ktrpzIntT2.f:5
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
intendndx2
integer npitbl real *nx dx real * intendndx2
Definition: Zcinippxc.h:10
pipsx
integer npitbl real *nx dx real pipsx
Definition: Zcinippxc.h:10
false
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
cinipipx
subroutine cinipipx
Definition: cinipipx.f:2
n
integer n
Definition: Zcinippxc.h:1