-
Notifications
You must be signed in to change notification settings - Fork 0
/
Energy.py
executable file
·164 lines (108 loc) · 5.47 KB
/
Energy.py
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
#!/bin/sh /cvmfs/icecube.opensciencegrid.org/py2-v2/icetray-start
#METAPROJECT icerec/V05-01-00
# This code generates histograms for the likelihood that describes the reconstructed signal neutrino energy
import numpy as np
import src/ING as ING #contains histograms
import h5py
import os
import argparse
import sys
#First get the signal spectra
parser = argparse.ArgumentParser(description='Process some parameters.')
parser.add_argument('-m', dest='massnumber', default=0, type=int, choices = range(13), metavar='massnumber', help='mass number')
args = parser.parse_args()
massnumber = args.massnumber
masslist = [100, 250, 350, 500, 750, 1000, 2500, 5000, 7500, 10000, 25000, 50000, 75000]
medmasslist = [1, 10, 100, 1000, 10000]
lifetimelist = [1, 10, 100, 1000, 10000]
channellist = [5, 8, 11, 13]
track = []
spectrum_mu = []
spectrum_amu = []
N_mu = []
N_amu = []
for mass in masslist:
trax = []
spect_mu = []
spect_amu = []
Nx_mu = []
Nx_amu = []
for medmass in medmasslist:
traxx = []
spec_mu = []
spec_amu = []
Nxx_mu = []
Nxx_amu = []
for lifetime in lifetimelist:
traxxx = []
spe_mu = []
spe_amu = []
Nxxx_mu = []
Nxxx_amu = []
for channel in channellist:
if os.path.isfile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_e.dat") and\
mass == masslist[massnumber]:
spe_mu.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_mu.dat"))
Nxxx_mu.append(spe_mu[-1].Integral(spe_mu[-1].points[0][0],spe_mu[-1].points[-1][0]))
spe_amu.append(ING.Graph.FromFile("/home/ctoennis/analyses/standard_analysis_framework/spectra_wimpsim/wa-sunsum-m" + str(mass) + "-ch" + str(channel) + "-x" + str(medmass) + "-l" + str(lifetime) + "-009080-000000-read-earth_amu.dat"))
Nxxx_amu.append(spe_amu[-1].Integral(spe_amu[-1].points[0][0],spe_amu[-1].points[-1][0]))
traxxx.append(ING.H1D.Empty(0.0,200000.0,2000))
else:
spe_mu.append("None")
Nxxx_mu.append("None")
spe_amu.append("None")
Nxxx_amu.append("None")
traxxx.append("None")
spec_mu.append(spe_mu)
spec_amu.append(spe_amu)
Nxx_mu.append(Nxxx_mu)
Nxx_amu.append(Nxxx_amu)
traxx.append(traxxx)
spect_mu.append(spec_mu)
spect_amu.append(spec_amu)
Nx_mu.append(Nxx_mu)
Nx_amu.append(Nxx_amu)
trax.append(traxx)
spectrum_mu.append(spect_mu)
spectrum_amu.append(spect_amu)
N_mu.append(Nx_mu)
N_amu.append(Nx_amu)
track.append(trax)
nevents = 50000
nfiles = 7000
# Get a histogram that gives the likelihood to encounter the sun at a certain elevation in the sky
hsun = ING.H1D.FromFile("/home/ctoennis/analyses/standard_analysis_framework/Signal_Ingredients/HSun.txt")
hsun.Scale(hsun.nbin*1.0/(hsun.integral))
for filename in os.listdir("/data/ana/analyses/northern_tracks/version-002-p00/"): # loop over MC files
if "MC.npy" in filename and "IC86" in filename: #Only take MC files
infile = np.load("/data/ana/analyses/northern_tracks/version-002-p00/"+filename,"r")
for entry in infile:
erec = entry["logE"]
nutype = entry["trueType"]
weight = entry["orig_OW"]
energy = entry["trueE"]
azt = entry["trueAzi"]
dect = entry["trueZen"]
azr = entry["azi"]
decr = entry["zen"]
k = massnumber
diff = (np.dot(np.array([np.sin(azt)*np.sin(dect),np.cos(azt)*np.sin(dect),np.cos(dect)]),np.array([np.sin(azr)*np.sin(decr),np.cos(azr)*np.sin(decr),np.cos(decr)]))) # Angular separation from the Sun
if diff > 0.998 and energy < masslist[k]: # only look at events in a small ROI around the Sun
for j in range(len(medmasslist)):
for i in range(len(lifetimelist)):
if masslist[k] < medmasslist[j]:
continue
for l in range(len(channellist)):
if not spectrum_mu[k][j][i][l] == "None" and N_mu[k][j][i][l]>0:
if nutype == 68 or nutype == 14:
track[k][j][i][l].Fill(np.power(10.0,erec),weight*spectrum_mu[k][j][i][l].Eval(energy/masslist[k]))
else:
track[k][j][i][l].Fill(np.power(10.0,erec),weight*spectrum_amu[k][j][i][l].Eval(energy/masslist[k]))
k = massnumber
#for k in range(len(masslist)):
for j in range(len(medmasslist)):
for i in range(len(lifetimelist)):
for l in range(len(channellist)):
if not spectrum_mu[k][j][i][l] == "None" and track[k][j][i][l].integral > 0.0:
track[k][j][i][l].Scale(len(track[k][j][i][l].content)*1.0/track[k][j][i][l].integral)
track[k][j][i][l].Write("/home/ctoennis/analyses/standard_analysis_framework/northern_ing/ENERGY_track_m" + str(masslist[k]) + "_med" + str(medmasslist[j]) + "_l" + str(lifetimelist[i]) + "_ch" + str(channellist[l]) + ".txt")