From 1346d373d1a03e631657dbfe4948b04dd6f02834 Mon Sep 17 00:00:00 2001 From: ibsafa Date: Tue, 16 Apr 2024 16:02:31 -0400 Subject: [PATCH] initial fixes --- taurunner/main.py | 69 +++++----- taurunner/track/utils/spline_column_depth.py | 2 +- taurunner/utils/__init__.py | 3 +- taurunner/utils/make_propagator.py | 125 ------------------- 4 files changed, 36 insertions(+), 163 deletions(-) delete mode 100644 taurunner/utils/make_propagator.py diff --git a/taurunner/main.py b/taurunner/main.py index d0bb355..d77558c 100755 --- a/taurunner/main.py +++ b/taurunner/main.py @@ -21,20 +21,20 @@ def initialize_parser(): # pragma: no cover parser.add_argument( '-s', dest='seed', - type=int, + type=int, default=None, help='just an integer seed to help with output file names' ) parser.add_argument( - '-n', - dest='nevents', - type=int, + '-n', + dest='nevents', + type=int, default=0, help='how many events do you want?' ) parser.add_argument( - '--flavor', - dest='flavor', + '--flavor', + dest='flavor', type=int, default=16, help='neutrino flavor (default is nutau): 12 for nue 14 for numu 16 for nutau' ) @@ -47,7 +47,7 @@ def initialize_parser(): # pragma: no cover ) # Energy arguments parser.add_argument( - '-e', + '-e', dest='energy', default='', help='Energy in GeV if numerical value greater than 0 passed.\n\ @@ -69,7 +69,7 @@ def initialize_parser(): # pragma: no cover # Angular arguments parser.add_argument( - '-t', + '-t', dest='theta', default='', help='nadir angle in degrees if numerical value(0 is through the core).\n\ @@ -78,23 +78,23 @@ def initialize_parser(): # pragma: no cover parser.add_argument( '--th_max', dest='th_max', - type=float, + type=float, default=90, help='If doing a theta range, maximum theta value. Default 90, i.e. skimming' ) parser.add_argument( - '--th_min', + '--th_min', type=float, default=0, help='If doing a theta range, minimum theta value. Default 0, i.e. through the core' ) - + # Saving arguments parser.add_argument( - '--save', - dest='save', - type=str, - default='', + '--save', + dest='save', + type=str, + default='', help="If saving output, provide a path here, if not specified, output will be printed" ) @@ -132,10 +132,9 @@ def initialize_parser(): # pragma: no cover default=0, help="Depth of the detector in km." ) - # Options parser.add_argument('--no_losses', - dest='no_losses', + dest='no_losses', default=False, action='store_true', help="Raise this flag if you want to turn off tau losses. In this case, taus will decay at rest." @@ -145,12 +144,12 @@ def initialize_parser(): # pragma: no cover action='store_true', help="Raise this flag to turn off secondaries" ) - parser.add_argument('-d', - dest='debug', - default=False, - action='store_true', + parser.add_argument('-d', + dest='debug', + default=False, + action='store_true', help='Do you want to print out debug statments? If so, raise this flag' - ) + ) parser.add_argument('--e_cut', dest='e_cut', default=0.0, @@ -162,10 +161,10 @@ def initialize_parser(): # pragma: no cover return args def run_MC( - eini: np.ndarray, + eini: np.ndarray, thetas: np.ndarray, - body: Body, - xs: CrossSections, + body: Body, + xs: CrossSections, seed: int = 0, no_secondaries: bool = False, flavor: int = 16, @@ -203,7 +202,7 @@ def run_MC( prev_th = thetas[0] prev_track = make_track(prev_th) secondary_basket = [] - idxx = [] + idxx = [] my_track = None prv_theta = np.nan # Run the algorithm @@ -218,16 +217,16 @@ def run_MC( if (cur_theta!=prv_theta and track_type=='chord') or my_track is None: # We need to make a new track my_track = getattr(track, track_type)(theta=cur_theta, depth=depth) particle = Particle( - particleIDs[i], - cur_e, + particleIDs[i], + cur_e, 0.0 , xs, - not no_secondaries, + not no_secondaries, no_losses ) - + out = Propagate(particle, my_track, body, clp, condition=condition) - + if (out.survived==False): #this muon/electron was absorbed. we record it in the output with outgoing energy 0 output.append((cur_e, 0., cur_theta, out.nCC, out.nNC, out.ID, i, out.position)) @@ -243,7 +242,7 @@ def run_MC( del out del particle idxx = np.asarray(idxx).astype(np.int32) - if not no_secondaries: + if not no_secondaries: #make muon propagator secondary_basket = np.concatenate(secondary_basket) for sec, i in zip(secondary_basket, idxx): @@ -318,7 +317,7 @@ def run_MC( savedir = '/'.join(TR_specs['save'].split('/')[:-1]) if not os.path.isdir(savedir): raise ValueError('Savedir %s does not exist' % TR_specs['save']) - + # Set up a random state rand = np.random.RandomState(TR_specs['seed']) @@ -330,7 +329,7 @@ def run_MC( # Make an array of injected energies from taurunner.utils.make_initial_e import make_initial_e eini = make_initial_e(TR_specs['nevents'], TR_specs['energy'], - e_min=TR_specs['e_min'], e_max=TR_specs['e_max'], rand=rand) + e_min=TR_specs['e_min'], e_max=TR_specs['e_max']) # Make an array of injected incident angles from taurunner.utils.make_initial_thetas import make_initial_thetas @@ -341,7 +340,7 @@ def run_MC( thetas = make_initial_thetas(TR_specs['nevents'], theta, rand=rand, track_type=TR_specs['track']) sorter = np.argsort(thetas) thetas = thetas[sorter] - + xs = CrossSections(getattr(XSModel, TR_specs['xs_model'].upper())) prop = make_propagator(TR_specs['flavor'] - np.sign(TR_specs["flavor"]), body, TR_specs['xs_model']) diff --git a/taurunner/track/utils/spline_column_depth.py b/taurunner/track/utils/spline_column_depth.py index 34a8a77..ee9a397 100644 --- a/taurunner/track/utils/spline_column_depth.py +++ b/taurunner/track/utils/spline_column_depth.py @@ -11,7 +11,7 @@ import os import pickle as pkl -from taurunner.utils import FileLock +from taurunner.utils.file_lock import FileLock import taurunner as tr diff --git a/taurunner/utils/__init__.py b/taurunner/utils/__init__.py index 9d80405..3973295 100644 --- a/taurunner/utils/__init__.py +++ b/taurunner/utils/__init__.py @@ -5,6 +5,5 @@ from .construct_body import construct_body from .setup_outdir import setup_outdir from .sample_powerlaw import sample_powerlaw -from .make_propagator import make_propagator from .make_initial_thetas import make_initial_thetas -from .make_initial_e import make_initial_e \ No newline at end of file +from .make_initial_e import make_initial_e diff --git a/taurunner/utils/make_propagator.py b/taurunner/utils/make_propagator.py deleted file mode 100644 index 7bc647d..0000000 --- a/taurunner/utils/make_propagator.py +++ /dev/null @@ -1,125 +0,0 @@ -import os -import numpy as np -import proposal as pp -from scipy.optimize import ridder -import sys -# In Python 3.6 and before importlib.resources is importlib_resources -if sys.version_info.major==3 and sys.version_info.minor<=6: - from importlib_resources import path -else: - from importlib.resources import path - -from scipy.integrate import quad -from taurunner.utils import units - -ID_2_name = {16: 'TauMinusDef', -16: 'TauPlusDef', - 14: 'MuMinusDef', -14: 'MuPlusDef', - 12: 'EMinusDef', -12: 'EPlusDef',} - -def segment_body(body, granularity=0.5): - descs = [] - for xi, xf in zip(body.layer_boundaries[:-1], body.layer_boundaries[1:]): - f_density = body.get_density(xf, right=True) - gran = np.abs(f_density - body.get_density(xi, right=False)) / body.get_density(xi, right=False) - if gran < granularity: # the percent difference within a layer is small - descs.append((xi, xf, body.get_average_density(0.5*(xi+xf)))) - else: - start = xi - s_density = body.get_density(start, right=True) - while np.abs(f_density-s_density)/s_density-granularity>0: - func = lambda x: np.abs(body.get_density(x, right=True)-s_density)/s_density-granularity - end = ridder(func, start, xf) - wrap = lambda x: body.get_density(x, right=True) - I = quad(wrap, start, end, full_output=1) - avg_density = I[0] / (end-start) - descs.append((start, end, avg_density)) - start = end - s_density = body.get_density(start) - return descs - -def make_propagator( - particle: Particle, - simulation_specs: dict, - path_dict: dict -) -> pp.Propagator: - """Make a PROPOSAL propagator - - params - ______ - particle: Prometheus particle for which we want a PROPOSAL propagator - simulation_specs: Dictionary specifying the configuration settings - path_dict: Dictionary specifying any required path variables - - returns - _______ - prop: PROPOSAL propagator for input Particle - """ - - pp.InterpolationSettings.tables_path = path_dict["tables path"] - pdef = make_particle_definition(particle) - utilities = make_propagation_utilities( - pdef, - path_dict["earth model location"], - simulation_specs - ) - geometries = make_geometries(path_dict["earth model location"]) - density_distrs = make_density_distributions(path_dict["earth model location"]) - prop = pp.Propagator(pdef, list(zip(geometries, utilities, density_distrs))) - - return prop - - -def make_particle_definition(particle: Particle) -> pp.particle.ParticleDef: - ''' - Builds a proposal particle definition - - Parameters - ---------- - particle: Prometheus particle you want a ParticleDef for - - Returns - ------- - pdef: PROPOSAL particle definition object corresponing - to input particle - ''' - if str(particle) not in 'MuMinus MuPlus EMinus EPlus TauMinus TauPlus'.split(): - raise ValueError(f"Particle string {particle} not recognized") - pdef = getattr(pp.particle, f'{particle}Def')() - return pdef - - -def make_sector(density, start, end, xs_model): - #Define a sector - sec_def = pp.SectorDefinition() - components = [ - pp.component.Hydrogen(2), - pp.component.Oxygen() - ] - sec_def.medium = pp.medium.Medium( - f'Ice_{density}', - 1, - 75.0, - -3.5017, - 0.09116, - 3.4773, - 0.2400, - 2.8004, - 0, - density, - components - ) - sec_def.geometry = pp.geometry.Sphere(pp.Vector3D(), end, start) - sec_def.particle_location = pp.ParticleLocation.inside_detector - sec_def.scattering_model = pp.scattering.ScatteringModel.Moliere - sec_def.crosssection_defs.brems_def.lpm_effect = True - sec_def.crosssection_defs.epair_def.lpm_effect = True - sec_def.cut_settings.ecut = -1.0 - sec_def.cut_settings.vcut = 1e-3 - sec_def.do_continuous_randomization = True - - if(xs_model=='dipole'): - sec_def.crosssection_defs.photo_def.parametrization = pp.parametrization.photonuclear.PhotoParametrization.BlockDurandHa - else: - sec_def.crosssection_defs.photo_def.parametrization = pp.parametrization.photonuclear.PhotoParametrization.AbramowiczLevinLevyMaor97 - - return sec_def