diff --git a/resources/examples/example2/DipolePortal_CCM.py b/resources/examples/example2/DipolePortal_CCM.py index 76c6d8d9..6bd5fa65 100644 --- a/resources/examples/example2/DipolePortal_CCM.py +++ b/resources/examples/example2/DipolePortal_CCM.py @@ -1,13 +1,12 @@ import os import numpy as np - import siren -from siren.SIREN_Controller import SIREN_Controller +from siren import utilities # Define a DarkNews model model_kwargs = { "m4": 0.0235, - "mu_tr_mu4": 6e-7, # GeV^-1 + "mu_tr_mu4": 6e-7, # GeV^-1 "UD4": 0, "Umu4": 0, "epsilon": 0.0, @@ -18,82 +17,94 @@ } # Number of events to inject -events_to_inject = 100000 +events_to_inject = 1 -# Expeirment to run +# Experiment to run experiment = "CCM" - -# Define the controller -controller = SIREN_Controller(events_to_inject, experiment) +detector_model = utilities.load_detector(experiment) # Particle to inject primary_type = siren.dataclasses.Particle.ParticleType.NuMu -xs_path = siren.utilities.get_cross_section_model_path(f"DarkNewsTables-v{siren.utilities.darknews_version()}", must_exist=False) -# Define DarkNews Model -table_dir = os.path.join( - xs_path, - "Dipole_M%2.2e_mu%2.2e" % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]), +# Load DarkNews processes +primary_processes, secondary_processes = utilities.load_processes( + f"DarkNewsTables-v{siren.utilities.darknews_version()}", + primary_type=primary_type, + detector_model = detector_model, + **model_kwargs, ) -controller.InputDarkNewsModel(primary_type, table_dir, **model_kwargs) -# Primary distributions -primary_injection_distributions = {} -primary_physical_distributions = {} +# Mass distribution +mass_ddist = siren.distributions.PrimaryMass(0) -# energy distribution +# Primary distributions nu_energy = 0.02965 # from pi+ DAR edist = siren.distributions.Monoenergetic(nu_energy) -primary_injection_distributions["energy"] = edist -primary_physical_distributions["energy"] = edist -# fill cross section tables at this energy -controller.DN_processes.FillCrossSectionTablesAtEnergy(nu_energy) - -# Flux normalization: -# using the number quoted in 2105.14020, 4.74e9 nu/m^2/s / (6.2e14 POT/s) * 4*pi*20m^2 to get nu/POT flux_units = siren.distributions.NormalizationConstant(3.76e-2) -primary_physical_distributions["flux_units"] = flux_units -# direction distribution: cone from lower W target +# Cone direction distribution opening_angle = np.arctan(5 / 23.0) -# slightly larger than CCM lower_target_origin = siren.math.Vector3D(0, 0, -0.241) detector_origin = siren.math.Vector3D(23, 0, -0.65) lower_dir = detector_origin - lower_target_origin lower_dir.normalize() lower_inj_ddist = siren.distributions.Cone(lower_dir, opening_angle) -phys_ddist = ( - siren.distributions.IsotropicDirection() -) # truly we are isotropic -primary_injection_distributions["direction"] = lower_inj_ddist -primary_physical_distributions["direction"] = phys_ddist +phys_ddist = siren.distributions.IsotropicDirection() -# Position distribution: consider neutrinos from a point source +# Position distribution max_dist = 25 lower_pos_dist = siren.distributions.PointSourcePositionDistribution( - lower_target_origin - detector_origin, max_dist, set(controller.GetDetectorModelTargets()[0]) + lower_target_origin - detector_origin, max_dist, ) -primary_injection_distributions["position"] = lower_pos_dist - -# SetProcesses -controller.SetProcesses( - primary_type, primary_injection_distributions, primary_physical_distributions -) - -controller.Initialize() +primary_injection_distributions = [ + mass_ddist, # Mass distribution + edist, # Energy distribution + lower_inj_ddist, # Direction distribution + lower_pos_dist # Position distribution +] + +primary_physical_distributions = [ + edist, # Energy distribution + phys_ddist, # Direction distribution +] + +fiducial_volume = siren.utilities.get_fiducial_volume(experiment) +secondary_injection_distributions = {} +for secondary_type in secondary_processes.keys(): + secondary_injection_distributions[secondary_type] = [ + siren.distributions.SecondaryBoundedVertexDistribution(fiducial_volume, max_dist) + ] + +# Define stopping condition for the injector def stop(datum, i): secondary_type = datum.record.signature.secondary_types[i] return secondary_type != siren.dataclasses.Particle.ParticleType.N4 -controller.injector.SetStoppingCondition(stop) +injector = siren.injection.Injector() +injector.number_of_events = 1 +injector.detector_model = detector_model +injector.primary_type = primary_type +injector.primary_interactions = primary_processes[primary_type] +injector.primary_injection_distributions = primary_injection_distributions +injector.secondary_interactions = secondary_processes +injector.secondary_injection_distributions = secondary_injection_distributions +injector.stopping_condition = stop + -events = controller.GenerateEvents(fill_tables_at_exit=False) +# Generate events +events = [injector.generate_event() for _ in range(events_to_inject)] +# Output the events os.makedirs("output", exist_ok=True) -controller.SaveEvents( - "output/CCM_Dipole_M%2.2e_mu%2.2e_example" - % (model_kwargs["m4"], model_kwargs["mu_tr_mu4"]), - fill_tables_at_exit=False -) +weighter = siren.injection.Weighter() +weighter.injectors = [injector] +weighter.detector_model = detector_model +weighter.primary_type = primary_type +weighter.primary_interactions = primary_processes[primary_type] +weighter.secondary_interactions = secondary_processes +weighter.primary_physical_distributions = primary_physical_distributions +weighter.secondary_physical_distributions = {} + +weights = [weighter(event) for event in events]