Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history
for more information, see https://pre-commit.ci
  • Loading branch information
pre-commit-ci[bot] committed Oct 2, 2024
1 parent 975ef57 commit c3163f5
Showing 1 changed file with 26 additions and 21 deletions.
47 changes: 26 additions & 21 deletions src/simweights/_genie_weighter.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,20 +8,21 @@
import numpy as np

from ._generation_surface import GenerationSurface, generation_surface
from ._nugen_weighter import nugen_spatial, nugen_spectrum
from ._powerlaw import PowerLaw
from ._spatial import CircleInjector
from ._utils import Column, Const, constcol, get_column, get_table, has_column, has_table
from ._nugen_weighter import nugen_spatial, nugen_spectrum
from ._weighter import Weighter

def genie_icetray_surface(mcweightdict: Iterable[Mapping[str, float]],
geniedict: Iterable[Mapping[str, float]],
nufraction: float = 0.7) -> GenerationSurface:

def genie_icetray_surface(
mcweightdict: Iterable[Mapping[str, float]], geniedict: Iterable[Mapping[str, float]], nufraction: float = 0.7
) -> GenerationSurface:
"""Inspect the rows of a GENIE-icetray's I3MCWeightDict table object to generate a surface object.
This is a bit of a pain: the oscillations group historically produced 4-5 energy bands with varying
generation parameters, then merged them all into one "file". Because of this, we need to check the
neutrino type, volume, and spectral index to get the correct surfaces. The type weight also isn't
neutrino type, volume, and spectral index to get the correct surfaces. The type weight also isn't
stored in the files: this was fixed to 70/30 for oscillation-produced genie-icetray files.
"""
pdgid = get_column(geniedict, "neu")
Expand All @@ -40,14 +41,15 @@ def genie_icetray_surface(mcweightdict: Iterable[Mapping[str, float]],
for pid, _, _ in unique_schemes:
mask = pid == pdgid

type_weight = nufraction if pid>0 else 1-nufraction
type_weight = nufraction if pid > 0 else 1 - nufraction
primary_type = int(constcol(geniedict, "neu", mask))
n_events = type_weight * constcol(mcweightdict, "NEvents", mask)
surfaces.append(n_events * generation_surface(primary_type, Column("wght"), spectrum, spatial))
ret = sum(surfaces)
assert isinstance(ret, GenerationSurface)
return ret


def genie_reader_surface(table: Iterable[Mapping[str, float]]) -> GenerationSurface:
"""Inspect the rows of a GENIE S-Frame table object to generate a surface object."""
surfaces = []
Expand Down Expand Up @@ -78,8 +80,7 @@ def genie_reader_surface(table: Iterable[Mapping[str, float]]) -> GenerationSurf
return retval


def GenieWeighter(file_obj: Any,
nfiles:float | None = None) -> Weighter: # noqa: N802
def GenieWeighter(file_obj: Any, nfiles: float | None = None) -> Weighter: # noqa: N802
# pylint: disable=invalid-name
"""Weighter for GENIE simulation.
Expand All @@ -89,10 +90,12 @@ def GenieWeighter(file_obj: Any,
if has_table(file_obj, "I3GenieInfo"):
# Branch for newer genie-reader files
if nfiles is not None:
msg = (f"GenieWeighter received an nfiles={nfiles}, but `{getattr(file_obj, 'filename', '<NONE>')}` "
"was produced with genie-reader instead of genie-icetray. We expect to read the number of "
"files from the number of observed S-frames in the file, so this is unnecessary. Do not pass "
"in a value for nfiles for genie-reader files.")
msg = (
f"GenieWeighter received an nfiles={nfiles}, but `{getattr(file_obj, 'filename', '<NONE>')}` "
"was produced with genie-reader instead of genie-icetray. We expect to read the number of "
"files from the number of observed S-frames in the file, so this is unnecessary. Do not pass "
"in a value for nfiles for genie-reader files."
)
raise RuntimeError(msg)

info_table = get_table(file_obj, "I3GenieInfo")
Expand All @@ -115,25 +118,27 @@ def GenieWeighter(file_obj: Any,
else:
# Branch for older genie-icetray files
if nfiles is None:
msg = (f"GenieWeighter received an nfiles={nfiles}, but `{getattr(file_obj, 'filename', '<NONE>')}` "
"was produced with genie-reader instead of genie-icetray. We expect to read the number of "
"files from the number of observed S-frames in the file, so this is unnecessary. Do not pass "
"in a value for nfiles for genie-reader files.")
msg = (
f"GenieWeighter received an nfiles={nfiles}, but `{getattr(file_obj, 'filename', '<NONE>')}` "
"was produced with genie-reader instead of genie-icetray. We expect to read the number of "
"files from the number of observed S-frames in the file, so this is unnecessary. Do not pass "
"in a value for nfiles for genie-reader files."
)
raise RuntimeError(msg)

weight_table = get_table(file_obj, "I3MCWeightDict")
result_table = get_table(file_obj, "I3GENIEResultDict")
surface = nfiles * genie_icetray_surface(weight_table, result_table)

momentum_vec = np.array([get_column(result_table, 'pxv'),
get_column(result_table, 'pyv'),
get_column(result_table, 'pzv')])
cos_zen = -1 * get_column(result_table, 'pzv') / np.sum(momentum_vec**2, axis=0)**0.5
momentum_vec = np.array(
[get_column(result_table, "pxv"), get_column(result_table, "pyv"), get_column(result_table, "pzv")]
)
cos_zen = -1 * get_column(result_table, "pzv") / np.sum(momentum_vec**2, axis=0) ** 0.5

weighter = Weighter([file_obj], surface)
weighter.add_weight_column("pdgid", get_column(result_table, "neu").astype(np.int32))
weighter.add_weight_column("energy", get_column(result_table, "Ev"))
weighter.add_weight_column("cos_zen", cos_zen)
weighter.add_weight_column("wght", get_column(result_table, "wght")*get_column(result_table, "_glbprbscale"))
weighter.add_weight_column("wght", get_column(result_table, "wght") * get_column(result_table, "_glbprbscale"))

return weighter

0 comments on commit c3163f5

Please sign in to comment.