aboutsummaryrefslogtreecommitdiff
path: root/macros/snogen
blob: 31e9a742ada56bd6adda8a297181748e4a4f5fe9 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#!/usr/bin/env python
"""
snogen is a script to generate SNOMAN command files.

Example:

    # simulate 100 1 GeV muons
    $ snogen -n 100 -e 1000 -p mu_minus | /path/to/snoman.exe

The output filename can be specified on the command line:

    $ snogen -o [output filename]

By default it will be [particle name]_[energy]_[number of events].zdab.
"""
from __future__ import print_function, division
import string

class MyTemplate(string.Template):
    delimiter = '@'

if __name__ == '__main__':
    import argparse

    parser = argparse.ArgumentParser("generate a SNOMAN command file")
    parser.add_argument("-p", "--particle-id", default="e_minus",
                        help="particle id")
    parser.add_argument("-e", "--mc-energy", default=500.0, type=float,
                        help="total energy")
    parser.add_argument("-n", "--num-events", default=100, type=int,
                        help="number of events")
    parser.add_argument("-o", "--output", default=None,
                        help="output file name")
    args = parser.parse_args()

    if args.output is None:
        output = "%s_%.0f_%i.zdab" % (args.particle_id,args.mc_energy,args.num_events)
    else:
        output = args.output

    with open("mc_run_10000_macro.cmd") as f:
        s = f.read()

    template = MyTemplate(s)
    print(template.safe_substitute(particle_id=args.particle_id,mc_energy=args.mc_energy,num_events=args.num_events,output=output))
class="kn">import matplotlib.pyplot as plt fig = plt.figure() colors = plt.rcParams["axes.prop_cycle"].by_key()["color"] linestyles = ['-','--'] def key(filename): head, tail = os.path.split(filename) k = 0 if tail.startswith('fmax'): k += 1 if 'nue' in tail: k += 10 elif 'nbe' in tail: k += 20 elif 'num' in tail: k += 30 elif 'nbm' in tail: k += 40 elif 'nut' in tail: k += 50 elif 'nbt' in tail: k += 60 return k for filename in sorted(args.filenames,key=key): head, tail = os.path.split(filename) print(filename) data = np.genfromtxt(filename) shape1 = len(np.unique(data[:,0])) x = data[:,0].reshape((-1,shape1)) y = data[:,1].reshape((-1,shape1)) z = data[:,2].reshape((-1,shape1)) # Convert to MeV x *= 1000.0 z /= 1000.0 zbins = np.linspace(-1,1,21) dz = zbins[1] - zbins[0] x = x[0] # Integrate over cos(theta) and multiply by 2*pi to convert 3D flux to # a total flux y = np.sum(z*dz,axis=0)*2*np.pi if 'sno_nue' in tail: plt.plot(x,y,color=colors[0],linestyle=linestyles[0],label=r'$\nu_e$') elif 'sno_nbe' in tail: plt.plot(x,y,color=colors[0],linestyle=linestyles[1],label=r'$\overline{\nu}_e$') elif 'sno_num' in tail: plt.plot(x,y,color=colors[1],linestyle=linestyles[0],label=r'$\nu_\mu$') elif 'sno_nbm' in tail: plt.plot(x,y,color=colors[1],linestyle=linestyles[1],label=r'$\overline{\nu}_\mu$') elif 'sno_nut' in tail: plt.plot(x,y,color=colors[2],linestyle=linestyles[0],label=r'$\nu_\tau$') elif 'sno_nbt' in tail: plt.plot(x,y,color=colors[2],linestyle=linestyles[1],label=r'$\overline{\nu}_\tau$') plt.gca().set_xscale("log") plt.gca().set_yscale("log") despine(fig,trim=True) plt.xlabel("$E$ (MeV)") plt.ylabel(r"$\mathrm{d}\Phi/\mathrm{d}E$ (1/$\mathrm{m}^2$/sec/MeV)") plt.legend() plt.tight_layout() if args.save: plt.savefig("irc01_atmospheric_flux.pdf") plt.savefig("irc01_atmospheric_flux.eps") plt.show()