"""
Script for simulating all 5 ARA stations, 
    triggering when 2 stations trigger on the same simulated event.

"""

import pyrex
import pyrex.custom.ara as ara

import numpy as np
import os.path
import logging

from datetime import datetime
import argparse


# Set logging
logging.basicConfig(level=logging.WARNING)


# Set command line arguments
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('energy', type=float,
                    help="neutrino energy (in GeV)")
parser.add_argument('-p', '--power_threshold', type=float, default=1,
                    help="power threshold for ant trigger (default 1)")
parser.add_argument('-n', '--number', type=int, default=10,
                    help="number of events to generate (default 10)")
parser.add_argument('--secondaries', action='store_true', 
                    help=("if present, secondary interactions will be "+
                          "considered"))
parser.add_argument('--noiseless', action='store_true', default=False,
                    help="if present, noise not added to signal")
parser.add_argument('--montecarlo', action='store_true',
                    help=("if present, suppress noise triggers based on "+
                          "MC truth"))
parser.add_argument('-o', '--output', required=True,
                    help="output file name")
args = parser.parse_args()


# Prepare output file
outfile = pyrex.File(
    args.output, 
    'w',
    write_particles=True,
    write_waveforms=False,
    write_triggers=True,
    write_antenna_triggers=False,
    write_rays=True,
    write_noise=False,
    # require_trigger=['waveforms'],
)
outfile.open()
outfile.create_analysis_metadataset("sim_parameters")

# Dictionary for metadata from the script for the output file
script_metadata = {}
script_metadata['secondaries'] = args.secondaries
script_metadata['noiselss']    = args.noiseless
script_metadata['monte_carlo'] = args.montecarlo


# Setup detectors from custom ARA objects
ara1 = ara.RegularStation(x=1e3,    y=1.732e3,        antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara2 = ara.RegularStation(x=0,      y=0, antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara3 = ara.RegularStation(x=-1e3,   y=1.732e3,        antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara4 = ara.RegularStation(x=0,      y=3.464e3,  antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara5 = ara.PhasedArrayStation(x=-1e3, y=-1.732e3, antennas_per_string = 4,
                              antenna_separation = [2, 18, 2])
stations = [ara1, ara2, ara3, ara4, ara5]

# VPol & HPol Antenna Trigger Threshold
power_threshold = -1 * abs(args.power_threshold)
script_metadata['power_threshold'] = power_threshold

# Build antennas for ARA 1-4 and set antenna names
antenna_prenames = ['B','B','T','T']
for i, station in enumerate(stations[:-1]): 
    station.build_antennas(power_threshold=power_threshold, 
                           noisy=not(args.noiseless))
    # Personalize the antenna names with station and string number
    for s, string in enumerate(station.subsets): 
        for a, antenna in enumerate(string):
            antenna.name = (f"ara{i+1}_string{s+1}_"
                            f"{antenna_prenames[a]}{antenna.name[:4]}" )
            
# Phased array (on ARA5) has a special set up, so we build it individually
stations[4].build_antennas(power_threshold=power_threshold, 
                           noisy=not(args.noiseless))
for s, string in enumerate(stations[4].subsets): 
    if s == 0:
        for antenna in string:
            antenna.name = f"ara5_PA_" + antenna.name
    else:
        for a, antenna in enumerate(string):
            antenna.name = (f"ara5_string{s}_"
                            f"{antenna_prenames[a]}{antenna.name[:4]}" )

# Compile all stations into one detector
detector = pyrex.detector.CombinedDetector(*stations)


# Individual station trigger settings
polarized_antenna_requirement = 3
require_mc_truth = args.montecarlo
beam_threshold = 2
beam_delays = None
beam_angles = None

# Set global trigger
def global_trigger(det):
    stations_triggered = 0
    for i, station in enumerate(det.subsets):
        if i+1 != 5: 
            if det.subsets[i].triggered(
                polarized_antenna_requirement=polarized_antenna_requirement,
                require_mc_truth=require_mc_truth,
            ):
                stations_triggered += 1
        elif i+1 == 5:
            if det.subsets[i].triggered(
                beam_threshold=beam_threshold,
                outrigger_antenna_requirement=polarized_antenna_requirement,
                beam_delays=beam_delays, 
                beam_angles=beam_angles, 
                require_mc_truth=require_mc_truth,
            ):
                stations_triggered += 1
    if stations_triggered > 1: 
        return True
    else: 
        return False

is_triggered = {
    "global": global_trigger,
}

# Add a triggering function for each station to is_triggered dictionary
for i, station in enumerate(detector.subsets[:-1]):
    is_triggered[f"ara{i+1}"] = (
        lambda detector, i=i:
        bool(detector.subsets[i].triggered(
             polarized_antenna_requirement=polarized_antenna_requirement,
             require_mc_truth=require_mc_truth,
        ))
    )
is_triggered[f"ara5"] = (
    lambda detector, i=4:
    bool(detector.subsets[4].triggered(
         beam_threshold=beam_threshold,
         outrigger_antenna_requirement=polarized_antenna_requirement,
         beam_delays=beam_delays, 
         beam_angles=beam_angles, 
         require_mc_truth=require_mc_truth,
    ))
)


# Set environmental and physics properties
ice = pyrex.ice_model.AntarcticIce()

earth = pyrex.earth_model.PREM()

signal_model = pyrex.askaryan.ARZAskaryanSignal

# Limit oncone angle range to reduce memory consumption
# (helpful when submitting to a cluster with 1GB of memory)
pyrex.AskaryanSignal.oncone_range = np.radians(1e-3)


# Set neutrino generator
# if 1e3<=args.energy<1e8:
#     extra = 1000
# elif 1e8<=args.energy<1e9:
#     extra = 2500
# elif 1e9<=args.energy<1e10:
#     extra = 5000
# elif 1e10<=args.energy:
#     extra = 10000
# else:
#     raise ValueError("Invalid energy: "+str(args.energy))
generator = pyrex.CylindricalGenerator(
#     dr=3000+extra, dz=3000,   # Setting for original 0930 Run
    dr=13000, dz=3000,          # Setting for finalized runs in November
    energy=args.energy,
    shadow=False,
    earth_model=earth,
    source="cosmogenic",
)

script_metadata['generator_size_r'] = generator.dr
script_metadata['generator_size_z'] = generator.dz 

# Set secondary interactions
if hasattr(generator, 'interaction_model'):
    generator.interaction_model.include_secondaries = args.secondaries

# Record additional metadata
if hasattr(generator, 'interaction_model'):
    script_metadata['interaction_model_class'] = str(
        generator.interaction_model )
    script_metadata['secondaries'] = (
        generator.interaction_model.include_secondaries )
if hasattr(generator, 'ratio'):
    script_metadata['flavor_ratio_e'] = generator.ratio[0]
    script_metadata['flavor_ratio_mu'] = generator.ratio[1]
    script_metadata['flavor_ratio_tau'] = generator.ratio[2]


# Add simulation parameters to outfile
outfile.add_analysis_metadata("sim_parameters", script_metadata)


# Prepare signal amplitude dataset
amp_dataset = outfile.create_analysis_dataset(
    "signal_amplitudes", shape=(0, len(detector), 2),
    dtype=np.float_, maxshape=(None, len(detector), 2)
)
amp_dataset.dims[0].label = "events"
amp_dataset.dims[1].label = "antennas"
amp_dataset.dims[2].label = "amplitudes"
amp_dataset.attrs['keys'] = [b"signal", b"signal+noise"]

# Prepare tunnel diode amplitude dataset
td_dataset = outfile.create_analysis_dataset(
    "tunnel_diode_data", shape=(0, len(detector), 2),
    dtype=np.float_, maxshape=(None, len(detector), 2)
)
td_dataset.dims[0].label = "events"
td_dataset.dims[1].label = "antennas"
td_dataset.dims[2].label = "min_amplitudes"
td_dataset.attrs['keys'] = [b"signal", b"signal+noise"]

# Prepare tunnel diode mean and rms dataset
td_thresholds = outfile.create_analysis_dataset(
    "tunnel_diode_thresholds", shape=(0, len(detector), 2),
    dtype=np.float_, maxshape=(None, len(detector), 2)
)
td_thresholds.dims[0].label = "events"
td_thresholds.dims[1].label = "antennas"
td_thresholds.dims[2].label = "diode_settings"
td_thresholds.attrs['keys'] = [b"mean", b"rms"]


# Build the simulation Kernel
kernel = pyrex.EventKernel(
    generator=generator,
    antennas=detector,
    ice_model=ice,
    event_writer=outfile,
    triggers=is_triggered,
    signal_model=signal_model,
    offcone_max=10,
    weight_min=(1e-5, 0),
)

# Simulate the events and save signal amplitudes from each
for i in range(args.number):
    print(i, end='..', flush=True)

    # Clear old signals from detector
    try:
        detector.clear(reset_noise=True)
    except TypeError:
        for ant in detector:
            ant.clear(reset_noise=True)

    # Process event
    event, triggered = kernel.event()

    # Store signal amplitudes
    row = amp_dataset.shape[0]
    max_waves = 2
    amp_dataset.resize(row+max_waves, axis=0)
    for col, ant in enumerate(detector):
        for j, (sig, wave) in enumerate(zip(ant.signals, ant.all_waveforms)):
            amp_dataset[row+j, col, 0] = np.max(np.abs(sig.values))
            amp_dataset[row+j, col, 1] = np.max(np.abs(wave.values))
    outfile.add_analysis_indices("signal_amplitudes", global_index=i,
                                 start_index=row, length=max_waves)
    
    # Store tunnel diode amplitudes compactly in outfile
    row = td_dataset.shape[0]
    max_waves = 2
    td_dataset.resize(row+max_waves, axis=0)
    for col, ant in enumerate(detector):
        for j, (sig, wave) in enumerate(zip(ant.signals, ant.all_waveforms)):
            td_signal = ant.tunnel_diode(sig)
            td_wave   = ant.tunnel_diode(wave)
            td_dataset[row+j, col, 0] = np.min(td_signal.values)
            td_dataset[row+j, col, 1] = np.min(td_wave.values)
    outfile.add_analysis_indices("tunnel_diode_data", global_index=i,
                                 start_index=row, length=max_waves)
    
    # Store tunnel diode thresholds compactly in outfile
    row = td_thresholds.shape[0]
    td_thresholds.resize(row+1, axis=0)
    for col, ant in enumerate(detector):
        if ant._power_mean == None:
            td_thresholds[row, col, 0] = np.min(-1e3)
            td_thresholds[row, col, 1] = np.min(-1e3)
        else:
            td_thresholds[row, col, 0] = np.min(ant._power_mean)
            td_thresholds[row, col, 1] = np.min(ant._power_std)
    outfile.add_analysis_indices("tunnel_diode_thresholds", global_index=i,
                                 start_index=row, length=1)

    # Store tunnel diode thresholds compactly in outfile
    row = 0
    for col, ant in enumerate(detector):
        if ant._power_mean == None:
            td_thresholds[row, col] = [-1000,-1000]
        else:
            td_thresholds[row, col] = [ant._power_mean, ant._power_std]
    outfile.add_analysis_indices("tunnel_diode_thresholds", global_index=i,
                                 start_index=row, length=max_waves)   
    
print()


# Gently finish simulations and close the outfile
outfile.add_file_metadata({
    "kernel_completed": datetime.now().strftime('%Y-%d-%m %H:%M:%S')
})
outfile.close()
"""
Script for simulating all 5 ARA stations, 
    triggering when 2 stations trigger on the same simulated event.

"""

import pyrex
import pyrex.custom.ara as ara

import numpy as np
import os.path
import logging

from datetime import datetime
import argparse


# Set logging
logging.basicConfig(level=logging.WARNING)


# Set command line arguments
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('energy', type=float,
                    help="neutrino energy (in GeV)")
parser.add_argument('-p', '--power_threshold', type=float, default=1,
                    help="power threshold for ant trigger (default 1)")
parser.add_argument('-n', '--number', type=int, default=10,
                    help="number of events to generate (default 10)")
parser.add_argument('--secondaries', action='store_true', 
                    help=("if present, secondary interactions will be "+
                          "considered"))
parser.add_argument('--noiseless', action='store_true', default=False,
                    help="if present, noise not added to signal")
parser.add_argument('--montecarlo', action='store_true',
                    help=("if present, suppress noise triggers based on "+
                          "MC truth"))
parser.add_argument('-o', '--output', required=True,
                    help="output file name")
args = parser.parse_args()


# Prepare output file
outfile = pyrex.File(
    args.output, 
    'w',
    write_particles=True,
    write_waveforms=False,
    write_triggers=True,
    write_antenna_triggers=False,
    write_rays=True,
    write_noise=False,
    # require_trigger=['waveforms'],
)
outfile.open()
outfile.create_analysis_metadataset("sim_parameters")

# Dictionary for metadata from the script for the output file
script_metadata = {}
script_metadata['secondaries'] = args.secondaries
script_metadata['noiselss']    = args.noiseless
script_metadata['monte_carlo'] = args.montecarlo


# Setup detectors from custom ARA objects
ara1 = ara.RegularStation(x=1e3,    y=1.732e3,        antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara2 = ara.RegularStation(x=0,      y=0, antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara3 = ara.RegularStation(x=-1e3,   y=1.732e3,        antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara4 = ara.RegularStation(x=0,      y=3.464e3,  antennas_per_string=4,
                         antenna_separation=[2, 18, 2])
ara5 = ara.PhasedArrayStation(x=-1e3, y=-1.732e3, antennas_per_string = 4,
                              antenna_separation = [2, 18, 2])
stations = [ara1, ara2, ara3, ara4, ara5]

# VPol & HPol Antenna Trigger Threshold
power_threshold = -1 * abs(args.power_threshold)
script_metadata['power_threshold'] = power_threshold

# Build antennas for ARA 1-4 and set antenna names
antenna_prenames = ['B','B','T','T']
for i, station in enumerate(stations[:-1]): 
    station.build_antennas(power_threshold=power_threshold, 
                           noisy=not(args.noiseless))
    # Personalize the antenna names with station and string number
    for s, string in enumerate(station.subsets): 
        for a, antenna in enumerate(string):
            antenna.name = (f"ara{i+1}_string{s+1}_"
                            f"{antenna_prenames[a]}{antenna.name[:4]}" )
            
# Phased array (on ARA5) has a special set up, so we build it individually
stations[4].build_antennas(power_threshold=power_threshold, 
                           noisy=not(args.noiseless))
for s, string in enumerate(stations[4].subsets): 
    if s == 0:
        for antenna in string:
            antenna.name = f"ara5_PA_" + antenna.name
    else:
        for a, antenna in enumerate(string):
            antenna.name = (f"ara5_string{s}_"
                            f"{antenna_prenames[a]}{antenna.name[:4]}" )

# Compile all stations into one detector
detector = pyrex.detector.CombinedDetector(*stations)


# Individual station trigger settings
polarized_antenna_requirement = 3
require_mc_truth = args.montecarlo
beam_threshold = 2
beam_delays = None
beam_angles = None

# Set global trigger
def global_trigger(det):
    stations_triggered = 0
    for i, station in enumerate(det.subsets):
        if i+1 != 5: 
            if det.subsets[i].triggered(
                polarized_antenna_requirement=polarized_antenna_requirement,
                require_mc_truth=require_mc_truth,
            ):
                stations_triggered += 1
        elif i+1 == 5:
            if det.subsets[i].triggered(
                beam_threshold=beam_threshold,
                outrigger_antenna_requirement=polarized_antenna_requirement,
                beam_delays=beam_delays, 
                beam_angles=beam_angles, 
                require_mc_truth=require_mc_truth,
            ):
                stations_triggered += 1
    if stations_triggered > 1: 
        return True
    else: 
        return False

is_triggered = {
    "global": global_trigger,
}

# Add a triggering function for each station to is_triggered dictionary
for i, station in enumerate(detector.subsets[:-1]):
    is_triggered[f"ara{i+1}"] = (
        lambda detector, i=i:
        bool(detector.subsets[i].triggered(
             polarized_antenna_requirement=polarized_antenna_requirement,
             require_mc_truth=require_mc_truth,
        ))
    )
is_triggered[f"ara5"] = (
    lambda detector, i=4:
    bool(detector.subsets[4].triggered(
         beam_threshold=beam_threshold,
         outrigger_antenna_requirement=polarized_antenna_requirement,
         beam_delays=beam_delays, 
         beam_angles=beam_angles, 
         require_mc_truth=require_mc_truth,
    ))
)


# Set environmental and physics properties
ice = pyrex.ice_model.AntarcticIce()

earth = pyrex.earth_model.PREM()

signal_model = pyrex.askaryan.ARZAskaryanSignal

# Limit oncone angle range to reduce memory consumption
# (helpful when submitting to a cluster with 1GB of memory)
pyrex.AskaryanSignal.oncone_range = np.radians(1e-3)


# Set neutrino generator
# if 1e3<=args.energy<1e8:
#     extra = 1000
# elif 1e8<=args.energy<1e9:
#     extra = 2500
# elif 1e9<=args.energy<1e10:
#     extra = 5000
# elif 1e10<=args.energy:
#     extra = 10000
# else:
#     raise ValueError("Invalid energy: "+str(args.energy))
generator = pyrex.CylindricalGenerator(
#     dr=3000+extra, dz=3000,   # Setting for original 0930 Run
    dr=13000, dz=3000,          # Setting for finalized runs in November
    energy=args.energy,
    shadow=False,
    earth_model=earth,
    source="cosmogenic",
)

script_metadata['generator_size_r'] = generator.dr
script_metadata['generator_size_z'] = generator.dz 

# Set secondary interactions
if hasattr(generator, 'interaction_model'):
    generator.interaction_model.include_secondaries = args.secondaries

# Record additional metadata
if hasattr(generator, 'interaction_model'):
    script_metadata['interaction_model_class'] = str(
        generator.interaction_model )
    script_metadata['secondaries'] = (
        generator.interaction_model.include_secondaries )
if hasattr(generator, 'ratio'):
    script_metadata['flavor_ratio_e'] = generator.ratio[0]
    script_metadata['flavor_ratio_mu'] = generator.ratio[1]
    script_metadata['flavor_ratio_tau'] = generator.ratio[2]


# Add simulation parameters to outfile
outfile.add_analysis_metadata("sim_parameters", script_metadata)


# Prepare signal amplitude dataset
amp_dataset = outfile.create_analysis_dataset(
    "signal_amplitudes", shape=(0, len(detector), 2),
    dtype=np.float_, maxshape=(None, len(detector), 2)
)
amp_dataset.dims[0].label = "events"
amp_dataset.dims[1].label = "antennas"
amp_dataset.dims[2].label = "amplitudes"
amp_dataset.attrs['keys'] = [b"signal", b"signal+noise"]

# Prepare tunnel diode amplitude dataset
td_dataset = outfile.create_analysis_dataset(
    "tunnel_diode_data", shape=(0, len(detector), 2),
    dtype=np.float_, maxshape=(None, len(detector), 2)
)
td_dataset.dims[0].label = "events"
td_dataset.dims[1].label = "antennas"
td_dataset.dims[2].label = "min_amplitudes"
td_dataset.attrs['keys'] = [b"signal", b"signal+noise"]

# Prepare tunnel diode mean and rms dataset
td_thresholds = outfile.create_analysis_dataset(
    "tunnel_diode_thresholds", shape=(0, len(detector), 2),
    dtype=np.float_, maxshape=(None, len(detector), 2)
)
td_thresholds.dims[0].label = "events"
td_thresholds.dims[1].label = "antennas"
td_thresholds.dims[2].label = "diode_settings"
td_thresholds.attrs['keys'] = [b"mean", b"rms"]


# Build the simulation Kernel
kernel = pyrex.EventKernel(
    generator=generator,
    antennas=detector,
    ice_model=ice,
    event_writer=outfile,
    triggers=is_triggered,
    signal_model=signal_model,
    offcone_max=10,
    weight_min=(1e-5, 0),
)

# Simulate the events and save signal amplitudes from each
for i in range(args.number):
    print(i, end='..', flush=True)

    # Clear old signals from detector
    try:
        detector.clear(reset_noise=True)
    except TypeError:
        for ant in detector:
            ant.clear(reset_noise=True)

    # Process event
    event, triggered = kernel.event()

    # Store signal amplitudes
    row = amp_dataset.shape[0]
    max_waves = 2
    amp_dataset.resize(row+max_waves, axis=0)
    for col, ant in enumerate(detector):
        for j, (sig, wave) in enumerate(zip(ant.signals, ant.all_waveforms)):
            amp_dataset[row+j, col, 0] = np.max(np.abs(sig.values))
            amp_dataset[row+j, col, 1] = np.max(np.abs(wave.values))
    outfile.add_analysis_indices("signal_amplitudes", global_index=i,
                                 start_index=row, length=max_waves)
    
    # Store tunnel diode amplitudes compactly in outfile
    row = td_dataset.shape[0]
    max_waves = 2
    td_dataset.resize(row+max_waves, axis=0)
    for col, ant in enumerate(detector):
        for j, (sig, wave) in enumerate(zip(ant.signals, ant.all_waveforms)):
            td_signal = ant.tunnel_diode(sig)
            td_wave   = ant.tunnel_diode(wave)
            td_dataset[row+j, col, 0] = np.min(td_signal.values)
            td_dataset[row+j, col, 1] = np.min(td_wave.values)
    outfile.add_analysis_indices("tunnel_diode_data", global_index=i,
                                 start_index=row, length=max_waves)
    
    # Store tunnel diode thresholds compactly in outfile
    row = td_thresholds.shape[0]
    td_thresholds.resize(row+1, axis=0)
    for col, ant in enumerate(detector):
        if ant._power_mean == None:
            td_thresholds[row, col, 0] = np.min(-1e3)
            td_thresholds[row, col, 1] = np.min(-1e3)
        else:
            td_thresholds[row, col, 0] = np.min(ant._power_mean)
            td_thresholds[row, col, 1] = np.min(ant._power_std)
    outfile.add_analysis_indices("tunnel_diode_thresholds", global_index=i,
                                 start_index=row, length=1)

    # Store tunnel diode thresholds compactly in outfile
    row = 0
    for col, ant in enumerate(detector):
        if ant._power_mean == None:
            td_thresholds[row, col] = [-1000,-1000]
        else:
            td_thresholds[row, col] = [ant._power_mean, ant._power_std]
    outfile.add_analysis_indices("tunnel_diode_thresholds", global_index=i,
                                 start_index=row, length=max_waves)   
    
print()


# Gently finish simulations and close the outfile
outfile.add_file_metadata({
    "kernel_completed": datetime.now().strftime('%Y-%d-%m %H:%M:%S')
})
outfile.close()
Script for simulating all 5 ARA stations,      triggering when 2 stations trigger on the same simulated event. # Set logging # Set command line arguments # Prepare output file # require_trigger=['waveforms'], # Dictionary for metadata from the script for the output file # Setup detectors from custom ARA objects # VPol & HPol Antenna Trigger Threshold # Build antennas for ARA 1-4 and set antenna names # Personalize the antenna names with station and string number # Phased array (on ARA5) has a special set up, so we build it individually # Compile all stations into one detector # Individual station trigger settings # Set global trigger # Add a triggering function for each station to is_triggered dictionary # Set environmental and physics properties # Limit oncone angle range to reduce memory consumption # (helpful when submitting to a cluster with 1GB of memory) # Set neutrino generator # if 1e3<=args.energy<1e8: #     extra = 1000 # elif 1e8<=args.energy<1e9: #     extra = 2500 # elif 1e9<=args.energy<1e10: #     extra = 5000 # elif 1e10<=args.energy: #     extra = 10000 # else: #     raise ValueError("Invalid energy: "+str(args.energy)) #     dr=3000+extra, dz=3000,   # Setting for original 0930 Run # Setting for finalized runs in November # Set secondary interactions # Record additional metadata # Add simulation parameters to outfile # Prepare signal amplitude dataset # Prepare tunnel diode amplitude dataset # Prepare tunnel diode mean and rms dataset # Build the simulation Kernel # Simulate the events and save signal amplitudes from each # Clear old signals from detector # Process event # Store signal amplitudes # Store tunnel diode amplitudes compactly in outfile # Store tunnel diode thresholds compactly in outfile # Store tunnel diode thresholds compactly in outfile # Gently finish simulations and close the outfile