Skip to content

Commit

Permalink
Merge pull request #48 from Harvard-Neutrino/dev/python_data
Browse files Browse the repository at this point in the history 
Python data
  • Loading branch information
@austinschneider
austinschneider authored Jan 24, 2024
2 parents d89559c + 6bbd8e2 commit de8c22e
Show file tree
Hide file tree
Showing 52 changed files with 735 additions and 913 deletions.
20 changes: 10 additions & 10 deletions projects/detector/private/DetectorModel.cxx
View file
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
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
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
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
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
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

0 comments on commit de8c22e

Please sign in to comment.