-
Notifications
You must be signed in to change notification settings - Fork 19
/
snowglobes_interface.py
249 lines (217 loc) · 10.9 KB
/
snowglobes_interface.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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
"""
The module ``snewpy.snowglobes_interface`` contains a low-level Python interface for SNOwGLoBES v1.2.
.. note::
Users should only use the high-level interface described above.
This low-level interface is not guaranteed to be stable and may change at
any time without warning, e.g. to support new SNOwGLoBES versions.
"""
from pathlib import Path
import pandas as pd
import numpy as np
import os
import logging
logger = logging.getLogger(__name__)
from tqdm.auto import tqdm
def guess_material(detector):
if detector.startswith('wc') or detector.startswith('ice'):
mat = 'water'
elif detector.startswith('d2O'):
mat = 'heavywater'
elif detector.startswith('ar'):
mat = 'argon'
elif detector.startswith('nova'):
mat = 'nova_soup'
elif detector.startswith('halo'):
mat = 'lead'
elif detector.startswith('scint'):
mat = 'scint'
else:
raise ValueError(f'Please provide material for {detector}')
if '_he' in detector:
mat = mat+'_he'
return mat
class SimpleRate():
def __init__(self, detectors:str="all", detector_effects=True, base_dir:Path=''):
"""Simple rate calculation interface.
Computes expected rate for a detector without using GLOBES. The formula for the rate is
Rate = Sum_i ([cross-section in 10^-38 cm^2] x 10^-38 x [fluence in cm^-2])_i x [smearing matrix]_{ij}
x [target mass in kton] x [Dalton per kton] x [energy bin size in GeV] x [efficiency]
with [target mass in kton] x [Dalton per kton] = number of reference targets in experiment.
On construction the code will read:
* detectors from `<base_dir>/detector_configurations.dat`,
* channels from `<base_dir>/channels/channel_*.dat`
if the detector_effects option is on, also read efficiencies and smearing:
* efficiencies from `<base_dir>/effic/effic_*.dat`,
* smearing matrices from `<base_dir>/smear/smear_*.dat`
After that use :meth:`SimpleRate.run` method to run the simulation for specific detector and flux file.
Parameters
----------
detectors: str
Name of detector. If ``"all"``, will use all detectors supported by SNOwGLoBES.
detector_effects: bool
If true, account for efficiency and smearing. If false, consider a perfect detector.
base_dir: Path or None
Path to the directory where the cross-section, detector, and channel files are located
If empty, try to get it from ``$SNOWGLOBES`` environment var
"""
if not base_dir:
base_dir = os.environ['SNOWGLOBES']
self.base_dir = Path(base_dir)
self._load_detectors(self.base_dir/'detector_configurations.dat', detectors)
self._load_channels(self.base_dir/'channels')
self.efficiencies = None
self.smearings = None
if detector_effects:
self._load_efficiency_vectors(self.base_dir/'effic')
self._load_smearing_matrices(self.base_dir/'smear')
def _load_detectors(self, path:Path, detectors:str):
df = pd.read_table(path,names=['name','mass','factor'], delim_whitespace=True, comment='#')
df['tgt_mass']=df.mass*df.factor
if detectors == 'all':
detectors = list(df['name'])
detectors.remove('d2O')
elif isinstance(detectors, str):
detectors = [detectors]
self.detectors = df.set_index('name').T.filter(items=detectors)
logger.info(f'read masses for detectors {list(self.detectors)}')
logger.debug(f'detectors: {self.detectors}')
def _load_channels(self, chan_dir):
def _read_binning(fname):
with open(fname) as f:
l = f.readline().strip()
if not l.startswith('%'):
l = '% 200 0.0005 0.100 200 0.0005 0.100'
tokens = l.split(' ')[1:]
return dict(zip(['nsamples','smin','smax','nbins','emin','emax'],tokens))
self.channels = {}
self.binning = {}
for f in chan_dir.glob('channels_*.dat'):
material = f.stem[len('channels_'):]
df = pd.read_table(f,delim_whitespace=True, index_col=1, comment='%', names=['name','n','parity','flavor','weight'])
self.channels[material] = df
self.binning[material] = _read_binning(f)
self.materials = list(self.channels.keys())
self.chan_dir = chan_dir
logger.info(f'read channels for materials: {self.materials}')
logger.debug(f'channels: {self.channels}')
def _load_efficiency_vectors(self, path):
result = {}
for detector in self.detectors:
res_det = {}
for file in path.glob(f'effic_*_{detector}.dat'):
channel = file.stem[len('effic_'):-len(detector)-1]
logger.debug(f'Reading file ({detector},{channel}): {file}')
with open(file) as f:
effs = np.fromiter(f.readlines()[0].split("{")[-1].split("}")[0].split(","), float)
res_det[channel]= effs
result[detector]=res_det
self.efficiencies = result
logger.info(f'read efficiencies for detectors: {list(self.efficiencies.keys())}')
logger.debug(f'efficiencies: {self.efficiencies}')
def _load_smearing_matrices(self, path):
result = {}
for detector in self.detectors:
res_det = {}
for file in path.glob(f'smear*_{detector}.dat'):
channel = file.stem[len('smear_'):-len(detector)-1]
with open(file) as f:
lines = f.readlines()[1:-1]
while not "{" in lines[0]: lines = lines[1:]
while not "{" in lines[-1]: lines = lines[:-1]
matrix = np.zeros((len(lines),len(lines)))
for i,l in enumerate(lines):
elements = np.fromiter(l.split("{")[-1].split("}")[0].split(","), float)
matrix[i, int(elements[0]+0.1):int(elements[1]+0.1)+1] = elements[2:]
res_det[channel]= matrix
result[detector]=res_det
self.smearings = result
logger.info(f'read smearing matrices for detectors: {list(self.smearings.keys())}')
logger.debug(f'smearing matrices: {self.smearings}')
def _compute_rates(self, detector, material, flux_file:Path):
flux_file = flux_file.resolve()
fluxes = np.loadtxt(flux_file)
TargetMass = self.detectors[detector].tgt_mass
data = {}
energies = np.linspace(7.49e-4, 9.975e-2, 200) # Use the same energy grid as SNOwGLoBES
if '_he' in detector:
energies = np.linspace(7.49e-4, 19.975e-2, 400) #SNOwGLoBES grid for he configurations
for channel in self.channels[material].itertuples():
xsec_path = f"xscns/xs_{channel.name}.dat"
xsec = np.loadtxt(self.base_dir/xsec_path)
flavor_index = 0 if 'e' in channel.flavor else (1 if 'm' in channel.flavor else 2)
flavor = flavor_index + (3 if channel.parity == '-' else 0)
flux = fluxes[:, (0,1+flavor)]
binsize = energies[1] - energies[0]
# Cross-section in 10^-38 cm^2
xsecs = np.interp(np.log(energies)/np.log(10), xsec[:, 0], xsec[:, 1+flavor], left=0, right=0) * energies
# Fluence (flux integrated over time bin) in cm^-2
# (must be divided by 0.2 MeV to compensate the multiplication in generate_time_series)
fluxs = np.interp(energies, flux[:, 0], flux[:, 1], left=0, right=0)/2e-4
# Rate computation
rates = xsecs * 1e-38 * fluxs * float(TargetMass) * 1./1.661e-33 * binsize
# Weighting
weighted_rates = rates * channel.weight
# Write to dictionary
data[(channel.name,'unsmeared','unweighted')] = rates
data[(channel.name,'unsmeared','weighted')] = weighted_rates
# Add detector effects
if self.smearings and self.efficiencies:
smear,effic = None,None
if channel.name in self.smearings[detector].keys():
smear = self.smearings[detector][channel.name]
else:
smear = np.eye(len(rates))
if channel.name in self.efficiencies[detector].keys():
effic = self.efficiencies[detector][channel.name]
else:
effic = np.ones(len(rates))
rates = np.dot(smear,rates) * effic
weighted_rates = rates * channel.weight
# Write to dictionary
data[(channel.name,'smeared','unweighted')] = rates
data[(channel.name,'smeared','weighted')] = weighted_rates
#collect everything to pandas DataFrame
df = pd.DataFrame(data, index = energies)
df.index.rename('E', inplace=True)
df.columns.rename(['channel','is_smeared','is_weighted'], inplace=True)
return df.reorder_levels([2,1,0], axis='columns')
def run(self, flux_files, detector:str, material:str=None):
""" Compute expected rates for given configuration,
collect the resulting data and return it in `pandas.DataFrame`
Parameters
-----------
flux_files: list(str) or str
An iterable of flux table filenames to process, or a single filename
detector: str
Detector name, SNOwGLoBES style
material: str or None
Material name, SNOwGLoBES style. If None, we'll try to guess it
Returns
--------
list(pd.DataFrame or Exception)
List with the data table for each flux_file, keeping the order.
Each table containing Energy (GeV) as index values,
and number of events for each energy bin, for all interaction channels.
Columns are hierarchical: (is_weighted, channel),
so one can easily access the desired final table.
If run failed with exception, this exception will be returned (not raised).
Raises
------
ValueError
if material or detector value is invalid
"""
if not detector in self.detectors:
raise ValueError(f'Detector "{detector}" is not in {list(self.detectors)}')
if material is None:
material = guess_material(detector)
if not material in self.materials:
raise ValueError(f'Material "{material}" is not in {self.materials}')
if isinstance(flux_files,str):
flux_files = [flux_files]
with tqdm(total=len(flux_files), leave=False, desc='Flux files') as progressbar:
results = []
for flux_file in flux_files:
result = self._compute_rates(detector,material,Path(flux_file))
progressbar.update()
results.append(result)
return results