From a7127a16b36fe0efd3eef910876fc504a86b7cca Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sat, 20 Jan 2024 16:57:41 -0700 Subject: [PATCH 01/14] Cmake should install the entire python directory --- CMakeLists.txt | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0dfafd561..5e742e8d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -226,11 +226,7 @@ if(DEFINED SKBUILD) LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) install(TARGETS injection LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) - install(FILES ${PROJECT_SOURCE_DIR}/python/__init__.py - DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) - install(FILES ${PROJECT_SOURCE_DIR}/python/LIDarkNews.py - DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) - install(FILES ${PROJECT_SOURCE_DIR}/python/LIController.py + install(DIRECTORY ${PROJECT_SOURCE_DIR}/python/ DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) endif() From f58fcfe31facc399433e6ccdd9d941261ec5748f Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 14:17:31 -0700 Subject: [PATCH 02/14] Install resources to python directory --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5e742e8d5..918862368 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -228,6 +228,8 @@ if(DEFINED SKBUILD) LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) install(DIRECTORY ${PROJECT_SOURCE_DIR}/python/ DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) + install(DIRECTORY ${PROJECT_SOURCE_DIR}/resources + DESTINATION ${CMAKE_INSTALL_LIBDIR}/leptoninjector) endif() # Export targets for use in downstream CMake projects From 6995d55fc2357e94f486a7957bddf2d0e10443b1 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:03:58 -0700 Subject: [PATCH 03/14] Use relative imports. Find the resources folder more robustly. Make it black. --- python/LIController.py | 342 ++++++++++++++++++++++++----------------- 1 file changed, 204 insertions(+), 138 deletions(-) diff --git a/python/LIController.py b/python/LIController.py index 4d23f7ff9..c79fb2383 100644 --- a/python/LIController.py +++ b/python/LIController.py @@ -2,21 +2,24 @@ import h5py import numpy as np -import leptoninjector as LI -from LIDarkNews import PyDarkNewsInteractionCollection +from . import utilities as _utilities +from . import detector as _detector +from . import injection as _injection +from . import distributions as _distributions +from . import dataclasses as _dataclasses +from . import interactions as _interactions + +from . import _util + +from .LIDarkNews import PyDarkNewsInteractionCollection # For determining fiducial volume of different experiments -fid_vol_dict = {'MiniBooNE':'fid_vol', - 'CCM':'ccm_inner_argon', - 'MINERvA':'fid_vol'} +fid_vol_dict = {"MiniBooNE": "fid_vol", "CCM": "ccm_inner_argon", "MINERvA": "fid_vol"} + # Parent python class for handling event generation class LIController: - - def __init__(self, - events_to_inject, - experiment, - seed=0): + def __init__(self, events_to_inject, experiment, seed=0): """ LI class constructor. :param int event_to_inject: number of events to generate @@ -24,9 +27,11 @@ def __init__(self, :param int seed: Optional random number generator seed """ + self.resources_dir = _util.resource_package_dir() + # Initialize a random number generator - self.random = LI.utilities.LI_random(seed) - + self.random = _utilities.LI_random(seed) + # Save number of events to inject self.events_to_inject = events_to_inject self.experiment = experiment @@ -35,32 +40,44 @@ def __init__(self, self.events = [] # Find the density and materials files - self.LI_SRC = os.environ.get('LEPTONINJECTOR_SRC') - materials_file = self.LI_SRC + '/resources/DetectorParams/materials/%s.dat'%experiment - if experiment in ['ATLAS','dune']: - detector_model_file = self.LI_SRC + '/resources/DetectorParams/densities/PREM_%s.dat'%experiment + 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 = self.LI_SRC + '/resources/DetectorParams/densities/%s.dat'%experiment - - self.detector_model = LI.detector.DetectorModel() + detector_model_file = os.path.join( + self.resources_dir, "DetectorParams", "densities", f"{experiment}.dat" + ) + + print(detector_model_file) + + self.detector_model = _detector.DetectorModel() self.detector_model.LoadMaterialModel(materials_file) self.detector_model.LoadDetectorModel(detector_model_file) # Define the primary injection and physical process - self.primary_injection_process = LI.injection.InjectionProcess() - self.primary_physical_process = LI.injection.PhysicalProcess() - + self.primary_injection_process = _injection.InjectionProcess() + self.primary_physical_process = _injection.PhysicalProcess() + # Define lists for the secondary injection and physical processes self.secondary_injection_processes = [] self.secondary_physical_processes = [] - def SetProcesses(self, - primary_type, - primary_injection_distributions, - primary_physical_distributions, - secondary_types=[], - secondary_injection_distributions=[], - secondary_physical_distributions=[]): + def SetProcesses( + self, + primary_type, + primary_injection_distributions, + primary_physical_distributions, + secondary_types=[], + secondary_injection_distributions=[], + secondary_physical_distributions=[], + ): """ LI process setter. :param ParticleType primary_type: The primary particle being generated @@ -76,81 +93,97 @@ def SetProcesses(self, self.primary_physical_process.primary_type = primary_type # Add all injection distributions - for _,idist in primary_injection_distributions.items(): + for _, idist in primary_injection_distributions.items(): self.primary_injection_process.AddInjectionDistribution(idist) # Add all physical distributions - for _,pdist in primary_physical_distributions.items(): + for _, pdist in primary_physical_distributions.items(): self.primary_physical_process.AddPhysicalDistribution(pdist) - + # Default injection distributions - if 'target' not in primary_injection_distributions.keys(): - self.primary_injection_process.AddInjectionDistribution(LI.distributions.TargetAtRest()) - if 'helicity' not in primary_injection_distributions.keys(): - self.primary_injection_process.AddInjectionDistribution(LI.distributions.PrimaryNeutrinoHelicityDistribution()) - + if "target" not in primary_injection_distributions.keys(): + self.primary_injection_process.AddInjectionDistribution( + _distributions.TargetAtRest() + ) + if "helicity" not in primary_injection_distributions.keys(): + self.primary_injection_process.AddInjectionDistribution( + _distributions.PrimaryNeutrinoHelicityDistribution() + ) + # Default injection distributions - if 'target' not in primary_physical_distributions.keys(): - self.primary_physical_process.AddPhysicalDistribution(LI.distributions.TargetAtRest()) - if 'helicity' not in primary_physical_distributions.keys(): - self.primary_physical_process.AddPhysicalDistribution(LI.distributions.PrimaryNeutrinoHelicityDistribution()) + if "target" not in primary_physical_distributions.keys(): + self.primary_physical_process.AddPhysicalDistribution( + _distributions.TargetAtRest() + ) + if "helicity" not in primary_physical_distributions.keys(): + self.primary_physical_process.AddPhysicalDistribution( + _distributions.PrimaryNeutrinoHelicityDistribution() + ) # Loop through possible secondary interactions - for i_sec,secondary_type in enumerate(secondary_types): - secondary_injection_process = LI.injection.InjectionProcess() - secondary_physical_process = LI.injection.PhysicalProcess() + for i_sec, secondary_type in enumerate(secondary_types): + secondary_injection_process = _injection.InjectionProcess() + secondary_physical_process = _injection.PhysicalProcess() secondary_injection_process.primary_type = secondary_type secondary_physical_process.primary_type = secondary_type - # Add all injection distributions + # Add all injection distributions for idist in secondary_injection_distributions[i_sec]: secondary_injection_process.AddInjectionDistribution(idist) # Add all physical distributions for pdist in secondary_physical_distributions[i_sec]: secondary_physical_process.AddPhysicalDistribution(pdist) - + # Add the position distribution fid_vol = self.GetFiducialVolume() if fid_vol is not None: - secondary_injection_process.AddInjectionDistribution(LI.distributions.SecondaryPositionDistribution(fid_vol)) + secondary_injection_process.AddInjectionDistribution( + _distributions.SecondaryPositionDistribution(fid_vol) + ) else: - secondary_injection_process.AddInjectionDistribution(LI.distributions.SecondaryPositionDistribution()) - + secondary_injection_process.AddInjectionDistribution( + _distributions.SecondaryPositionDistribution() + ) + self.secondary_injection_processes.append(secondary_injection_process) self.secondary_physical_processes.append(secondary_physical_process) - def InputDarkNewsModel(self, - primary_type, - table_dir, - model_kwargs): + def InputDarkNewsModel(self, primary_type, table_dir, model_kwargs): """ Sets up the relevant processes and cross section/decay objects related to a provided DarkNews model dictionary. Will handle the primary cross section collection as well as the entire list of secondary processes - :param LI.dataclasses.Particle.ParticleType primary_type: primary particle to be generated + :param _dataclasses.Particle.ParticleType primary_type: primary particle to be generated :param string table_dir: Directory for storing cross section and decay tables :param dict model_kwargs: The dict of DarkNews model parameters """ # Add nuclear targets to the model arguments - model_kwargs['nuclear_targets'] = self.GetDetectorModelTargets()[1] + model_kwargs["nuclear_targets"] = self.GetDetectorModelTargets()[1] # Initialize DarkNews cross sections and decays - self.DN_processes = PyDarkNewsInteractionCollection(table_dir=table_dir, - **model_kwargs) - + self.DN_processes = PyDarkNewsInteractionCollection( + table_dir=table_dir, **model_kwargs + ) + # Initialize primary InteractionCollection # Loop over available cross sections and save those which match primary type primary_cross_sections = [] for cross_section in self.DN_processes.cross_sections: - if primary_type == LI.dataclasses.Particle.ParticleType(cross_section.ups_case.nu_projectile.pdgid): + if primary_type == _dataclasses.Particle.ParticleType( + cross_section.ups_case.nu_projectile.pdgid + ): primary_cross_sections.append(cross_section) - primary_interaction_collection = LI.interactions.InteractionCollection(primary_type, primary_cross_sections) - + primary_interaction_collection = _interactions.InteractionCollection( + primary_type, primary_cross_sections + ) + # Initialize secondary processes and define secondary InteractionCollection objects secondary_decays = {} # Also keep track of the minimum decay width for defining the position distribution later self.DN_min_decay_width = 0 # Loop over available decays, group by parent type for decay in self.DN_processes.decays: - secondary_type = LI.dataclasses.Particle.ParticleType(decay.dec_case.nu_parent.pdgid) + secondary_type = _dataclasses.Particle.ParticleType( + decay.dec_case.nu_parent.pdgid + ) if secondary_type not in secondary_decays.keys(): secondary_decays[secondary_type] = [] secondary_decays[secondary_type].append(decay) @@ -159,30 +192,37 @@ def InputDarkNewsModel(self, self.DN_min_decay_width = total_decay_width # Now make the list of secondary cross section collections # Add new secondary injection and physical processes at the same time - fid_vol = self.GetFiducialVolume() # find fiducial volume for secondary position distirbutions + fid_vol = ( + self.GetFiducialVolume() + ) # find fiducial volume for secondary position distirbutions secondary_interaction_collections = [] - for secondary_type,decay_list in secondary_decays.items(): - + for secondary_type, decay_list in secondary_decays.items(): # Define a sedcondary injection distribution - secondary_injection_process = LI.injection.InjectionProcess() - secondary_physical_process = LI.injection.PhysicalProcess() + secondary_injection_process = _injection.InjectionProcess() + secondary_physical_process = _injection.PhysicalProcess() secondary_injection_process.primary_type = secondary_type secondary_physical_process.primary_type = secondary_type - + # Add the secondary position distribution if fid_vol is not None: - secondary_injection_process.AddInjectionDistribution(LI.distributions.SecondaryPositionDistribution(fid_vol)) + secondary_injection_process.AddInjectionDistribution( + _distributions.SecondaryPositionDistribution(fid_vol) + ) else: - secondary_injection_process.AddInjectionDistribution(LI.distributions.SecondaryPositionDistribution()) - + secondary_injection_process.AddInjectionDistribution( + _distributions.SecondaryPositionDistribution() + ) + self.secondary_injection_processes.append(secondary_injection_process) self.secondary_physical_processes.append(secondary_physical_process) - secondary_interaction_collections.append(LI.interactions.InteractionCollection(secondary_type, decay_list)) + secondary_interaction_collections.append( + _interactions.InteractionCollection(secondary_type, decay_list) + ) - self.SetCrossSections(primary_interaction_collection,secondary_interaction_collections) - - + self.SetInteractions( + primary_interaction_collection, secondary_interaction_collections + ) def GetFiducialVolume(self): """ @@ -191,7 +231,7 @@ def GetFiducialVolume(self): fid_vol = None for sector in self.detector_model.Sectors: if self.experiment in fid_vol_dict.keys(): - if sector.name==fid_vol_dict[self.experiment]: + if sector.name == fid_vol_dict[self.experiment]: fid_vol = sector.geo return fid_vol @@ -208,33 +248,40 @@ def GetDetectorModelTargets(self): for _target in self.detector_model.Materials.GetMaterialTargets(count): if _target not in targets: targets.append(_target) - if str(_target).find('Nucleus') == -1: + if str(_target).find("Nucleus") == -1: continue else: - target_str = str(_target)[str(_target).find('Type')+5:str(_target).find('Nucleus')] - if target_str == 'H': target_str = 'H1' + target_str = str(_target)[ + str(_target).find("Type") + 5 : str(_target).find("Nucleus") + ] + if target_str == "H": + target_str = "H1" if target_str not in target_strs: target_strs.append(target_str) count += 1 return targets, target_strs - - def SetCrossSections(self, - primary_interaction_collection, - secondary_interaction_collections): + + def SetInteractions( + self, primary_interaction_collection, secondary_interaction_collections=None + ): """ Set cross sections for the primary and secondary processes :param InteractionCollection primary_interaction_collection: The cross section collection for the primary process :param list secondary_interaction_collections: The list of cross section collections for the primary process """ + if secondary_interaction_collections is None: + secondary_interaction_collections = [] + # Set primary cross sections self.primary_injection_process.interactions = primary_interaction_collection self.primary_physical_process.interactions = primary_interaction_collection - + # Loop through secondary processes - for sec_inj,sec_phys in zip(self.secondary_injection_processes, - self.secondary_physical_processes): - assert(sec_inj.primary_type == sec_phys.primary_type) - record = LI.dataclasses.InteractionRecord() + for sec_inj, sec_phys in zip( + self.secondary_injection_processes, self.secondary_physical_processes + ): + assert sec_inj.primary_type == sec_phys.primary_type + record = _dataclasses.InteractionRecord() record.signature.primary_type = sec_inj.primary_type found_collection = False # Loop through possible seconday cross sections @@ -244,78 +291,97 @@ def SetCrossSections(self, sec_inj.interactions = sec_xs sec_phys.interactions = sec_xs found_collection = True - if(not found_collection): - print('Couldn\'t find cross section collection for secondary particle %s; Exiting'%record.primary_type) + if not found_collection: + print( + "Couldn't find cross section collection for secondary particle %s; Exiting" + % record.primary_type + ) exit(0) - - - def Initialize(self): - # Define stopping condition # TODO: make this more general def StoppingCondition(datum): return True - + # Define the injector object - self.injector = LI.injection.Injector(self.events_to_inject, - self.detector_model, - self.primary_injection_process, - self.secondary_injection_processes, - self.random) - + self.injector = _injection.Injector( + self.events_to_inject, + self.detector_model, + self.primary_injection_process, + self.secondary_injection_processes, + self.random, + ) + self.injector.SetStoppingCondition(StoppingCondition) # Define the weighter object - self.weighter = LI.injection.LeptonTreeWeighter([self.injector], - self.detector_model, - self.primary_physical_process, - self.secondary_physical_processes) - - def GenerateEvents(self,N=None): + self.weighter = _injection.LeptonTreeWeighter( + [self.injector], + self.detector_model, + self.primary_physical_process, + self.secondary_physical_processes, + ) + + def GenerateEvents(self, N=None): if N is None: N = self.events_to_inject count = 0 - while (self.injector.InjectedEvents() < self.events_to_inject) \ - and (count < N): - print('Injecting Event %d/%d '%(count,N),end='\r') + while (self.injector.InjectedEvents() < self.events_to_inject) and (count < N): + print("Injecting Event %d/%d " % (count, N), end="\r") tree = self.injector.GenerateEvent() self.events.append(tree) count += 1 - DN_processes = getattr(self,"DN_processes") - if DN_processes is not None: - DN_processes.SaveCrossSectionTables() + if hasattr(self, "DN_processes"): + self.DN_processes.SaveCrossSectionTables() return self.events - def SaveEvents(self,filename): - fout = h5py.File(filename,'w') - fout.attrs['num_events'] = len(self.events) - for ie,event in enumerate(self.events): - print('Saving Event %d/%d '%(ie,len(self.events)),end='\r') - event_group = fout.require_group("event%d"%ie) - event_group.attrs['event_weight'] = self.weighter.EventWeight(event) - event_group.attrs['num_interactions'] = len(event.tree) - for id,datum in enumerate(event.tree): - interaction_group = event_group.require_group("interaction%d"%id) - + def SaveEvents(self, filename): + fout = h5py.File(filename, "w") + fout.attrs["num_events"] = len(self.events) + for ie, event in enumerate(self.events): + print("Saving Event %d/%d " % (ie, len(self.events)), end="\r") + event_group = fout.require_group("event%d" % ie) + event_group.attrs["event_weight"] = self.weighter.EventWeight(event) + event_group.attrs["num_interactions"] = len(event.tree) + for id, datum in enumerate(event.tree): + interaction_group = event_group.require_group("interaction%d" % id) + # Add metadata on interaction signature - interaction_group.attrs['primary_type'] = str(datum.record.signature.primary_type) - interaction_group.attrs['target_type'] = str(datum.record.signature.target_type) - for isec,secondary in enumerate(datum.record.signature.secondary_types): - interaction_group.attrs['secondary_type%d'%isec] = str(secondary) + interaction_group.attrs["primary_type"] = str( + datum.record.signature.primary_type + ) + interaction_group.attrs["target_type"] = str( + datum.record.signature.target_type + ) + for isec, secondary in enumerate( + datum.record.signature.secondary_types + ): + interaction_group.attrs["secondary_type%d" % isec] = str(secondary) # Save vertex as dataset - interaction_group.create_dataset('vertex',data=np.array(datum.record.interaction_vertex,dtype=float)) - + interaction_group.create_dataset( + "vertex", + data=np.array(datum.record.interaction_vertex, dtype=float), + ) + # Save each four-momenta as a dataset - interaction_group.create_dataset('primary_momentum',data=np.array(datum.record.primary_momentum,dtype=float)) - interaction_group.create_dataset('target_momentum',data=np.array(datum.record.primary_momentum,dtype=float)) - for isec_momenta,sec_momenta in enumerate(datum.record.secondary_momenta): - interaction_group.create_dataset('secondary_momentum%d'%isec_momenta,data=np.array(sec_momenta,dtype=float)) + interaction_group.create_dataset( + "primary_momentum", + data=np.array(datum.record.primary_momentum, dtype=float), + ) + interaction_group.create_dataset( + "target_momentum", + data=np.array(datum.record.primary_momentum, dtype=float), + ) + for isec_momenta, sec_momenta in enumerate( + datum.record.secondary_momenta + ): + interaction_group.create_dataset( + "secondary_momentum%d" % isec_momenta, + data=np.array(sec_momenta, dtype=float), + ) fout.close() - DN_processes = getattr(self,"DN_processes") - if DN_processes is not None: - DN_processes.SaveCrossSectionTables() - + if hasattr(self, "DN_processes"): + self.DN_processes.SaveCrossSectionTables() From d107f89ac738ea67931e53cff7ee3d3da0b1ac1e Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:04:23 -0700 Subject: [PATCH 04/14] Add argument annotations to the DISFromSpline constructors --- .../private/pybindings/DISFromSpline.h | 55 +++++++++++++++++-- 1 file changed, 49 insertions(+), 6 deletions(-) diff --git a/projects/interactions/private/pybindings/DISFromSpline.h b/projects/interactions/private/pybindings/DISFromSpline.h index 566d9edfc..01ca7990d 100644 --- a/projects/interactions/private/pybindings/DISFromSpline.h +++ b/projects/interactions/private/pybindings/DISFromSpline.h @@ -21,13 +21,56 @@ void register_DISFromSpline(pybind11::module_ & m) { class_, CrossSection> disfromspline(m, "DISFromSpline"); disfromspline + .def(init<>()) - .def(init, std::vector, int, double, double, std::set, std::set, std::string>()) - .def(init, std::vector, int, double, double, std::vector, std::vector, std::string>()) - .def(init, std::set, std::string>()) - .def(init, std::set, std::string>()) - .def(init, std::vector, std::string>()) - .def(init, std::vector, std::string>()) + .def(init, std::vector, int, double, double, std::set, std::set, std::string>(), + arg("total_xs_data"), + arg("differential_xs_data"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::vector, int, double, double, std::vector, std::vector, std::string>(), + arg("total_xs_data"), + arg("differential_xs_data"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::set, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::set, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::vector, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("interaction"), + arg("target_mass"), + arg("minimum_Q2"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) + .def(init, std::vector, std::string>(), + arg("total_xs_filename"), + arg("differential_xs_filename"), + arg("primary_types"), + arg("target_types"), + arg("units") = std::string("cm")) .def(self == self) .def("TotalCrossSection",overload_cast(&DISFromSpline::TotalCrossSection, const_)) .def("TotalCrossSection",overload_cast(&DISFromSpline::TotalCrossSection, const_)) From 9950257a57659b52ab858c0f40e596887a7b3ff0 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:05:20 -0700 Subject: [PATCH 05/14] Get IceCube example working --- resources/Examples/Example1/DIS_IceCube.py | 57 ++++++++++++---------- 1 file changed, 30 insertions(+), 27 deletions(-) diff --git a/resources/Examples/Example1/DIS_IceCube.py b/resources/Examples/Example1/DIS_IceCube.py index d53dc909a..71b6b3d9d 100644 --- a/resources/Examples/Example1/DIS_IceCube.py +++ b/resources/Examples/Example1/DIS_IceCube.py @@ -3,64 +3,67 @@ import leptoninjector as LI import sys -sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python') -from LIController import LIController - -LI_SRC = os.environ.get('LEPTONINJECTOR_SRC') +from leptoninjector import _util +from leptoninjector.LIController import LIController + +resources_dir = _util.resource_package_dir() # Number of events to inject events_to_inject = 1000 # Expeirment to run -experiment = 'IceCube' +experiment = "IceCube" # Define the controller -controller = LIController(events_to_inject, - experiment) +controller = LIController(events_to_inject, experiment) # Particle to inject primary_type = LI.dataclasses.Particle.ParticleType.NuMu # Cross Section Model -xsfiledir = LI_SRC+'resources/CrossSectionTables/DISSplines/' +xsfiledir = os.path.join(resources_dir, "CrossSectionTables", "DISSplines") target_type = LI.dataclasses.Particle.ParticleType.Nucleon -DIS_xs = LI.interactions.DISFromSpline(xsfiledir+'test_xs.fits', - xsfiledir+'test_xs_total.fits', - [primary_type], - [target_type]) +DIS_xs = LI.interactions.DISFromSpline( + os.path.join(xsfiledir, "test_xs.fits"), + os.path.join(xsfiledir, "test_xs_total.fits"), + [primary_type], + [target_type], +) -primary_xs = LI.interactions.CrossSectionCollection(primary_type, [DIS_xs]) -controller.SetCrossSections(primary_xs) +primary_xs = LI.interactions.InteractionCollection(primary_type, [DIS_xs]) +controller.SetInteractions(primary_xs) # Primary distributions primary_injection_distributions = {} primary_physical_distributions = {} # energy distribution -edist = LI.distributions.PowerLaw(2,1e3,1e6) -primary_injection_distributions['energy'] = edist -primary_physical_distributions['energy'] = edist +edist = LI.distributions.PowerLaw(2, 1e3, 1e6) +primary_injection_distributions["energy"] = edist +primary_physical_distributions["energy"] = edist # direction distribution direction_distribution = LI.distributions.IsotropicDirection() -primary_injection_distributions['direction'] = direction_distribution -primary_physical_distributions['direction'] = direction_distribution +primary_injection_distributions["direction"] = direction_distribution +primary_physical_distributions["direction"] = direction_distribution # position distribution muon_range_func = LI.distributions.LeptonDepthFunction() -position_distribution = LI.distributions.ColumnDepthPositionDistribution(600, 600., - muon_range_func, - set(controller.GetEarthModelTargets()[0])) -primary_injection_distributions['position'] = position_distribution +position_distribution = LI.distributions.ColumnDepthPositionDistribution( + 600, 600.0, muon_range_func, set(controller.GetDetectorModelTargets()[0]) +) +primary_injection_distributions["position"] = position_distribution # SetProcesses -controller.SetProcesses(primary_type, - primary_injection_distributions, - primary_physical_distributions) +controller.SetProcesses( + primary_type, primary_injection_distributions, primary_physical_distributions +) controller.Initialize() events = controller.GenerateEvents() -controller.SaveEvents('output/MINERvA_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4'])) \ No newline at end of file +controller.SaveEvents( + "IceCube_DIS.hdf5" +) From cf7d36e5c7c2c50c0dd2294ef118607efb950ad7 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:15:34 -0700 Subject: [PATCH 06/14] Add utilities for locating a good resources directory --- python/_util.py | 164 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 164 insertions(+) create mode 100644 python/_util.py diff --git a/python/_util.py b/python/_util.py new file mode 100644 index 000000000..3a3626521 --- /dev/null +++ b/python/_util.py @@ -0,0 +1,164 @@ +import os +import sys + +THIS_DIR = os.path.abspath(os.path.dirname(__file__)) + +# From pyzolib/paths.py (https://bitbucket.org/pyzo/pyzolib/src/tip/paths.py) +def appdata_dir(appname=None, roaming=False): + """appdata_dir(appname=None, roaming=False) + + Get the path to the application directory, where applications are allowed + to write user specific files (e.g. configurations). For non-user specific + data, consider using common_appdata_dir(). + If appname is given, a subdir is appended (and created if necessary). + If roaming is True, will prefer a roaming directory (Windows Vista/7). + """ + + # Define default user directory + userDir = os.getenv("LEPTONINJECTOR_USERDIR", None) + if userDir is None: + userDir = os.path.expanduser("~") + if not os.path.isdir(userDir): # pragma: no cover + userDir = "/var/tmp" # issue #54 + + # Get system app data dir + path = None + if sys.platform.startswith("win"): + path1, path2 = os.getenv("LOCALAPPDATA"), os.getenv("APPDATA") + path = (path2 or path1) if roaming else (path1 or path2) + elif sys.platform.startswith("darwin"): + path = os.path.join(userDir, "Library", "Application Support") + # On Linux and as fallback + if not (path and os.path.isdir(path)): + path = userDir + + # Maybe we should store things local to the executable (in case of a + # portable distro or a frozen application that wants to be portable) + prefix = sys.prefix + if getattr(sys, "frozen", None): + prefix = os.path.abspath(os.path.dirname(sys.executable)) + for reldir in ("settings", "../settings"): + localpath = os.path.abspath(os.path.join(prefix, reldir)) + if os.path.isdir(localpath): # pragma: no cover + try: + open(os.path.join(localpath, "test.write"), "wb").close() + os.remove(os.path.join(localpath, "test.write")) + except IOError: + pass # We cannot write in this directory + else: + path = localpath + break + + # Get path specific for this app + if appname: + if path == userDir: + appname = "." + appname.lstrip(".") # Make it a hidden directory + path = os.path.join(path, appname) + if not os.path.isdir(path): # pragma: no cover + os.makedirs(path, exist_ok=True) + + # Done + return path + +# from imageio +# https://github.com/imageio/imageio/blob/65d79140018bb7c64c0692ea72cb4093e8d632a0/imageio/core/util.py +def resource_dirs(): + """resource_dirs() + + Get a list of directories where leptoninjector resources may be located. + The first directory in this list is the "resources" directory in + the package itself. The second directory is the appdata directory + (~/.leptoninjector on Linux). The list further contains the application + directory (for frozen apps), and may include additional directories + in the future. + """ + dirs = [resource_package_dir()] + # Resource dir baked in the package. + # Appdata directory + try: + dirs.append(appdata_dir("leptoninjector")) + except Exception: # pragma: no cover + pass # The home dir may not be writable + # Directory where the app is located (mainly for frozen apps) + if getattr(sys, "frozen", None): + dirs.append(os.path.abspath(os.path.dirname(sys.executable))) + elif sys.path and sys.path[0]: + dirs.append(os.path.abspath(sys.path[0])) + return dirs + + +# from imageio +# https://github.com/imageio/imageio/blob/65d79140018bb7c64c0692ea72cb4093e8d632a0/imageio/core/util.py +def resource_package_dir(): + """package_dir + + Get the resources directory in the leptoninjector package installation + directory. + + Notes + ----- + This is a convenience method that is used by `resource_dirs` and + leptoninjector entry point scripts. + """ + # Make pkg_resources optional if setuptools is not available + try: + # Avoid importing pkg_resources in the top level due to how slow it is + # https://github.com/pypa/setuptools/issues/510 + import pkg_resources + except ImportError: + pkg_resources = None + + if pkg_resources: + # The directory returned by `pkg_resources.resource_filename` + # also works with eggs. + pdir = pkg_resources.resource_filename("leptoninjector", "resources") + else: + # If setuptools is not available, use fallback + pdir = os.path.abspath(os.path.join(THIS_DIR, "resources")) + return pdir + + +# from imageio +# https://github.com/imageio/imageio/blob/65d79140018bb7c64c0692ea72cb4093e8d632a0/imageio/core/util.py +def get_platform(): + """get_platform() + + Get a string that specifies the platform more specific than + sys.platform does. The result can be: linux32, linux64, win32, + win64, osx32, osx64. Other platforms may be added in the future. + """ + # Get platform + if sys.platform.startswith("linux"): + plat = "linux%i" + elif sys.platform.startswith("win"): + plat = "win%i" + elif sys.platform.startswith("darwin"): + plat = "osx%i" + elif sys.platform.startswith("freebsd"): + plat = "freebsd%i" + else: # pragma: no cover + return None + + return plat % (struct.calcsize("P") * 8) # 32 or 64 bits + + +# from imageio +# https://github.com/imageio/imageio/blob/65d79140018bb7c64c0692ea72cb4093e8d632a0/imageio/core/util.py +def has_module(module_name): + """Check to see if a python module is available.""" + if sys.version_info > (3, 4): + import importlib + + name_parts = module_name.split(".") + for i in range(len(name_parts)): + if importlib.util.find_spec(".".join(name_parts[: i + 1])) is None: + return False + return True + else: # pragma: no cover + import imp + + try: + imp.find_module(module_name) + except ImportError: + return False + return True From b003f5748e3ac1be23dca58adec88be3e38a2452 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:18:52 -0700 Subject: [PATCH 07/14] Make it black --- .../Examples/Example2/DipolePortal_MINERvA.py | 87 ++++++++++-------- .../Example2/DipolePortal_MiniBooNE.py | 88 ++++++++++--------- 2 files changed, 96 insertions(+), 79 deletions(-) diff --git a/resources/Examples/Example2/DipolePortal_MINERvA.py b/resources/Examples/Example2/DipolePortal_MINERvA.py index 7334fe668..6c1cff09b 100644 --- a/resources/Examples/Example2/DipolePortal_MINERvA.py +++ b/resources/Examples/Example2/DipolePortal_MINERvA.py @@ -2,80 +2,89 @@ import leptoninjector as LI import sys -sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python') + +sys.path.insert( + 1, + "/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python", +) from LIController import LIController - -LI_SRC = os.environ.get('LEPTONINJECTOR_SRC') -# Define a DarkNews model +LI_SRC = os.environ.get("LEPTONINJECTOR_SRC") + +# Define a DarkNews model model_kwargs = { - 'm4': 0.47,#0.140, - 'mu_tr_mu4': 2.5e-6, #1e-6, # GeV^-1 - 'UD4': 0, - 'Umu4': 0, - 'epsilon': 0.0, - 'gD': 0.0, - 'decay_product':'photon', - 'noHC':True, - 'HNLtype':"dirac", + "m4": 0.47, # 0.140, + "mu_tr_mu4": 2.5e-6, # 1e-6, # GeV^-1 + "UD4": 0, + "Umu4": 0, + "epsilon": 0.0, + "gD": 0.0, + "decay_product": "photon", + "noHC": True, + "HNLtype": "dirac", } # Number of events to inject events_to_inject = 1000 # Expeirment to run -experiment = 'MINERvA' +experiment = "MINERvA" # Define the controller -controller = LIController(events_to_inject, - experiment) +controller = LIController(events_to_inject, experiment) # Particle to inject primary_type = LI.dataclasses.Particle.ParticleType.NuMu # Define DarkNews Model -table_dir = LI_SRC+'/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']) -controller.InputDarkNewsModel(primary_type, - table_dir, - model_kwargs) +table_dir = ( + LI_SRC + + "/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/" + % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]) +) +controller.InputDarkNewsModel(primary_type, table_dir, model_kwargs) # Primary distributions primary_injection_distributions = {} primary_physical_distributions = {} # energy distribution -flux_file = LI_SRC + '/resources/Fluxes/NUMI_Flux_Tables/ME_FHC_numu.txt' +flux_file = LI_SRC + "/resources/Fluxes/NUMI_Flux_Tables/ME_FHC_numu.txt" edist = LI.distributions.TabulatedFluxDistribution(flux_file, True) -edist_gen = LI.distributions.TabulatedFluxDistribution(1.05*model_kwargs['m4'], 10, flux_file, False) -primary_injection_distributions['energy'] = edist_gen -primary_physical_distributions['energy'] = edist +edist_gen = LI.distributions.TabulatedFluxDistribution( + 1.05 * model_kwargs["m4"], 10, flux_file, False +) +primary_injection_distributions["energy"] = edist_gen +primary_physical_distributions["energy"] = edist # Flux normalization: go from cm^-2 to m^-2 flux_units = LI.distributions.NormalizationConstant(1e4) -primary_physical_distributions['normalization'] = flux_units +primary_physical_distributions["normalization"] = flux_units # direction distribution direction_distribution = LI.distributions.FixedDirection(LI.math.Vector3D(0, 0, 1.0)) -primary_injection_distributions['direction'] = direction_distribution -primary_physical_distributions['direction'] = direction_distribution +primary_injection_distributions["direction"] = direction_distribution +primary_physical_distributions["direction"] = direction_distribution # position distribution -decay_range_func = LI.distributions.DecayRangeFunction(model_kwargs['m4'], - controller.DN_min_decay_width, - 3, - 240) -position_distribution = LI.distributions.RangePositionDistribution(1.24, 5., - decay_range_func, - set(controller.GetDetectorModelTargets()[0])) -primary_injection_distributions['position'] = position_distribution +decay_range_func = LI.distributions.DecayRangeFunction( + model_kwargs["m4"], controller.DN_min_decay_width, 3, 240 +) +position_distribution = LI.distributions.RangePositionDistribution( + 1.24, 5.0, decay_range_func, set(controller.GetEarthModelTargets()[0]) +) +primary_injection_distributions["position"] = position_distribution # SetProcesses -controller.SetProcesses(primary_type, - primary_injection_distributions, - primary_physical_distributions) +controller.SetProcesses( + primary_type, primary_injection_distributions, primary_physical_distributions +) controller.Initialize() events = controller.GenerateEvents() -controller.SaveEvents('output/MINERvA_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4'])) \ No newline at end of file +controller.SaveEvents( + "output/MINERvA_Dipole_M%2.2f_mu%2.2e_example.hdf5" + % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]) +) diff --git a/resources/Examples/Example2/DipolePortal_MiniBooNE.py b/resources/Examples/Example2/DipolePortal_MiniBooNE.py index d13ede9b9..7cc207996 100644 --- a/resources/Examples/Example2/DipolePortal_MiniBooNE.py +++ b/resources/Examples/Example2/DipolePortal_MiniBooNE.py @@ -2,81 +2,89 @@ import leptoninjector as LI import sys -sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python') + +sys.path.insert( + 1, + "/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python", +) from LIController import LIController - -LI_SRC = os.environ.get('LEPTONINJECTOR_SRC') -# Define a DarkNews model +LI_SRC = os.environ.get("LEPTONINJECTOR_SRC") + +# Define a DarkNews model model_kwargs = { - 'm4': 0.47,#0.140, - 'mu_tr_mu4': 2.5e-6, #1e-6, # GeV^-1 - 'UD4': 0, - 'Umu4': 0, - 'epsilon': 0.0, - 'gD': 0.0, - 'decay_product':'photon', - 'noHC':True, - 'HNLtype':"dirac", + "m4": 0.47, # 0.140, + "mu_tr_mu4": 2.5e-6, # 1e-6, # GeV^-1 + "UD4": 0, + "Umu4": 0, + "epsilon": 0.0, + "gD": 0.0, + "decay_product": "photon", + "noHC": True, + "HNLtype": "dirac", } # Number of events to inject events_to_inject = 1000 # Expeirment to run -experiment = 'MiniBooNE' +experiment = "MiniBooNE" # Define the controller -controller = LIController(events_to_inject, - experiment) +controller = LIController(events_to_inject, experiment) # Particle to inject primary_type = LI.dataclasses.Particle.ParticleType.NuMu # Define DarkNews Model -table_dir = LI_SRC+'/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']) -controller.InputDarkNewsModel(primary_type, - table_dir, - model_kwargs) +table_dir = ( + LI_SRC + + "/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/" + % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]) +) +controller.InputDarkNewsModel(primary_type, table_dir, model_kwargs) # Primary distributions primary_injection_distributions = {} primary_physical_distributions = {} # energy distribution -flux_file = LI_SRC + '/resources/Fluxes/BNB_Flux_Tables/BNB_numu_flux.txt' +flux_file = LI_SRC + "/resources/Fluxes/BNB_Flux_Tables/BNB_numu_flux.txt" edist = LI.distributions.TabulatedFluxDistribution(flux_file, True) -edist_gen = LI.distributions.TabulatedFluxDistribution(1.05*model_kwargs['m4'], 10, flux_file, False) -primary_injection_distributions['energy'] = edist_gen -primary_physical_distributions['energy'] = edist +edist_gen = LI.distributions.TabulatedFluxDistribution( + 1.05 * model_kwargs["m4"], 10, flux_file, False +) +primary_injection_distributions["energy"] = edist_gen +primary_physical_distributions["energy"] = edist # Flux normalization: go from cm^-2 to m^-2 flux_units = LI.distributions.NormalizationConstant(1e4) -primary_physical_distributions['flux_units'] = flux_units +primary_physical_distributions["flux_units"] = flux_units # direction distribution direction_distribution = LI.distributions.FixedDirection(LI.math.Vector3D(0, 0, 1.0)) -primary_injection_distributions['direction'] = direction_distribution -primary_physical_distributions['direction'] = direction_distribution +primary_injection_distributions["direction"] = direction_distribution +primary_physical_distributions["direction"] = direction_distribution # position distribution -decay_range_func = LI.distributions.DecayRangeFunction(model_kwargs['m4'], - controller.DN_min_decay_width, - 3, - 541) -position_distribution = LI.distributions.RangePositionDistribution(6.2, 6.2, - decay_range_func, - set(controller.GetDetectorModelTargets()[0])) -primary_injection_distributions['position'] = position_distribution +decay_range_func = LI.distributions.DecayRangeFunction( + model_kwargs["m4"], controller.DN_min_decay_width, 3, 541 +) +position_distribution = LI.distributions.RangePositionDistribution( + 6.2, 6.2, decay_range_func, set(controller.GetEarthModelTargets()[0]) +) +primary_injection_distributions["position"] = position_distribution # SetProcesses -controller.SetProcesses(primary_type, - primary_injection_distributions, - primary_physical_distributions) +controller.SetProcesses( + primary_type, primary_injection_distributions, primary_physical_distributions +) controller.Initialize() events = controller.GenerateEvents() -controller.SaveEvents('output/MiniBooNE_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4'])) - +controller.SaveEvents( + "output/MiniBooNE_Dipole_M%2.2f_mu%2.2e_example.hdf5" + % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]) +) From 9f11908575bb742b0109916ce4186ef724391398 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:19:21 -0700 Subject: [PATCH 08/14] Re-order imports --- resources/Examples/Example1/DIS_IceCube.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/resources/Examples/Example1/DIS_IceCube.py b/resources/Examples/Example1/DIS_IceCube.py index 71b6b3d9d..d03962ff0 100644 --- a/resources/Examples/Example1/DIS_IceCube.py +++ b/resources/Examples/Example1/DIS_IceCube.py @@ -1,8 +1,8 @@ -import numpy as np import os +import sys +import numpy as np import leptoninjector as LI -import sys from leptoninjector import _util from leptoninjector.LIController import LIController From 62402b542592c3e57c9eecfa6e5d10bf7c4be292 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:41:58 -0700 Subject: [PATCH 09/14] Make it black. --- python/LIDarkNews.py | 703 +++++++++++++++++++++++++------------------ 1 file changed, 417 insertions(+), 286 deletions(-) diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py index cb7012bf4..826cb8845 100644 --- a/python/LIDarkNews.py +++ b/python/LIDarkNews.py @@ -5,11 +5,11 @@ import json import ntpath import pickle -from scipy.interpolate import CloughTocher2DInterpolator,CubicSpline +from scipy.interpolate import CloughTocher2DInterpolator, CubicSpline # LeptonInjector methods import leptoninjector as LI -from leptoninjector.interactions import DarkNewsCrossSection,DarkNewsDecay +from leptoninjector.interactions import DarkNewsCrossSection, DarkNewsDecay from leptoninjector.dataclasses import Particle # DarkNews methods @@ -19,110 +19,122 @@ from DarkNews.nuclear_tools import NuclearTarget from DarkNews.integrands import get_decay_momenta_from_vegas_samples + # Class containing all upscattering and decay modes available in DarkNews class PyDarkNewsInteractionCollection: - - def __init__(self, - table_dir=None, - param_file=None, - tolerance=1e-6, - interp_tolerance=5e-2, - **kwargs): + def __init__( + self, + table_dir=None, + param_file=None, + tolerance=1e-6, + interp_tolerance=5e-2, + **kwargs, + ): # Defines a series of upscattering and decay objects # Each derive from the respective LeptonInjector classes - + # Get our model container with all ups_case and dec_case DarkNews objects - self.models = ModelContainer(param_file,**kwargs) + self.models = ModelContainer(param_file, **kwargs) self.table_dir = table_dir - + # Default table_dir settings if self.table_dir is None: - self.table_dir = os.environ.get('LEPTONINJECTOR_SRC') + '/resources/CrossSectionTables/DarkNewsTables/' - ct = datetime.datetime.now().strftime('%Y_%m_%d__%H:%M/') + self.table_dir = ( + os.environ.get("LEPTONINJECTOR_SRC") + + "/resources/CrossSectionTables/DarkNewsTables/" + ) + ct = datetime.datetime.now().strftime("%Y_%m_%d__%H:%M/") self.table_dir += ct - + # Make the table directory where will we store cross section integrators table_dir_exists = False - if(os.path.exists(self.table_dir)): - #print("Directory '%s' already exists"%self.table_dir) + if os.path.exists(self.table_dir): + # print("Directory '%s' already exists"%self.table_dir) table_dir_exists = True else: try: - os.makedirs(self.table_dir, exist_ok = False) - print("Directory '%s' created successfully"%self.table_dir) + os.makedirs(self.table_dir, exist_ok=False) + print("Directory '%s' created successfully" % self.table_dir) except OSError as error: - print("Directory '%s' cannot be created"%self.table_dir) + print("Directory '%s' cannot be created" % self.table_dir) exit(0) if table_dir_exists: # Ensure that the model requested matches the model file already in the dictionary if param_file is not None: # ensure the param filename already exists - param_filename = ntpath.basename(param_file) # should be OS-independent - assert(os.path.isfile(self.table_dir + param_filename)) + param_filename = ntpath.basename(param_file) # should be OS-independent + assert os.path.isfile(self.table_dir + param_filename) # Make sure the model arguments agree - with open(self.table_dir+'model_parameters.json',) as f: + with open( + self.table_dir + "model_parameters.json", + ) as f: _model_args_dict = json.load(f) - assert(self.models.model_args_dict==_model_args_dict) + assert self.models.model_args_dict == _model_args_dict else: # Write a file to the directory containing infomration on the parameters used to create the model if param_file is not None: # Copy the param_file to the folder - command = 'scp ' + param_file + ' ' + self.table_dir + command = "scp " + param_file + " " + self.table_dir os.system(command) # Dump the model arguments - with open(self.table_dir+'model_parameters.json','w') as f: - json.dump(self.models.model_args_dict,f) - + with open(self.table_dir + "model_parameters.json", "w") as f: + json.dump(self.models.model_args_dict, f) + # Save all unique scattering processes self.cross_sections = [] - for ups_key,ups_case in self.models.ups_cases.items(): - table_subdirs = 'CrossSection_' + for ups_key, ups_case in self.models.ups_cases.items(): + table_subdirs = "CrossSection_" for x in ups_key: - if(type(x)==NuclearTarget): + if type(x) == NuclearTarget: x = x.name - table_subdirs += '%s_'%str(x) - table_subdirs += '/' - self.cross_sections.append(PyDarkNewsCrossSection(ups_case, - table_dir = self.table_dir + table_subdirs, - tolerance=tolerance, - interp_tolerance=interp_tolerance)) - # Save all unique decay processes + table_subdirs += "%s_" % str(x) + table_subdirs += "/" + self.cross_sections.append( + PyDarkNewsCrossSection( + ups_case, + table_dir=self.table_dir + table_subdirs, + tolerance=tolerance, + interp_tolerance=interp_tolerance, + ) + ) + # Save all unique decay processes self.decays = [] - for dec_key,dec_case in self.models.dec_cases.items(): - table_subdirs = 'Decay_' + for dec_key, dec_case in self.models.dec_cases.items(): + table_subdirs = "Decay_" for x in dec_key: - table_subdirs += '%s_'%str(x) - table_subdirs += '/' - self.decays.append(PyDarkNewsDecay(dec_case, - table_dir = self.table_dir + table_subdirs)) - - def SaveCrossSectionTables(self,fill_tables_at_exit=False): + table_subdirs += "%s_" % str(x) + table_subdirs += "/" + self.decays.append( + PyDarkNewsDecay(dec_case, table_dir=self.table_dir + table_subdirs) + ) + + def SaveCrossSectionTables(self, fill_tables_at_exit=True): if not fill_tables_at_exit: - print('WARNING: Saving tables without filling PyDarkNewsCrossSection interpolation tables. Future updates to DarkNews can lead to inconsistent behavior if new entries are ever added to this table') + print( + "WARNING: Saving tables without filling PyDarkNewsCrossSection interpolation tables. Future updates to DarkNews can lead to inconsistent behavior if new entries are ever added to this table" + ) for cross_section in self.cross_sections: if fill_tables_at_exit: - print('Filling cross section table at %s'%cross_section.table_dir) + print("Filling cross section table at %s" % cross_section.table_dir) num = cross_section.FillInterpolationTables() - print('Added %d points'%num) + print("Added %d points" % num) cross_section.SaveInterpolationTables() - - # A class representing a single ups_case DarkNews class # Only handles methods concerning the upscattering part class PyDarkNewsCrossSection(DarkNewsCrossSection): + def __init__( + self, + ups_case, # DarkNews UpscatteringProcess instance + table_dir=None, # table to store + tolerance=1e-6, # supposed to represent machine epsilon + interp_tolerance=5e-2, # relative interpolation tolerance + interpolate_differential=False, + ): + DarkNewsCrossSection.__init__(self) # C++ constructor - def __init__(self, - ups_case, # DarkNews UpscatteringProcess instance - table_dir=None, # table to store - tolerance=1e-6, # supposed to represent machine epsilon - interp_tolerance=5e-2, # relative interpolation tolerance - interpolate_differential=False): - - DarkNewsCrossSection.__init__(self) # C++ constructor - self.ups_case = ups_case self.tolerance = tolerance self.interp_tolerance = interp_tolerance @@ -130,58 +142,71 @@ def __init__(self, self.interpolate_differential = interpolate_differential # 2D table in E, sigma - self.total_cross_section_table = np.empty((0,2),dtype=float) + self.total_cross_section_table = np.empty((0, 2), dtype=float) self.total_cross_section_interpolator = None # 3D table in E, z, dsigma/dQ2 where z = (Q2 - Q2min) / (Q2max - Q2min) - self.differential_cross_section_table = np.empty((0,3),dtype=float) + self.differential_cross_section_table = np.empty((0, 3), dtype=float) self.differential_cross_section_interpolator = None - - if table_dir is None: - print('No table_dir specified; disabling interpolation\nWARNING: this will siginficantly slow down event generation') + + if table_dir is None: + print( + "No table_dir specified; disabling interpolation\nWARNING: this will siginficantly slow down event generation" + ) return # Make the table directory where will we store cross section integrators table_dir_exists = False - if(os.path.exists(self.table_dir)): - #print("Directory '%s' already exists"%self.table_dir) + if os.path.exists(self.table_dir): + # print("Directory '%s' already exists"%self.table_dir) table_dir_exists = True else: try: - os.makedirs(self.table_dir, exist_ok = False) - print("Directory '%s' created successfully"%self.table_dir) + os.makedirs(self.table_dir, exist_ok=False) + print("Directory '%s' created successfully" % self.table_dir) except OSError as error: - print("Directory '%s' cannot be created"%self.table_dir) + print("Directory '%s' cannot be created" % self.table_dir) exit(0) - + # Look in table dir and check whether total/differential xsec tables exist if table_dir_exists: - total_xsec_file = self.table_dir + 'total_cross_sections.npy' + total_xsec_file = 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 = self.table_dir + 'differential_cross_sections.npy' + diff_xsec_file = self.table_dir + "differential_cross_sections.npy" if os.path.exists(diff_xsec_file): self.differential_cross_section_table = np.load(diff_xsec_file) - - self._redefine_interpolation_objects(total=True,diff=True) + + self._redefine_interpolation_objects(total=True, diff=True) # Sorts and redefines scipy interpolation objects - def _redefine_interpolation_objects(self,total=False,diff=False): + 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) + self.total_cross_section_table = np.sort( + self.total_cross_section_table, axis=0 + ) 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]) + 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) + self.differential_cross_section_table = np.sort( + self.differential_cross_section_table, axis=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) + 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): + def _interpolation_flags(self, inputs, mode): # # returns UseSinglePoint,Interpolate,closest_idx # UseSinglePoint: whether to use a single point in table @@ -189,141 +214,152 @@ def _interpolation_flags(self,inputs,mode): # closest_idx: index of closest point in table (for UseSinglePoint) # Determine which table we are using - if mode=='total': + if mode == "total": interp_table = self.total_cross_section_table - elif mode=='differential': + elif mode == "differential": interp_table = self.differential_cross_section_table else: - print('Invalid interpolation table mode %s'%mode) + print("Invalid interpolation table mode %s" % mode) exit(0) # first check if we have saved table points already - if len(interp_table) == 0: - return False,False,-1 + if len(interp_table) == 0: + return False, False, -1 # bools to keep track of whether to use a single point or interpolate UseSinglePoint = True Interpolate = True # First check whether we have a close-enough single point - closest_idx = np.argmin(np.sum(np.abs(interp_table[:,:-1] - inputs))) - diff = (interp_table[closest_idx,:-1] - inputs)/inputs - if np.all(np.abs(diff)= self.interp_tolerance: Interpolate = False else: # closest existing input is within interpolation range - if diff>=0: + if diff >= 0: # closest existing input is above or equal to requested input idx_above = closest_idx # check if we are at the boundary if closest_idx == 0: Interpolate = False else: - idx_below = closest_idx-1 - diff_below = (input - interp_table[idx_below,i])/input + idx_below = closest_idx - 1 + diff_below = (input - interp_table[idx_below, i]) / input # check if the node below is also within the interpolation tolerance - if diff_below>=self.interp_tolerance: - Interpolate = False - elif diff<0 and -diff= self.interp_tolerance: + Interpolate = False + elif diff < 0 and -diff < self.interp_tolerance: # closest existing input is below requested input idx_below = closest_idx # check if we are at boundary - if closest_idx >= len(interp_table)-1: + if closest_idx >= len(interp_table) - 1: Interpolate = False else: - idx_above = closest_idx+1 - diff_above = (interp_table[idx_above,i] - input)/input + idx_above = closest_idx + 1 + diff_above = (interp_table[idx_above, i] - input) / input # check if the node above is also within the interpolation tolerance - if diff_above>=self.interp_tolerance: + if diff_above >= self.interp_tolerance: Interpolate = False return UseSinglePoint, Interpolate, closest_idx - - # return entries in interpolation table if we have inputs - def _query_interpolation_table(self,inputs,mode): - # + + # return entries in interpolation table if we have inputs + def _query_interpolation_table(self, inputs, mode): + # # returns: # 0 if we are not close enough to any points in the interpolation table # otherwise, returns the desired interpolated value - + # Determine which table we are using - if mode=='total': + if mode == "total": interp_table = self.total_cross_section_table interpolator = self.total_cross_section_interpolator - elif mode=='differential': + elif mode == "differential": interp_table = self.differential_cross_section_table interpolator = self.differential_cross_section_interpolator else: - print('Invalid interpolation table mode %s'%mode) + print("Invalid interpolation table mode %s" % mode) exit(0) - UseSinglePoint, Interpolate, closest_idx = self._interpolation_flags(inputs,mode) + UseSinglePoint, Interpolate, closest_idx = self._interpolation_flags( + inputs, mode + ) if UseSinglePoint: - return interp_table[closest_idx,-1] + return interp_table[closest_idx, -1] elif Interpolate: return interpolator(inputs) else: - return 0 - + return 0 + # Fills the total and differential cross section tables within interp_tolerance - def FillInterpolationTables(self,total=True,diff=True): - Emin = (1.0+self.tolerance)*self.ups_case.Ethreshold - Emax = np.max(self.total_cross_section_table[:,0]) + def FillInterpolationTables(self, total=True, diff=True): + Emin = (1.0 + self.tolerance) * self.ups_case.Ethreshold + Emax = np.max(self.total_cross_section_table[:, 0]) num_added_points = 0 if total: E = Emin - while E 0: # we have recovered the differential cross section from the interpolation table return val - + # If we have reached this block, we must compute the differential cross section using DarkNews dxsec = self.ups_case.diff_xsec_Q2(energy, Q2).item() - self.differential_cross_section_table = np.append(self.differential_cross_section_table,[[energy,z,dxsec]],axis=0) - if self.interpolate_differential: self._redefine_interpolation_objects(diff=True) + self.differential_cross_section_table = np.append( + self.differential_cross_section_table, [[energy, z, dxsec]], axis=0 + ) + if self.interpolate_differential: + self._redefine_interpolation_objects(diff=True) return dxsec - + def SetUpscatteringMasses(self, interaction): interaction.primary_mass = 0 interaction.target_mass = self.ups_case.MA @@ -417,22 +488,24 @@ def SetUpscatteringMasses(self, interaction): def SetUpscatteringHelicities(self, interaction): secondary_helicities = [] - secondary_helicities.append(self.ups_case.h_upscattered * interaction.primary_helicity) + secondary_helicities.append( + self.ups_case.h_upscattered * interaction.primary_helicity + ) secondary_helicities.append(interaction.target_helicity) interaction.secondary_helicity = secondary_helicities self.h_ups = self.ups_case.m_ups self.h_target = self.ups_case.MA - + def TotalCrossSection(self, arg1, energy=None, target=None): # Handle overloaded arguments - if type(arg1==LI.dataclasses.InteractionRecord): + if type(arg1 == LI.dataclasses.InteractionRecord): primary = arg1.signature.primary_type energy = arg1.primary_momentum[0] target = arg1.signature.target_type elif energy is not None and target is not None: primary = arg1 else: - print('Incorrect function call to TotalCrossSection!') + print("Incorrect function call to TotalCrossSection!") exit(0) if primary != self.ups_case.nu_projectile: return 0 @@ -442,40 +515,47 @@ def TotalCrossSection(self, arg1, energy=None, target=None): interaction.primary_momentum[0] = energy if energy < self.InteractionThreshold(interaction): return 0 - + # Check if we can interpolate - val = self._query_interpolation_table([energy],mode='total') + val = self._query_interpolation_table([energy], mode="total") if val > 0: # we have recovered the cross section from the interpolation table return val - + # If we have reached this block, we must compute the cross section using DarkNews xsec = self.ups_case.scalar_total_xsec(energy) - self.total_cross_section_table = np.append(self.total_cross_section_table,[[energy,xsec]],axis=0) + self.total_cross_section_table = np.append( + self.total_cross_section_table, [[energy, xsec]], axis=0 + ) self._redefine_interpolation_objects(total=True) return xsec - def InteractionThreshold(self, interaction): return self.ups_case.Ethreshold def Q2Min(self, interaction): - return phase_space.upscattering_Q2min(interaction.primary_momentum[0], self.ups_case.m_ups, interaction.target_mass) + return phase_space.upscattering_Q2min( + interaction.primary_momentum[0], + self.ups_case.m_ups, + interaction.target_mass, + ) def Q2Max(self, interaction): - return phase_space.upscattering_Q2max(interaction.primary_momentum[0], self.ups_case.m_ups, interaction.target_mass) + return phase_space.upscattering_Q2max( + interaction.primary_momentum[0], + self.ups_case.m_ups, + interaction.target_mass, + ) + # A class representing a single decay_case DarkNews class # Only handles methods concerning the decay part class PyDarkNewsDecay(DarkNewsDecay): - - def __init__(self, - dec_case, - table_dir=None): - DarkNewsDecay.__init__(self) # C++ constructor + def __init__(self, dec_case, table_dir=None): + DarkNewsDecay.__init__(self) # C++ constructor self.dec_case = dec_case self.table_dir = table_dir - + # Some variables for storing the decay phase space integrator self.decay_integrator = None self.decay_norm = None @@ -484,39 +564,43 @@ def __init__(self, self.PS_weights_CDF = None self.total_width = None - if table_dir is None: - print('No table_dir specified; will sample from new VEGAS integrator for each decay') - print('WARNING: this will siginficantly slow down event generation') + if table_dir is None: + print( + "No table_dir specified; will sample from new VEGAS integrator for each decay" + ) + print("WARNING: this will siginficantly slow down event generation") return # Make the table directory where will we store cross section integrators table_dir_exists = False - if(os.path.exists(self.table_dir)): - #print("Directory '%s' already exists"%self.table_dir) + if os.path.exists(self.table_dir): + # print("Directory '%s' already exists"%self.table_dir) table_dir_exists = True else: try: - os.makedirs(self.table_dir, exist_ok = False) - print("Directory '%s' created successfully"%self.table_dir) + os.makedirs(self.table_dir, exist_ok=False) + print("Directory '%s' created successfully" % self.table_dir) except OSError as error: - print("Directory '%s' cannot be created"%self.table_dir) + print("Directory '%s' cannot be created" % self.table_dir) exit(0) - + # Look in table dir for existing decay integrator + normalization info if table_dir_exists: self.SetIntegratorAndNorm() def SetIntegratorAndNorm(self): # Try to find the decay integrator - int_file = self.table_dir + 'decay_integrator.pkl' + int_file = self.table_dir + "decay_integrator.pkl" if os.path.isfile(int_file): - with open(int_file,'rb') as ifile: - _,self.decay_integrator = pickle.load(ifile) + with open(int_file, "rb") as ifile: + _, self.decay_integrator = pickle.load(ifile) # Try to find the normalization information - norm_file = self.table_dir + 'decay_norm.json' + norm_file = self.table_dir + "decay_norm.json" if os.path.isfile(norm_file): - with open(norm_file,) as nfile: - self.decay_norm = json.load(nfile) + with open( + norm_file, + ) as nfile: + self.decay_norm = json.load(nfile) ##### START METHODS FOR SERIALIZATION ######### # def get_initialized_dict(config): @@ -525,7 +609,7 @@ def SetIntegratorAndNorm(self): # return pddn.__dict__ # # return the conent of __dict__ for PyDerivedDarkNews - # @staticmethod + # @staticmethod # def get_config(self): # return self.config ##### END METHODS FOR SERIALIZATION ######### @@ -542,12 +626,16 @@ def GetPossibleSignatures(self): return [signature] def GetPossibleSignaturesFromParent(self, primary_type): - if (Particle.ParticleType(self.dec_case.nu_parent.pdgid) == primary_type): + if Particle.ParticleType(self.dec_case.nu_parent.pdgid) == primary_type: signature = LI.dataclasses.InteractionSignature() - signature.primary_type = Particle.ParticleType(self.dec_case.nu_parent.pdgid) + signature.primary_type = Particle.ParticleType( + self.dec_case.nu_parent.pdgid + ) signature.target_type = Particle.ParticleType.Decay secondary_types = [] - secondary_types.append(Particle.ParticleType(self.dec_case.nu_daughter.pdgid)) + secondary_types.append( + Particle.ParticleType(self.dec_case.nu_daughter.pdgid) + ) for secondary in self.dec_case.secondaries: secondary_types.append(Particle.ParticleType(secondary.pdgid)) signature.secondary_types = secondary_types @@ -555,178 +643,221 @@ def GetPossibleSignaturesFromParent(self, primary_type): return [] def DifferentialDecayWidth(self, record): - # Momentum variables of HNL necessary for calculating decay phase space PN = np.array(record.primary_momentum) - - if type(self.dec_case)==FermionSinglePhotonDecay: + + if type(self.dec_case) == FermionSinglePhotonDecay: gamma_idx = 0 for secondary in record.signature.secondary_types: if secondary == LI.dataclasses.Particle.ParticleType.Gamma: break gamma_idx += 1 if gamma_idx >= len(record.signature.secondary_types): - print('No gamma found in the list of secondaries!') + print("No gamma found in the list of secondaries!") exit(0) - + Pgamma = np.array(record.secondary_momenta[gamma_idx]) - momenta = np.expand_dims(PN,0),np.expand_dims(Pgamma,0) + momenta = np.expand_dims(PN, 0), np.expand_dims(Pgamma, 0) - elif type(self.dec_case)==FermionDileptonDecay: + elif type(self.dec_case) == FermionDileptonDecay: lepminus_idx = -1 lepplus_idx = -1 nu_idx = -1 - for idx,secondary in enumerate(record.signature.secondary_types): - if secondary in [LI.dataclasses.Particle.ParticleType.EMinus, - LI.dataclasses.Particle.ParticleType.MuMinus, - LI.dataclasses.Particle.ParticleType.TauMinus]: + for idx, secondary in enumerate(record.signature.secondary_types): + if secondary in [ + LI.dataclasses.Particle.ParticleType.EMinus, + LI.dataclasses.Particle.ParticleType.MuMinus, + LI.dataclasses.Particle.ParticleType.TauMinus, + ]: lepminus_idx = idx - elif secondary in [LI.dataclasses.Particle.ParticleType.EPlus, - LI.dataclasses.Particle.ParticleType.MuPlus, - LI.dataclasses.Particle.ParticleType.TauPlus]: + elif secondary in [ + LI.dataclasses.Particle.ParticleType.EPlus, + LI.dataclasses.Particle.ParticleType.MuPlus, + LI.dataclasses.Particle.ParticleType.TauPlus, + ]: lepplus_idx = idx else: nu_idx = idx - if -1 in [lepminus_idx,lepplus_idx,nu_idx]: - print('Couldn\'t find two leptons and a neutrino in the final state!') + if -1 in [lepminus_idx, lepplus_idx, nu_idx]: + print("Couldn't find two leptons and a neutrino in the final state!") exit(0) Pnu = np.array(record.secondary_momenta[nu_idx]) Plepminus = np.array(record.secondary_momenta[lepminus_idx]) Plepplus = np.array(record.secondary_momenta[lepplus_idx]) - momenta = np.expand_dims(PN,0),np.expand_dims(Plepminus,0),np.expand_dims(Plepplus,0),np.expand_dims(Pnu,0) + momenta = ( + np.expand_dims(PN, 0), + np.expand_dims(Plepminus, 0), + np.expand_dims(Plepplus, 0), + np.expand_dims(Pnu, 0), + ) else: - print('%s is not a valid decay class type!'%type(self.dec_case)) + print("%s is not a valid decay class type!" % type(self.dec_case)) exit(0) return self.dec_case.differential_width(momenta) - + def TotalDecayWidth(self, arg1): - if type(arg1)==LI.dataclasses.InteractionRecord: + if type(arg1) == LI.dataclasses.InteractionRecord: primary = arg1.signature.primary_type - elif type(arg1)==LI.dataclasses.Particle.ParticleType: + elif type(arg1) == LI.dataclasses.Particle.ParticleType: primary = arg1 else: - print('Incorrect function call to TotalDecayWidth!') + print("Incorrect function call to TotalDecayWidth!") exit(0) if primary != self.dec_case.nu_parent: return 0 if self.total_width is None: # Need to set the total width - if type(self.dec_case) == FermionDileptonDecay and \ - (self.dec_case.vector_off_shell and self.dec_case.scalar_off_shell): + if type(self.dec_case) == FermionDileptonDecay and ( + self.dec_case.vector_off_shell and self.dec_case.scalar_off_shell + ): # total width calculation requires evaluating an integral - if (self.decay_integrator is None or self.decay_norm is None): - # We need to initialize a new VEGAS integrator in DarkNews - int_file = self.table_dir + 'decay_integrator.pkl' - norm_file = self.table_dir + 'decay_norm.json' - self.total_width = self.dec_case.total_width(savefile_norm=norm_file,savefile_dec=int_file) - self.SetIntegratorAndNorm() + if self.decay_integrator is None or self.decay_norm is None: + # We need to initialize a new VEGAS integrator in DarkNews + int_file = self.table_dir + "decay_integrator.pkl" + norm_file = self.table_dir + "decay_norm.json" + self.total_width = self.dec_case.total_width( + savefile_norm=norm_file, savefile_dec=int_file + ) + self.SetIntegratorAndNorm() else: - self.total_width = self.decay_integrator["diff_decay_rate_0"].mean * self.decay_norm["diff_decay_rate_0"] + self.total_width = ( + self.decay_integrator["diff_decay_rate_0"].mean + * self.decay_norm["diff_decay_rate_0"] + ) else: self.total_width = self.dec_case.total_width() return self.total_width - - - def TotalDecayWidthForFinalState(self,record): + + def TotalDecayWidthForFinalState(self, record): sig = self.GetPossibleSignatures()[0] - if (record.signature.primary_type != sig.primary_type) or \ - (record.signature.target_type != sig.target_type) or \ - (len(record.signature.secondary_types) != len(sig.secondary_types)) or \ - (np.any([record.signature.secondary_types[i] != sig.secondary_types[i] for i in range(len(sig.secondary_types))])): + if ( + (record.signature.primary_type != sig.primary_type) + or (record.signature.target_type != sig.target_type) + or (len(record.signature.secondary_types) != len(sig.secondary_types)) + or ( + np.any( + [ + record.signature.secondary_types[i] != sig.secondary_types[i] + for i in range(len(sig.secondary_types)) + ] + ) + ) + ): return 0 return self.dec_case.total_width() - + def DensityVariables(self): - if type(self.dec_case)==FermionSinglePhotonDecay: + if type(self.dec_case) == FermionSinglePhotonDecay: return "cost" - elif type(self.dec_case)==FermionDileptonDecay: + elif type(self.dec_case) == FermionDileptonDecay: if self.dec_case.vector_on_shell and self.dec_case.scalar_on_shell: - print('Can\'t have both the scalar and vector on shell') + print("Can't have both the scalar and vector on shell") exit(0) - elif (self.dec_case.vector_on_shell and self.dec_case.scalar_off_shell) or \ - (self.dec_case.vector_off_shell and self.dec_case.scalar_on_shell): + elif (self.dec_case.vector_on_shell and self.dec_case.scalar_off_shell) or ( + self.dec_case.vector_off_shell and self.dec_case.scalar_on_shell + ): return "cost" elif self.dec_case.vector_off_shell and self.dec_case.scalar_off_shell: return "t,u,c3,phi34" else: - print('%s is not a valid decay class type!'%type(self.dec_case)) + print("%s is not a valid decay class type!" % type(self.dec_case)) exit(0) return "" - + def GetPSSample(self, random): - # Make the PS weight CDF if that hasn't been done if self.PS_weights_CDF is None: self.PS_weights_CDF = np.cumsum(self.PS_weights) - - # Random number to determine - x = random.Uniform(0,self.PS_weights_CDF[-1]) + + # Random number to determine + x = random.Uniform(0, self.PS_weights_CDF[-1]) # find first instance of a CDF entry greater than x PSidx = np.argmax(x - self.PS_weights_CDF <= 0) - return self.PS_samples[:,PSidx] + return self.PS_samples[:, PSidx] - def SampleRecordFromDarkNews(self,record,random): - + def SampleRecordFromDarkNews(self, record, random): # First, make sure we have PS samples and weights if self.PS_samples is None or self.PS_weights is None: # We need to generate new PS samples if self.decay_integrator is None or self.decay_norm is None: # We need to initialize a new VEGAS integrator in DarkNews - int_file = self.table_dir + 'decay_integrator.pkl' - norm_file = self.table_dir + 'decay_norm.json' - self.PS_samples, PS_weights_dict = self.dec_case.SamplePS(savefile_norm=norm_file,savefile_dec=int_file) - self.PS_weights = PS_weights_dict['diff_decay_rate_0'] + int_file = self.table_dir + "decay_integrator.pkl" + norm_file = self.table_dir + "decay_norm.json" + self.PS_samples, PS_weights_dict = self.dec_case.SamplePS( + savefile_norm=norm_file, savefile_dec=int_file + ) + self.PS_weights = PS_weights_dict["diff_decay_rate_0"] self.SetIntegratorAndNorm() else: # We already have an integrator, we just need new PS samples - self.PS_samples, PS_weights_dict = self.dec_case.SamplePS(existing_integrator=self.decay_integrator) - self.PS_weights = PS_weights_dict['diff_decay_rate_0'] - + self.PS_samples, PS_weights_dict = self.dec_case.SamplePS( + existing_integrator=self.decay_integrator + ) + self.PS_weights = PS_weights_dict["diff_decay_rate_0"] + # Now we must sample an PS point on the hypercube PS = self.GetPSSample(random) # Find the four-momenta associated with this point # Expand dims required to call DarkNews function on signle sample - four_momenta = get_decay_momenta_from_vegas_samples(np.expand_dims(PS,0),self.dec_case,np.expand_dims(np.array(record.primary_momentum),0)) - - if type(self.dec_case)==FermionSinglePhotonDecay: + four_momenta = get_decay_momenta_from_vegas_samples( + np.expand_dims(PS, 0), + self.dec_case, + np.expand_dims(np.array(record.primary_momentum), 0), + ) + + if type(self.dec_case) == FermionSinglePhotonDecay: gamma_idx = 0 for secondary in record.signature.secondary_types: if secondary == LI.dataclasses.Particle.ParticleType.Gamma: break gamma_idx += 1 if gamma_idx >= len(record.signature.secondary_types): - print('No gamma found in the list of secondaries!') + print("No gamma found in the list of secondaries!") exit(0) nu_idx = 1 - gamma_idx secondary_momenta = [] - secondary_momenta.insert(gamma_idx,list(np.squeeze(four_momenta["P_decay_photon"]))) - secondary_momenta.insert(nu_idx,list(np.squeeze(four_momenta["P_decay_N_daughter"]))) + secondary_momenta.insert( + gamma_idx, list(np.squeeze(four_momenta["P_decay_photon"])) + ) + secondary_momenta.insert( + nu_idx, list(np.squeeze(four_momenta["P_decay_N_daughter"])) + ) record.secondary_momenta = secondary_momenta - elif type(self.dec_case)==FermionDileptonDecay: + elif type(self.dec_case) == FermionDileptonDecay: lepminus_idx = -1 lepplus_idx = -1 nu_idx = -1 - for idx,secondary in enumerate(record.signature.secondary_types): - if secondary in [LI.dataclasses.Particle.ParticleType.EMinus, - LI.dataclasses.Particle.ParticleType.MuMinus, - LI.dataclasses.Particle.ParticleType.TauMinus]: + for idx, secondary in enumerate(record.signature.secondary_types): + if secondary in [ + LI.dataclasses.Particle.ParticleType.EMinus, + LI.dataclasses.Particle.ParticleType.MuMinus, + LI.dataclasses.Particle.ParticleType.TauMinus, + ]: lepminus_idx = idx - elif secondary in [LI.dataclasses.Particle.ParticleType.EPlus, - LI.dataclasses.Particle.ParticleType.MuPlus, - LI.dataclasses.Particle.ParticleType.TauPlus]: + elif secondary in [ + LI.dataclasses.Particle.ParticleType.EPlus, + LI.dataclasses.Particle.ParticleType.MuPlus, + LI.dataclasses.Particle.ParticleType.TauPlus, + ]: lepplus_idx = idx else: nu_idx = idx - if -1 in [lepminus_idx,lepplus_idx,nu_idx]: - print([lepminus_idx,lepplus_idx,nu_idx]) + if -1 in [lepminus_idx, lepplus_idx, nu_idx]: + print([lepminus_idx, lepplus_idx, nu_idx]) print(record.signature.secondary_types) - print('Couldn\'t find two leptons and a neutrino in the final state!') + print("Couldn't find two leptons and a neutrino in the final state!") exit(0) secondary_momenta = [] - secondary_momenta.insert(lepminus_idx,list(np.squeeze(four_momenta["P_decay_ell_minus"]))) - secondary_momenta.insert(lepplus_idx,list(np.squeeze(four_momenta["P_decay_ell_plus"]))) - secondary_momenta.insert(nu_idx,list(np.squeeze(four_momenta["P_decay_N_daughter"]))) + secondary_momenta.insert( + lepminus_idx, list(np.squeeze(four_momenta["P_decay_ell_minus"])) + ) + secondary_momenta.insert( + lepplus_idx, list(np.squeeze(four_momenta["P_decay_ell_plus"])) + ) + secondary_momenta.insert( + nu_idx, list(np.squeeze(four_momenta["P_decay_N_daughter"])) + ) record.secondary_momenta = secondary_momenta return record From 018df0d92d4b864a21bcddd267f20397babd8e7c Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:42:34 -0700 Subject: [PATCH 10/14] Use relative imports and the resources_dir --- .../Examples/Example2/DipolePortal_MINERvA.py | 31 ++++++++++--------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/resources/Examples/Example2/DipolePortal_MINERvA.py b/resources/Examples/Example2/DipolePortal_MINERvA.py index 6c1cff09b..95bab0902 100644 --- a/resources/Examples/Example2/DipolePortal_MINERvA.py +++ b/resources/Examples/Example2/DipolePortal_MINERvA.py @@ -1,15 +1,12 @@ import os - -import leptoninjector as LI import sys +import numpy as np -sys.path.insert( - 1, - "/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python", -) -from LIController import LIController +import leptoninjector as LI +from leptoninjector import _util +from leptoninjector.LIController import LIController -LI_SRC = os.environ.get("LEPTONINJECTOR_SRC") +resources_dir = _util.resource_package_dir() # Define a DarkNews model model_kwargs = { @@ -37,10 +34,11 @@ primary_type = LI.dataclasses.Particle.ParticleType.NuMu # Define DarkNews Model -table_dir = ( - LI_SRC - + "/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/" - % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]) +table_dir = os.path.join( + resources_dir, + "CrossSectionTables", + "DarkNewsTables", + "Dipole_M%2.2f_mu%2.2e" % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]), ) controller.InputDarkNewsModel(primary_type, table_dir, model_kwargs) @@ -49,7 +47,12 @@ primary_physical_distributions = {} # energy distribution -flux_file = LI_SRC + "/resources/Fluxes/NUMI_Flux_Tables/ME_FHC_numu.txt" +flux_file = os.path.join( + resources_dir, + "Fluxes", + "NUMI_Flux_Tables", + "ME_FHC_numu.txt", +) edist = LI.distributions.TabulatedFluxDistribution(flux_file, True) edist_gen = LI.distributions.TabulatedFluxDistribution( 1.05 * model_kwargs["m4"], 10, flux_file, False @@ -71,7 +74,7 @@ model_kwargs["m4"], controller.DN_min_decay_width, 3, 240 ) position_distribution = LI.distributions.RangePositionDistribution( - 1.24, 5.0, decay_range_func, set(controller.GetEarthModelTargets()[0]) + 1.24, 5.0, decay_range_func, set(controller.GetDetectorModelTargets()[0]) ) primary_injection_distributions["position"] = position_distribution From 66f16d712b1d846142ecaff55c8513616915dfaf Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:53:43 -0700 Subject: [PATCH 11/14] Use resources dir --- python/LIDarkNews.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py index 826cb8845..1d9b337ca 100644 --- a/python/LIDarkNews.py +++ b/python/LIDarkNews.py @@ -11,6 +11,7 @@ import leptoninjector as LI from leptoninjector.interactions import DarkNewsCrossSection, DarkNewsDecay from leptoninjector.dataclasses import Particle +from leptoninjector import _util # DarkNews methods from DarkNews import phase_space @@ -20,6 +21,9 @@ from DarkNews.integrands import get_decay_momenta_from_vegas_samples +resources_dir = _util.resource_package_dir() + + # Class containing all upscattering and decay modes available in DarkNews class PyDarkNewsInteractionCollection: def __init__( @@ -39,12 +43,12 @@ def __init__( # Default table_dir settings if self.table_dir is None: - self.table_dir = ( - os.environ.get("LEPTONINJECTOR_SRC") - + "/resources/CrossSectionTables/DarkNewsTables/" + self.table_dir = os.path.join( + resources_dir, + "CrossSectionTables", + "DarkNewsTables", + datetime.datetime.now().strftime("%Y_%m_%d__%H:%M"), ) - ct = datetime.datetime.now().strftime("%Y_%m_%d__%H:%M/") - self.table_dir += ct # Make the table directory where will we store cross section integrators table_dir_exists = False From 7e9fb46a1abe4958f4644f496eafbeb2304a589d Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 18:54:02 -0700 Subject: [PATCH 12/14] Use resources dir --- .../Example2/DipolePortal_MiniBooNE.py | 29 ++++++++++--------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/resources/Examples/Example2/DipolePortal_MiniBooNE.py b/resources/Examples/Example2/DipolePortal_MiniBooNE.py index 7cc207996..34d0e7592 100644 --- a/resources/Examples/Example2/DipolePortal_MiniBooNE.py +++ b/resources/Examples/Example2/DipolePortal_MiniBooNE.py @@ -1,15 +1,10 @@ import os import leptoninjector as LI -import sys +from leptoninjector import _util +from leptoninjector.LIController import LIController -sys.path.insert( - 1, - "/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python", -) -from LIController import LIController - -LI_SRC = os.environ.get("LEPTONINJECTOR_SRC") +resources_dir = _util.resource_package_dir() # Define a DarkNews model model_kwargs = { @@ -37,10 +32,11 @@ primary_type = LI.dataclasses.Particle.ParticleType.NuMu # Define DarkNews Model -table_dir = ( - LI_SRC - + "/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/" - % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]) +table_dir = os.path.join( + resources_dir, + "CrossSectionTables", + "DarkNewsTables", + "Dipole_M%2.2f_mu%2.2e" % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]), ) controller.InputDarkNewsModel(primary_type, table_dir, model_kwargs) @@ -49,7 +45,12 @@ primary_physical_distributions = {} # energy distribution -flux_file = LI_SRC + "/resources/Fluxes/BNB_Flux_Tables/BNB_numu_flux.txt" +flux_file = os.path.join( + resources_dir, + "Fluxes", + "BNB_Flux_Tables", + "BNB_numu_flux.txt", +) edist = LI.distributions.TabulatedFluxDistribution(flux_file, True) edist_gen = LI.distributions.TabulatedFluxDistribution( 1.05 * model_kwargs["m4"], 10, flux_file, False @@ -71,7 +72,7 @@ model_kwargs["m4"], controller.DN_min_decay_width, 3, 541 ) position_distribution = LI.distributions.RangePositionDistribution( - 6.2, 6.2, decay_range_func, set(controller.GetEarthModelTargets()[0]) + 6.2, 6.2, decay_range_func, set(controller.GetDetectorModelTargets()[0]) ) primary_injection_distributions["position"] = position_distribution From f110c65814ac9eca5ca6d944c98cd82f4fa9dabe Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 19:04:23 -0700 Subject: [PATCH 13/14] os.path.join --- python/LIDarkNews.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py index 1d9b337ca..c05a4cb6b 100644 --- a/python/LIDarkNews.py +++ b/python/LIDarkNews.py @@ -68,10 +68,10 @@ def __init__( if param_file is not None: # ensure the param filename already exists param_filename = ntpath.basename(param_file) # should be OS-independent - assert os.path.isfile(self.table_dir + param_filename) + assert os.path.isfile(os.path.join(self.table_dir, param_filename)) # Make sure the model arguments agree with open( - self.table_dir + "model_parameters.json", + os.path.join(self.table_dir, "model_parameters.json"), ) as f: _model_args_dict = json.load(f) assert self.models.model_args_dict == _model_args_dict @@ -82,7 +82,7 @@ def __init__( command = "scp " + param_file + " " + self.table_dir os.system(command) # Dump the model arguments - with open(self.table_dir + "model_parameters.json", "w") as f: + with open(os.path.join(self.table_dir, "model_parameters.json"), "w") as f: json.dump(self.models.model_args_dict, f) # Save all unique scattering processes @@ -97,7 +97,7 @@ def __init__( self.cross_sections.append( PyDarkNewsCrossSection( ups_case, - table_dir=self.table_dir + table_subdirs, + table_dir=os.path.join(self.table_dir, table_subdirs), tolerance=tolerance, interp_tolerance=interp_tolerance, ) @@ -110,7 +110,7 @@ def __init__( table_subdirs += "%s_" % str(x) table_subdirs += "/" self.decays.append( - PyDarkNewsDecay(dec_case, table_dir=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): @@ -173,10 +173,10 @@ def __init__( # Look in table dir and check whether total/differential xsec tables exist if table_dir_exists: - total_xsec_file = self.table_dir + "total_cross_sections.npy" + 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 = 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) @@ -357,11 +357,11 @@ def FillInterpolationTables(self, total=True, diff=True): def SaveInterpolationTables(self, total=True, diff=True): if total: self._redefine_interpolation_objects(total=True) - with open(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(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 ######### @@ -594,12 +594,12 @@ def __init__(self, dec_case, table_dir=None): def SetIntegratorAndNorm(self): # Try to find the decay integrator - int_file = self.table_dir + "decay_integrator.pkl" + int_file = os.path.join(self.table_dir, "decay_integrator.pkl") if os.path.isfile(int_file): with open(int_file, "rb") as ifile: _, self.decay_integrator = pickle.load(ifile) # Try to find the normalization information - norm_file = self.table_dir + "decay_norm.json" + norm_file = os.path.join(self.table_dir, "decay_norm.json") if os.path.isfile(norm_file): with open( norm_file, @@ -717,8 +717,8 @@ def TotalDecayWidth(self, arg1): # total width calculation requires evaluating an integral if self.decay_integrator is None or self.decay_norm is None: # We need to initialize a new VEGAS integrator in DarkNews - int_file = self.table_dir + "decay_integrator.pkl" - norm_file = self.table_dir + "decay_norm.json" + int_file = os.path.join(self.table_dir, "decay_integrator.pkl") + norm_file = os.path.join(self.table_dir, "decay_norm.json") self.total_width = self.dec_case.total_width( savefile_norm=norm_file, savefile_dec=int_file ) @@ -786,8 +786,8 @@ def SampleRecordFromDarkNews(self, record, random): # We need to generate new PS samples if self.decay_integrator is None or self.decay_norm is None: # We need to initialize a new VEGAS integrator in DarkNews - int_file = self.table_dir + "decay_integrator.pkl" - norm_file = self.table_dir + "decay_norm.json" + int_file = os.path.join(self.table_dir, "decay_integrator.pkl") + norm_file = os.path.join(self.table_dir, "decay_norm.json") self.PS_samples, PS_weights_dict = self.dec_case.SamplePS( savefile_norm=norm_file, savefile_dec=int_file ) From cdd29f47cb1188b42ac23f715f712e57ff0fe109 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 22 Jan 2024 16:55:22 -0700 Subject: [PATCH 14/14] Remove printouts --- projects/injection/private/TreeWeighter.cxx | 4 ---- 1 file changed, 4 deletions(-) diff --git a/projects/injection/private/TreeWeighter.cxx b/projects/injection/private/TreeWeighter.cxx index 41072d433..30a5b5622 100644 --- a/projects/injection/private/TreeWeighter.cxx +++ b/projects/injection/private/TreeWeighter.cxx @@ -202,8 +202,6 @@ void LeptonProcessWeighter::Initialize() { } unique_gen_distributions = inj_process->GetInjectionDistributions(); unique_phys_distributions = phys_process->GetPhysicalDistributions(); - std::cout << "Num gen distributions: " << unique_gen_distributions.size() << std::endl; - std::cout << "Num phys distributions: " << unique_phys_distributions.size() << std::endl; for(std::vector>::reverse_iterator gen_it = unique_gen_distributions.rbegin(); gen_it != unique_gen_distributions.rend(); ++gen_it) { for(std::vector>::reverse_iterator phys_it = unique_phys_distributions.rbegin(); @@ -215,8 +213,6 @@ void LeptonProcessWeighter::Initialize() { } } } - std::cout << "Num gen distributions: " << unique_gen_distributions.size() << std::endl; - std::cout << "Num phys distributions: " << unique_phys_distributions.size() << std::endl; } double LeptonProcessWeighter::InteractionProbability(std::pair const & bounds, LI::dataclasses::InteractionRecord const & record) const {