7
def frange(startValue, endValue, steps):
8
incrementValue = (float(endValue) - startValue) / steps
10
yield float(startValue) + i*incrementValue
14
parent_energy = 3.0 # GeV
17
E_max = parent_energy # GeV
20
mass_muon = 0.1056 # GeV
26
distance = 1 # km, GLoBES normalization
28
flux_file_plus = open('flux_file_plus.dat', 'w')
29
flux_file_minus = open('flux_file_minus.dat', 'w')
31
for E in frange (E_min, E_max, E_steps):
33
y = float(E) / parent_energy
34
beta = math.sqrt(1 - mass_muon**2 / parent_energy**2)
37
N /= pi * distance**2 * mass_muon**6
38
N *= parent_energy**4 * y**2 * (1 - beta * math.cos(theta))
40
flux_mu = 4 * N * (3 * mass_muon**2 - 4 * parent_energy**2 * y * (1 - beta * math.cos(theta)))
41
flux_e = 24 * N * (1 * mass_muon**2 - 2 * parent_energy**2 * y * (1 - beta * math.cos(theta)))
44
flux_e *= (2*pi) / parent_energy # solid angle -> cos(theta), y -> E
45
flux_mu *= (2*pi) / parent_energy # solid angle -> cos(theta), y -> E
47
flux_file_plus.write('%0.9g %0.9g %0.9g 0 %0.9g %0.9g 0\n' % (E, flux_e, 0, 0, flux_mu))
48
flux_file_minus.write('%0.9g %0.9g %0.9g 0 %0.9g %0.9g 0\n' % (E, 0, flux_mu, flux_e, 0))
51
flux_file_plus.close()
52
flux_file_minus.close()