Welcome to Flarestack’s documentation!¶
Base PDFs¶
Time PDFs¶

class
flarestack.core.time_pdf.
Box
(t_pdf_dict, season)[source]¶ The simplest timedependent case for a Time PDF. Used for a source that is uniformly emitting for a fixed period of time. Requires arguments of PreWindow and Post_window, and gives a box from PreWindow days before the reference time to PostWindow days after the reference time.

effective_injection_time
(source)[source]¶ Calculates the effective injection time for the given PDF. The livetime is measured in days, but here is converted to seconds.
Parameters: source – Source to be considered Returns: Effective Livetime in seconds

flare_time_mask
(source)[source]¶ In this case, the interesting period for Flare Searches is the period of overlap of the flare and the box. Thus, for a given season, return the source and data
Returns: Start time (MJD) and End Time (MJD) for flare search period

raw_injection_time
(source)[source]¶ Calculates the ‘raw injection time’ which is the injection time assuming a detector with 100% uptime. Useful for calculating source emission times for sourceframe energy estimation.
Parameters: source – Source to be considered Returns: Time in seconds for 100% uptime

sig_t0
(source)[source]¶ Calculates the starting time for the window, equal to the source reference time in MJD minus the length of the prereferencetime window (in days).
Parameters: source – Source to be considered Returns: Time of Window Start

sig_t1
(source)[source]¶ Calculates the starting time for the window, equal to the source reference time in MJD plus the length of the postreferencetime window (in days).
Parameters: source – Source to be considered Returns: Time of Window End

signal_f
(t, source)[source]¶ In this case, the signal PDF is a uniform PDF for a fixed duration of time. It is normalised with the length of the box in LIVETIME rather than days, to give an integral of 1.
Parameters:  t – Time
 source – Source to be considered
Returns: Value of normalised box function at t

signal_integral
(t, source)[source]¶ In this case, the signal PDF is a uniform PDF for a fixed duration of time. Thus, the integral is simply a linear function increasing between t0 (box start) and t1 (box end). After t1, the integral is equal to 1, while it is equal to 0 for t < t0.
Parameters:  t – Time
 source – Source to be considered
Returns: Value of normalised box function at t


class
flarestack.core.time_pdf.
FixedEndBox
(t_pdf_dict, season)[source]¶ The simplest timedependent case for a Time PDF. Used for a source that is uniformly emitting for a fixed period of time. In this case, the start and end time for the box is unique for each source. The sources must have a field “Start Time (MJD)” and another “End Time (MJD)”, specifying the period of the Time PDF.

class
flarestack.core.time_pdf.
FixedRefBox
(t_pdf_dict, season)[source]¶ The simplest timedependent case for a Time PDF. Used for a source that is uniformly emitting for a fixed period of time. In this case, the start and end time for the box is unique for each source. The sources must have a field “Start Time (MJD)” and another “End Time (MJD)”, specifying the period of the Time PDF.

class
flarestack.core.time_pdf.
Steady
(t_pdf_dict, season)[source]¶ The timeindependent case for a Time PDF. Requires no additional arguments in the dictionary for __init__. Used for a steady source that is continuously emitting.

effective_injection_time
(source)[source]¶ Calculates the effective injection time for the given PDF. The livetime is measured in days, but here is converted to seconds.
Parameters: source – Source to be considered Returns: Effective Livetime in seconds

flare_time_mask
(source)[source]¶ In this case, the interesting period for Flare Searches is the entire season. Thus returns the start and end times for the season.
Returns: Start time (MJD) and End Time (MJD) for flare search period

raw_injection_time
(source)[source]¶ Calculates the ‘raw injection time’ which is the injection time assuming a detector with 100% uptime. Useful for calculating source emission times for sourceframe energy estimation.
Parameters: source – Source to be considered Returns: Time in seconds for 100% uptime

sig_t0
(source)[source]¶ Calculates the starting time for the window, equal to the source reference time in MJD minus the length of the prereferencetime window (in days).
Parameters: source – Source to be considered Returns: Time of Window Start

sig_t1
(source)[source]¶ Calculates the starting time for the window, equal to the source reference time in MJD plus the length of the postreferencetime window (in days).
Parameters: source – Source to be considered Returns: Time of Window End

signal_f
(t, source)[source]¶ In the case of a steady source, the signal PDF is a uniform PDF in time. It is thus simply equal to the season_f, normalised with the length of the season to give an integral of 1. It is thus equal to the background PDF.
Parameters:  t – Time
 source – Source to be considered
Returns: Value of normalised box function at t

signal_integral
(t, source)[source]¶ In the case of a steady source, the signal PDF is a uniform PDF in time. Thus, the integral is simply a linear function increasing between t0 (box start) and t1 (box end). After t1, the integral is equal to 1, while it is equal to 0 for t < t0.
Parameters:  t – Time
 source – Source to be considered
Returns: Value of normalised box function at t


class
flarestack.core.time_pdf.
TimePDF
(t_pdf_dict, season)[source]¶ 

background_f
(t, source)[source]¶ In all cases, we assume that the background is uniform in time. Thus, the background PDF is just a normalised version of the season_f box function.
Parameters:  t – Time
 source – Source to be considered
Returns: Value of normalised box function at t

inverse_interpolate
(source)[source]¶ Calculates the values for the integral of the signal PDF within the season. Then rescales these values, such that the start of the season yields 0, and then end of the season yields 1. Creates a function to interpolate between these values. Then, for a number between 0 and 1, the interpolated function will return the MJD time at which that fraction of the cumulative distribution was reached.
Parameters: source – Source to be considered Returns: Interpolated function

product_integral
(t, source)[source]¶ Calculates the product of the given signal PDF with the season box function. Thus gives 0 everywhere outside the season, and otherwise the value of the normalised integral. The season function is offset by 1e9, to ensure that f(t1) is equal to 1. (i.e the function is equal to 1 at the end of the box).
Parameters:  t – Time
 source – Source to be considered
Returns: Product of signal integral and season

simulate_times
(source, n_s)[source]¶ Randomly draws times for n_s events for a given source, all lying within the current season. The values are based on an interpolation of the integrated time PDF.
Parameters:  source – Source being considered
 n_s – Number of event times to be simulated
Returns: Array of times in MJD for a given source

subclasses
= {'Box': <class 'flarestack.core.time_pdf.Box'>, 'FixedEndBox': <class 'flarestack.core.time_pdf.FixedEndBox'>, 'FixedRefBox': <class 'flarestack.core.time_pdf.FixedRefBox'>, 'Steady': <class 'flarestack.core.time_pdf.Steady'>}¶

Energy PDFs¶
This script contains the EnergyPDF classes, that are used for weighting events based on a given energy PDF.

class
flarestack.core.energy_pdf.
EnergyPDF
(e_pdf_dict)[source]¶ 

fluence_integral
()[source]¶ Performs an integral for fluence over a given energy range. This is gives the total energy per unit area per second that is radiated.

integrate_over_E
(f, lower=None, upper=None)[source]¶ Uses Newton’s method to integrate function f over the energy range. By default, uses 100GeV to 10PeV, unless otherwise specified. Uses 1000 logarithmicallyspaced bins to calculate integral.
Parameters: f – Function to be integrated Returns: Integral of function

classmethod
register_subclass
(energy_pdf_name)[source]¶ Adds a new subclass of EnergyPDF, with class name equal to “energy_pdf_name”.

subclasses
= {'PowerLaw': <class 'flarestack.core.energy_pdf.PowerLaw'>, 'Spline': <class 'flarestack.core.energy_pdf.Spline'>}¶


class
flarestack.core.energy_pdf.
PowerLaw
(e_pdf_dict=None)[source]¶ A Power Law energy PDF. Takes an argument of gamma in the dictionary for the init function, where gamma is the spectral index of the Power Law.

__init__
(e_pdf_dict=None)[source]¶ Creates a PowerLaw object, which is an energy PDF based on a power law. The power law is generated from e_pdf_dict, which can specify a spectral index (Gamma), as well as an optional minimum energy (E Min) and a maximum energy (E Max)
Parameters: e_pdf_dict – Dictionary containing parameters

fluence_integral
()[source]¶ Performs an integral for fluence over a given energy range. This is gives the total energy per unit area per second that is radiated.

weight_mc
(mc, gamma=None)[source]¶ Returns an array containing the weights for each MC event, given that the spectral index gamma has been chosen. Weights each event as (E/GeV)^gamma, and multiplies this by the preexisting MC oneweight value, to give the overall oneweight.
Parameters:  mc – Monte Carlo
 gamma – Spectral Index (default is value in e_pdf_dict)
Returns: Weights Array


class
flarestack.core.energy_pdf.
Spline
(e_pdf_dict={})[source]¶ A Power Law energy PDF. Takes an argument of gamma in the dictionary for the init function, where gamma is the spectral index of the Power Law.

__init__
(e_pdf_dict={})[source]¶ Creates a PowerLaw object, which is an energy PDF based on a power law. The power law is generated from e_pdf_dict, which can specify a spectral index (Gamma), as well as an optional minimum energy (E Min) and a maximum energy (E Max)
Parameters: e_pdf_dict – Dictionary containing parameters

weight_mc
(mc)[source]¶ Returns an array containing the weights for each MC event, given that the spectral index gamma has been chosen. Weights each event using the energy spline, and multiplies this by the preexisting MC oneweight value, to give the overall oneweight.
Parameters: mc – Monte Carlo Returns: Weights Array

Spatial PDFs¶

class
flarestack.core.spatial_pdf.
CircularGaussian
(spatial_pdf_dict)[source]¶ 
static
signal_spatial
(source, cut_data)[source]¶ Calculates the angular distance between the source and the coincident dataset. Uses a Gaussian PDF function, centered on the source. Returns the value of the Gaussian at the given distances.
Parameters:  source – Single Source
 cut_data – Subset of Dataset with coincident events
Returns: Array of Spatial PDF values

static

class
flarestack.core.spatial_pdf.
SpatialPDF
(spatial_pdf_dict)[source]¶ 

classmethod
register_subclass
(spatial_pdf_name)[source]¶ Adds a new subclass of SpatialPDF, with class name equal to “spatial_pdf_name”.

static
rotate
(ra1, dec1, ra2, dec2, ra3, dec3)[source]¶ Rotate ra1 and dec1 in a way that ra2 and dec2 will exactly map onto ra3 and dec3, respectively. All angles are treated as radians. Essentially rotates the events, so that they behave as if they were originally incident on the source.
Parameters:  ra1 – Event Right Ascension
 dec1 – Event Declination
 ra2 – True Event Right Ascension
 dec2 – True Event Declination
 ra3 – Source Right Ascension
 dec3 – Source Declination
Returns: Returns new Right Ascensions and Declinations

rotate_to_position
(ev, ra, dec)[source]¶ Modifies the events by reassigning the Right Ascension and Declination of the events. Rotates the events, so that they are distributed as if they originated from the source. Removes the additional Monte Carlo information from sampled events, so that they appear like regular data.
 The fields removed are:
 True Right Ascension, True Declination, True Energy, OneWeight
Parameters:  ev – Events
 ra – Source Right Ascension (radians)
 dec – Source Declination (radians)
Returns: Events (modified)

subclasses
= {'circular_gaussian': <class 'flarestack.core.spatial_pdf.CircularGaussian'>}¶

classmethod
Composite PDF Objects¶
Injector¶

class
flarestack.core.injector.
BaseInjector
(season, sources, **kwargs)[source]¶ Base Injector Class

__init__
(season, sources, **kwargs)[source]¶ Initialize self. See help(type(self)) for accurate signature.

create_dataset
(scale, pull_corrector)[source]¶ Create a dataset based on scrambled data for background, and Monte Carlo simulation for signal. Returns the composite dataset. The source flux can be scaled by the scale parameter.
Parameters: scale – Ratio of Injected Flux to source flux Returns: Simulated dataset

classmethod
register_subclass
(inj_name)[source]¶ Adds a new subclass of EnergyPDF, with class name equal to “energy_pdf_name”.

subclasses
= {}¶


class
flarestack.core.injector.
EffectiveAreaInjector
(season, sources, **kwargs)[source]¶ Class for injecting signal events by relying on effective areas rather than preexisting Monte Carlo simulation. This Injector should be used for analysing public data, as no MC is provided.

class
flarestack.core.injector.
LowMemoryInjector
(season, sources, **kwargs)[source]¶ For large numbers of sources O(~100), saving MC masks becomes increasingly burdensome. As a solution, the LowMemoryInjector should be used instead. It will be somewhat slower, but will have much more reasonable memory consumption.

class
flarestack.core.injector.
MCInjector
(season, sources, **kwargs)[source]¶ Core Injector Class, returns a dataset on which calculations can be performed. This base class is tailored for injection of MC into mock background. This can be either MC background, or scrambled real data.

__init__
(season, sources, **kwargs)[source]¶ Initialize self. See help(type(self)) for accurate signature.

calculate_fluence
(source, scale, source_mc, band_mask, omega)[source]¶ Function to calculate the fluence for a given source, and multiply the oneweights by this. After this step, the oneweight sum is equal to the expected neutrino number.
Parameters:  source – Source to be calculated
 scale – Flux scale
 source_mc – MC that is close to source
 band_mask – Closeness mask for MC
 omega – Solid angle covered by MC mask
Returns: Modified source MC

calculate_single_source
(source, scale)[source]¶ Calculate the weighted MC for a single source, given a flux scale and a distance scale.
Parameters:  source –
 scale –
Returns:

inject_signal
(scale)[source]¶ Randomly select simulated events from the Monte Carlo dataset to simulate a signal for each source. The source flux can be scaled by the scale parameter.
Parameters: scale – Ratio of Injected Flux to source flux. Returns: Set of signal events for the given IC Season.

select_mc_band
(source)[source]¶ For a given source, selects MC events within a declination band of width +/ 5 degrees that contains the source. Then returns the MC data subset containing only those MC events.
Parameters: source – Source to be simulated Returns: mc (cut): Simulated events which lie within the band Returns: omega: Solid Angle of the chosen band Returns: band_mask: The mask which removes events outside band

subclasses
= {'low_memory_injector': <class 'flarestack.core.injector.LowMemoryInjector'>}¶


class
flarestack.core.injector.
MockUnblindedInjector
(season, sources=nan, **kwargs)[source]¶ If the data is not really to be unblinded, then MockUnblindedInjector should be called. In this case, the create_dataset function simply returns one background scramble.

class
flarestack.core.injector.
TrueUnblindedInjector
(season, sources, **kwargs)[source]¶ If the data is unblinded, then UnblindedInjector should be called. In this case, the create_dataset function simply returns the unblinded dataset.
Log Likelihood¶

class
flarestack.core.llh.
FixedEnergyLLH
(season, sources, llh_dict)[source]¶ 
__init__
(season, sources, llh_dict)[source]¶ Initialize self. See help(type(self)) for accurate signature.

calculate_test_statistic
(params, weights, **kwargs)[source]¶ Calculates the test statistic, given the parameters. Uses numexpr for faster calculations.
Parameters:  params – Parameters from Minimisation
 weights – Normalised fraction of n_s allocated to each source
Returns: 2 * llh value (Equal to Test Statistic)

create_energy_functions
()[source]¶ Creates the acceptance function, which parameterises signal acceptance as a function of declination, and the energy weighting function, which gives the energy signaloverbackground ratio
Returns: Acceptance function, energy_weighting_function

create_kwargs
(data, pull_corrector, weight_f=None)[source]¶ Creates a likelihood function to minimise, based on the dataset.
Parameters: data – Dataset Returns: LLH function that can be minimised

fit_energy
= False¶


class
flarestack.core.llh.
LLH
(season, sources, llh_dict)[source]¶ Base class LLH.

__init__
(season, sources, llh_dict)[source]¶ Initialize self. See help(type(self)) for accurate signature.

static
assume_background
(n_s, n_coincident, n_all)[source]¶ To save time with likelihood calculation, it can be assumed that all events defined as “noncoincident”, because of distance in space and time to the source, are in fact background events. This is equivalent to setting S=0 for all noncoincident events. IN this case, the likelihood can be calculated as the product of the number of noncoincident events, and the likelihood of an event which has S=0.
Parameters:  n_s – Array of expected number of events
 n_coincident – Number of events that were not assumed to have S=0
 n_all – The total number of events
Returns: Log Likelihood value for the given

background_pdf
(source, cut_data)[source]¶ Calculates the value of the background spatial PDF for a given source for each event in the coincident data subsample. Thus is done by calling the self.bkg_spline spline function, which was fitted to the Sin(Declination) distribution of the data.
If there is a signal Time PDF given, then the background time PDF is also calculated for each event. This is assumed to be a normalised uniform distribution for the season.
Returns either the background spatial PDF values, or the product of the background spatial and time PDFs.
Parameters:  source – Source to be considered
 cut_data – Subset of Dataset with coincident events
Returns: Array of Background Spacetime PDF values

create_energy_functions
()[source]¶ Creates the acceptance function, which parameterises signal acceptance as a function of declination, and the energy weighting function, which gives the energy signaloverbackground ratio
Returns: Acceptance function, energy_weighting_function

create_llh_function
(data, pull_corrector, weight_f=None)[source]¶ Creates a likelihood function to minimise, based on the dataset.
Parameters: data – Dataset Returns: LLH function that can be minimised

classmethod
register_subclass
(llh_name)[source]¶ Adds a new subclass of EnergyPDF, with class name equal to “energy_pdf_name”.

select_spatially_coincident_data
(data, sources)[source]¶ Checks each source, and only identifies events in data which are both spatially and timecoincident with the source. Spatial coincidence is defined as a +/ 5 degree box centered on the given source. Time coincidence is determined by the parameters of the LLH Time PDF. Produces a mask for the dataset, which removes all events which are not coincident with at least one source.
Parameters:  data – Dataset to be tested
 sources – Sources to be tested
Returns: Mask to remove

signal_pdf
(source, cut_data)[source]¶ Calculates the value of the signal spatial PDF for a given source for each event in the coincident data subsample. If there is a Time PDF given, also calculates the value of the signal Time PDF for each event. Returns either the signal spatial PDF values, or the product of the signal spatial and time PDFs.
Parameters:  source – Source to be considered
 cut_data – Subset of Dataset with coincident events
Returns: Array of Signal Spacetime PDF values

subclasses
= {'fixed_energy': <class 'flarestack.core.llh.FixedEnergyLLH'>, 'spatial': <class 'flarestack.core.llh.SpatialLLH'>, 'standard': <class 'flarestack.core.llh.StandardLLH'>, 'standard_matrix': <class 'flarestack.core.llh.StandardMatrixLLH'>, 'standard_overlapping': <class 'flarestack.core.llh.StandardOverlappingLLH'>}¶


class
flarestack.core.llh.
SpatialLLH
(season, sources, llh_dict)[source]¶ Most basic LLH, in which only spatial, and optionally also temporal, information is included. No Energy PDF is used, and no energy weighting is applied.

__init__
(season, sources, llh_dict)[source]¶ Initialize self. See help(type(self)) for accurate signature.

calculate_test_statistic
(params, weights, **kwargs)[source]¶ Calculates the test statistic, given the parameters. Uses numexpr for faster calculations.
Parameters:  params – Parameters from Minimisation
 weights – Normalised fraction of n_s allocated to each source
Returns: 2 * llh value (Equal to Test Statistic)

create_energy_function
()[source]¶ In the most simple case of spatialonly weighting, you would neglect the energy weighting of events. Then, you can simply assume that the detector acceptance is roughly proportional to the data rate, i.e assuming that the incident background atmospheric neutrino flux is uniform. Thus the acceptance of the detector is simply the background spatial PDF (which is a spline fitted to data as a function of declination). This method does, admittedly neglect the fact that background in the southern hemisphere is mainly composed of muon bundles, rather than atmospheric neutrinos. Still, it’s slighty better than assuming a uniform detector acceptance
Returns: 1D linear interpolation

create_llh_function
(data, pull_corrector, weight_f=None)[source]¶ Creates a likelihood function to minimise, based on the dataset.
Parameters:  data – Dataset
 pull_corrector – pull_corrector
Returns: LLH function that can be minimised

fit_energy
= False¶


class
flarestack.core.llh.
StandardLLH
(season, sources, llh_dict)[source]¶ 
__init__
(season, sources, llh_dict)[source]¶ Initialize self. See help(type(self)) for accurate signature.

calculate_test_statistic
(params, weights, **kwargs)[source]¶ Calculates the test statistic, given the parameters. Uses numexpr for faster calculations.
Parameters:  params – Parameters from Minimisation
 weights – Normalised fraction of n_s allocated to each source
Returns: 2 * llh value (Equal to Test Statistic)

create_SoB_energy_cache
(cut_data)[source]¶ Evaluates the Log(Signal/Background) values for all coincident data. For each value of gamma in self.gamma_support_points, calculates the Log(Signal/Background) values for the coincident data. Then saves each weight array to a dictionary.
Parameters: cut_data – Subset of the data containing only coincident events Returns: Dictionary containing SoB values for each event for each gamma value.

create_acceptance_function
()[source]¶ Creates a 2D linear interpolation of the acceptance of the detector for the given season, as a function of declination and gamma. Returns this interpolation function.
Returns: 2D linear interpolation

create_energy_functions
()[source]¶ Creates the acceptance function, which parameterises signal acceptance as a function of declination, and the energy weighting function, which gives the energy signaloverbackground ratio
Returns: Acceptance function, energy_weighting_function

create_kwargs
(data, pull_corrector, weight_f=None)[source]¶ Creates a likelihood function to minimise, based on the dataset.
Parameters: data – Dataset Returns: LLH function that can be minimised

estimate_energy_weights
(gamma, energy_SoB_cache)[source]¶ Quickly estimates the value of Signal/Background for Gamma. Uses precalculated values for first and second derivatives. Uses a Taylor series to estimate S(gamma), unless SoB has already been calculated for a given gamma.
Parameters:  gamma – Spectral Index
 energy_SoB_cache – Weight cache
Returns: Estimated value for S(gamma)

fit_energy
= True¶

new_acceptance
(source, params=None)[source]¶ Calculates the detector acceptance for a given source, using the 2D interpolation of the acceptance as a function of declination and gamma. If gamma IS NOT being fit, uses the default value of gamma for weighting (determined in __init__). If gamma IS being fit, it will be the last entry in the parameter array, and is the acceptance uses this value.
Parameters:  source – Source to be considered
 params – Parameter array
Returns: Value for the acceptance of the detector, in the given
season, for the source


class
flarestack.core.llh.
StandardOverlappingLLH
(season, sources, llh_dict)[source]¶ 
calculate_test_statistic
(params, weights, **kwargs)[source]¶ Calculates the test statistic, given the parameters. Uses numexpr for faster calculations.
Parameters:  params – Parameters from Minimisation
 weights – Normalised fraction of n_s allocated to each source
Returns: 2 * llh value (Equal to Test Statistic)

Utils¶
IceCube Utils¶

flarestack.icecube_utils.dataset_loader.
data_loader
(data_path, floor=True, cut_fields=True)[source]¶ Helper function to load data for a given season/set of season. Adds sinDec field if this is not available, and combines multiple years of data is appropriate (different sets of data from the same icecube configuration should be given as a list)
Parameters:  data_path – Path to data or list of paths to data
 cut_fields – Boolean to remove unused fields from datasets on loading
Returns: Loaded Dataset (experimental or MC)