Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Python data #48

Merged
merged 15 commits into from
Jan 24, 2024
Merged

Python data #48

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
d28de58
Add version manipulation and model finding methods to _util.py
austinschneider Jan 24, 2024
127fc43
Move resources
austinschneider Jan 24, 2024
224103f
Fix path
austinschneider Jan 24, 2024
8ae6db9
Use methods from _util to get latest detector model files
austinschneider Jan 24, 2024
3a92c6c
Fix path
austinschneider Jan 24, 2024
d0ec2be
Fix Example2
austinschneider Jan 24, 2024
1b9598a
Make it black
austinschneider Jan 24, 2024
d29746f
Remove bad DIS splines
austinschneider Jan 24, 2024
8aa190e
Add DIS splines from DUNEAtmo
austinschneider Jan 24, 2024
fe71e62
Use versioned cross sections in DIS example
austinschneider Jan 24, 2024
f3128b7
Add more options when searching for models
austinschneider Jan 24, 2024
ef3af3c
Use versioned cross section paths in Example2
austinschneider Jan 24, 2024
ab24bea
Fix typo
austinschneider Jan 24, 2024
b2e193b
Add comments and fix dev version comparison
austinschneider Jan 24, 2024
6bbd8e2
Fix bugs in LIDarkNews impacting saving of CrossSection tables
nickkamp1 Jan 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 10 additions & 10 deletions projects/detector/private/DetectorModel.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -267,11 +267,11 @@ void DetectorModel::LoadDetectorModel(std::string const & detector_model) {
else if(fexists(path_ + "/densities/" + detector_model + ".dat")) {
fname = path_ + "/densities/" + detector_model + ".dat";
}
else if(fexists(path_ + "/DetectorParams/" + detector_model)) {
fname = path_ + "/DetectorParams/" + detector_model;
else if(fexists(path_ + "/Detectors/" + detector_model)) {
fname = path_ + "/Detectors/" + detector_model;
}
else if(fexists(path_ + "/DetectorParams/" + detector_model + ".dat")) {
fname = path_ + "/DetectorParams/" + detector_model + ".dat";
else if(fexists(path_ + "/Detectors/" + detector_model + ".dat")) {
fname = path_ + "/Detectors/" + detector_model + ".dat";
}
else if(fexists(path_ + "/" + detector_model)) {
fname = path_ + "/" + detector_model;
Expand All @@ -287,7 +287,7 @@ void DetectorModel::LoadDetectorModel(std::string const & detector_model) {

// if the detectormodel file doesn't exist, stop simulation
if(in.fail()){
throw(std::runtime_error("Failed to open " + fname + " Set correct DetectorParamsPath."));
throw(std::runtime_error("Failed to open " + fname + " Set correct DetectorsPath."));
}

ClearSectors();
Expand Down Expand Up @@ -1482,11 +1482,11 @@ void DetectorModel::LoadConcentricShellsFromLegacyFile(std::string model_fname,
else if(fexists(path_ + "/densities/" + model_fname + ".dat")) {
fname = path_ + "/densities/" + model_fname + ".dat";
}
else if(fexists(path_ + "/DetectorParams/" + model_fname)) {
fname = path_ + "/DetectorParams/" + model_fname;
else if(fexists(path_ + "/Detectors/" + model_fname)) {
fname = path_ + "/Detectors/" + model_fname;
}
else if(fexists(path_ + "/DetectorParams/" + model_fname + ".dat")) {
fname = path_ + "/DetectorParams/" + model_fname + ".dat";
else if(fexists(path_ + "/Detectors/" + model_fname + ".dat")) {
fname = path_ + "/Detectors/" + model_fname + ".dat";
}
else if(fexists(path_ + "/" + model_fname)) {
fname = path_ + "/" + model_fname;
Expand All @@ -1502,7 +1502,7 @@ void DetectorModel::LoadConcentricShellsFromLegacyFile(std::string model_fname,

// if the detectormodel file doesn't exist, stop simulation
if(in.fail()){
throw(std::runtime_error("Failed to open " + fname + " Set correct DetectorParamsPath."));
throw(std::runtime_error("Failed to open " + fname + " Set correct DetectorsPath."));
}

ClearSectors();
Expand Down
4 changes: 2 additions & 2 deletions projects/injection/private/test/CCM_HNL_TEST.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,8 +60,8 @@ bool inelastic = true;
std::string tot_xsec_table_path = "/home/nwkamp/Research/Pheno/Neutrissimos2/Sandbox/xsec_tables/tot_xsec_Enu/";
std::string diff_xsec_table_path = "/home/nwkamp/Research/Pheno/Neutrissimos2/Sandbox/xsec_tables/";

std::string material_file = "/home/nwkamp/Research/CCM/DipoleAnalysis/sources/LeptonInjectorDevPrivate/resources/DetectorParams/materials/CCM.dat";
std::string detector_file = "/home/nwkamp/Research/CCM/DipoleAnalysis/sources/LeptonInjectorDevPrivate/resources/DetectorParams/densities/PREM_ccm.dat";
std::string material_file = "/home/nwkamp/Research/CCM/DipoleAnalysis/sources/LeptonInjectorDevPrivate/resources/Detectors/materials/CCM.dat";
std::string detector_file = "/home/nwkamp/Research/CCM/DipoleAnalysis/sources/LeptonInjectorDevPrivate/resources/Detectors/densities/PREM_ccm.dat";

double hnl_mass = 0.01375; // in GeV; The HNL mass we are injecting
double dipole_coupling = 1.0e-6; // in GeV^-1; the effective dipole coupling strength
Expand Down
12 changes: 6 additions & 6 deletions projects/injection/private/test/Injector_TEST.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,18 +212,18 @@ TEST(Injector, Generation)
using ParticleType = Particle::ParticleType;

#ifdef AUSTIN
std::string material_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/DetectorParams/materials/Minerva.dat";
std::string detector_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/DetectorParams/densities/PREM_minerva.dat";
std::string material_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/Detectors/materials/Minerva.dat";
std::string detector_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/Detectors/densities/PREM_minerva.dat";
std::string flux_file = "/home/austin/nu-dipole/fluxes/LE_FHC_numu.txt";
z_samp = false;
in_invGeV = false;
#else
std::string material_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/DetectorParams/materials/Minerva.dat";
std::string detector_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/DetectorParams/densities/PREM_minerva.dat";
std::string material_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/Detectors/materials/Minerva.dat";
std::string detector_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/Detectors/densities/PREM_minerva.dat";
std::string flux_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/Sandbox/NUMI_Flux_Tables/ME_FHC_numu.txt";
if(miniboone) {
material_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/DetectorParams/materials/MiniBooNE.dat";
detector_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/DetectorParams/densities/PREM_miniboone.dat";
material_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/Detectors/materials/MiniBooNE.dat";
detector_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/Detectors/densities/PREM_miniboone.dat";
flux_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/Sandbox/BNB_Flux_Tables/BNB_numu_flux.txt";
inelastic = true;
}
Expand Down
8 changes: 4 additions & 4 deletions projects/interactions/private/test/ElasticScattering_TEST.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,12 +58,12 @@ TEST(ElasticScattering, Generation)
using ParticleType = Particle::ParticleType;

#ifdef AUSTIN
std::string material_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/DetectorParams/materials/Minerva.dat";
std::string detector_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/DetectorParams/densities/PREM_minerva.dat";
std::string material_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/Detectors/materials/Minerva.dat";
std::string detector_file = "/home/austin/programs/LIDUNE/sources/LeptonInjectorDUNE/resources/Detectors/densities/PREM_minerva.dat";
std::string flux_file = "/home/austin/nu-dipole/fluxes/LE_FHC_numu.txt";
#else
std::string material_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/DetectorParams/materials/Minerva.dat";
std::string detector_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/DetectorParams/densities/PREM_minerva.dat";
std::string material_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/Detectors/materials/Minerva.dat";
std::string detector_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/sources/LeptonInjectorDUNE/resources/Detectors/densities/PREM_minerva.dat";
std::string flux_file = "/home/nwkamp/Research/Pheno/Neutrissimos2/Sandbox/NUMI_Flux_Tables/ME_FHC_numu.txt";
#endif

Expand Down
17 changes: 2 additions & 15 deletions python/LIController.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,21 +40,8 @@ def __init__(self, events_to_inject, experiment, seed=0):
self.events = []

# Find the density and materials files
materials_file = os.path.join(
self.resources_dir, "DetectorParams", "materials", f"{experiment}.dat"
)
if experiment in ["ATLAS", "dune", "IceCube"]:
detector_model_file = os.path.join(
self.resources_dir,
"DetectorParams",
"densities",
f"PREM_{experiment}.dat",
)
else:
detector_model_file = os.path.join(
self.resources_dir, "DetectorParams", "densities", f"{experiment}.dat"
)

materials_file = _util.get_material_model_path(experiment)
detector_model_file = _util.get_detector_model_path(experiment)
print(detector_model_file)

self.detector_model = _detector.DetectorModel()
Expand Down
68 changes: 39 additions & 29 deletions python/LIDarkNews.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ def __init__(
if self.table_dir is None:
self.table_dir = os.path.join(
resources_dir,
"CrossSectionTables",
"CrossSections",
"DarkNewsTables",
datetime.datetime.now().strftime("%Y_%m_%d__%H:%M"),
)
Expand Down Expand Up @@ -110,7 +110,9 @@ def __init__(
table_subdirs += "%s_" % str(x)
table_subdirs += "/"
self.decays.append(
PyDarkNewsDecay(dec_case, table_dir=os.path.join(self.table_dir, table_subdirs))
PyDarkNewsDecay(
dec_case, table_dir=os.path.join(self.table_dir, table_subdirs)
)
)

def SaveCrossSectionTables(self, fill_tables_at_exit=True):
Expand Down Expand Up @@ -176,7 +178,9 @@ def __init__(
total_xsec_file = os.path.join(self.table_dir, "total_cross_sections.npy")
if os.path.exists(total_xsec_file):
self.total_cross_section_table = np.load(total_xsec_file)
diff_xsec_file = os.path.join(self.table_dir, "differential_cross_sections.npy")
diff_xsec_file = os.path.join(
self.table_dir, "differential_cross_sections.npy"
)
if os.path.exists(diff_xsec_file):
self.differential_cross_section_table = np.load(diff_xsec_file)

Expand All @@ -185,29 +189,31 @@ def __init__(
# Sorts and redefines scipy interpolation objects
def _redefine_interpolation_objects(self, total=False, diff=False):
if total:
self.total_cross_section_table = np.sort(
self.total_cross_section_table, axis=0
if len(self.total_cross_section_table) <= 1: return
idxs = np.argsort(
self.total_cross_section_table[:,0]
)
self.total_cross_section_table = self.total_cross_section_table[idxs]
self.total_cross_section_interpolator = CubicSpline(
self.total_cross_section_table[:, 0],
self.total_cross_section_table[:, 1],
)
if len(self.total_cross_section_table) > 1:
self.total_cross_section_interpolator = CubicSpline(
self.total_cross_section_table[:, 0],
self.total_cross_section_table[:, 1],
)
if diff:
self.differential_cross_section_table = np.sort(
self.differential_cross_section_table, axis=0
if len(self.differential_cross_section_table) <= 1: return
idxs = np.lexsort(
(self.differential_cross_section_table[:,1],
self.differential_cross_section_table[:,0])
)
if len(self.differential_cross_section_table) > 1:
# If we only have two energy points, don't try to construct interpolator
if len(np.unique(self.differential_cross_section_table[:, 0])) <= 2:
return
self.differential_cross_section_interpolator = (
CloughTocher2DInterpolator(
self.differential_cross_section_table[:, :2],
self.differential_cross_section_table[:, 2],
rescale=True,
)
self.differential_cross_section_table = self.differential_cross_section_table[idxs]
# If we only have two energy points, don't try to construct interpolator
if len(np.unique(self.differential_cross_section_table[:, 0])) <= 2: return
self.differential_cross_section_interpolator = (
CloughTocher2DInterpolator(
self.differential_cross_section_table[:, :2],
self.differential_cross_section_table[:, 2],
rescale=True,
)
)

# Check whether we have close-enough entries in the intrepolation tables
def _interpolation_flags(self, inputs, mode):
Expand All @@ -231,10 +237,10 @@ def _interpolation_flags(self, inputs, mode):
return False, False, -1

# bools to keep track of whether to use a single point or interpolate
UseSinglePoint = True
UseSinglePoint = False
Interpolate = True
# First check whether we have a close-enough single point
closest_idx = np.argmin(np.sum(np.abs(interp_table[:, :-1] - inputs)))
closest_idx = np.argmin(np.sum(np.abs(interp_table[:, :-1] - inputs),axis=-1))
diff = (interp_table[closest_idx, :-1] - inputs) / inputs
if np.all(np.abs(diff) < self.tolerance):
UseSinglePoint = True
Expand Down Expand Up @@ -316,7 +322,7 @@ def FillInterpolationTables(self, total=True, diff=True):
self.total_cross_section_table, [[E, xsec]], axis=0
)
num_added_points += 1
E *= 1 + self.interp_tolerance
E *= (1 + self.interp_tolerance)
if diff:
# interaction record to calculate Q2 bounds
interaction = LI.dataclasses.InteractionRecord()
Expand Down Expand Up @@ -348,20 +354,24 @@ def FillInterpolationTables(self, total=True, diff=True):
axis=0,
)
num_added_points += 1
z *= 1 + self.interp_tolerance
E *= 1 + self.interp_tolerance
z *= (1 + self.interp_tolerance)
E *= (1 + self.interp_tolerance)
self._redefine_interpolation_objects(total=total, diff=diff)
return num_added_points

# Saves the tables for the scipy interpolation objects
def SaveInterpolationTables(self, total=True, diff=True):
if total:
self._redefine_interpolation_objects(total=True)
with open(os.path.join(self.table_dir, "total_cross_sections.npy"), "wb") as f:
with open(
os.path.join(self.table_dir, "total_cross_sections.npy"), "wb"
) as f:
np.save(f, self.total_cross_section_table)
if diff:
self._redefine_interpolation_objects(diff=True)
with open(os.path.join(self.table_dir, "differential_cross_sections.npy"), "wb") as f:
with open(
os.path.join(self.table_dir, "differential_cross_sections.npy"), "wb"
) as f:
np.save(f, self.differential_cross_section_table)

##### START METHODS FOR SERIALIZATION #########
Expand Down
Loading
Loading