*------- このプログラムは、mu中間子が空から降ってくる量を計算するものである。
double precision xmuon(41), xwd(41), xnum(10000)
double precision s, yr, time, hei , xp , torr ,xn ,thres, n, k
double precision fac, xl, phm, tri, prob, pro
open (10, file="cer2res.dat", status="unknown", form="formatted")
data xmuon / 7.7e-09, 1.7e-08, 2.6e-08, 3.3e-08, 4.4e-08,
& 5.3e-08, 6.3e-08, 8.4e-08, 9.0e-08,
& 1.0e-07, 2.1e-07, 3.4e-07, 4.6e-07, 5.2e-07,
& 6.7e-07, 8.1e-07, 9.2e-07, 1.1e-06,
& 1.2e-06, 2.6e-06, 3.6e-06, 4.8e-06, 5.8e-06,
& 5.3e-06, 4.8e-06, 4.0e-06, 3.2e-06,
& 2.3e-06, 3.7e-07, 1.1e-07, 4.0e-08, 1.4e-08,
& 5.5e-09, 2.6e-09, 1.3e-09, 6.2e-10,
& 2.2e-10, 2.0e-11, 9.0e-12, 3.0e-12, 6.0e-13 /
data xwd / 1e-03 , 4e-03 , 9e-03 , 1.6e-02, 2.5e-02,
& 3.6e-02, 4.9e-02, 6.4e-02, 8.1e-02,
& 0.1 , 0.4 , 0.9 , 1.6 , 2.5 ,
& 3.6 , 4.9 , 6.4 , 8.1 ,
& 10 , 40 , 90 , 160 , 250 ,
& 360 , 490 , 640 , 810 ,
& 1000 , 4000 , 9000 , 16000 , 25000 ,
& 36000 , 49000 , 64000 , 81000 ,
& 1e5 , 2e5 , 3e5 , 5e5 , 1e6 /
print *, 'detectorの面積を入力しなさい。(cm^2)'
read *, s
print *, '観察期間を入力しなさい。(年)'
read *, yr
time = yr*31556926
print *, '2つのdetectorの間の距離を入力しなさい。(cm)'
read *, hei
print *, 'detectorの厚さを入力しなさい。(cm)'
read *, xl
print *, '分割幅を入力しなさい。'
read *, n
print *, 'photonのトリガーの数をいくつにするか入力せよ。'
read *, tri
jt = int(tri)
n1 = int(n)
fac = s*s*time*4/(hei*hei)
*------ remark! この実験では集光率を10%、量子効率を20%としている。
do 10 i=1, n1+1
xp = -10*(i-1)/n+5
torr=exp( xp*2.3025851)
xn = 3.52e-7 * torr
phm = 980* xl* xn
thres = 105.65 / sqrt(2*xn)
*------- remark ! グラフのエネルギーは運動エネルギーを表している。

prob=1
do 11 k=0, jt
call poisson (k, phm*0.02, pro)
prob=prob-pro
11 continue
do 20 j=1, 40
if (xwd(j).ge.thres) then
xnum(i)=xnum(i)+xmuon(j+1)*(xwd(j+1)-xwd(j))*prob*prob*fac
*------- remark! 2つのdetectorの効果を考えたので確率は(prob)^2で評価される！
endif
20 continue
print *, xp, xnum(i)
write (10,40) xp , xnum(i)
40 format(1x, 2f27.14)
10 continue
end
subroutine poisson (n, ave, res)
double precision ave, res, xxx, n
res = -ave + n* log(ave)
do 50 xxx=1, n
res = res - log(xxx)
50 continue
res = exp(res)
end