The spacecraft information is a list, for 30-s intervals, of the livetime and instrument $\theta$ of the source. We need a weighted effective area as a function of $\theta$. This requires a weighting function of energy, or a spectrum. Matthew used a simple power law with index -2.1. implemented in by the class KerrAeff. But here we have the actual source spectrum available. Also, the data here excludes low-energy Back events the class SourceAeff uses both.

Implemented here:

• time_bin_edges
• cell_exposure
• sc_data_selection
• binned_exposure

time_bin_edges[source]

time_bin_edges(config, exposure, tbin=None)

Return an interleaved array of start/stop values

tbin: an array (a,b,d), default config.time_bins

interpretation of a, b:

if > 50000: interpret as MJD
if < 0:     offset from stop
otherwise:  offset from start

d : if positive, the day bin size if 0; return contiguous bins

cell_exposure[source]

cell_exposure(config, exposure, time_bins, exposure_factor=1e-06, set_costh=False)

Determine the exposure per cell input:

• exposure -- dataframe with spacecraft data
• time_bins -- time binning (start, stop, step) tuple

Return: a dict of

• edges -- list of cell edges, an interleaved array of start/stop values
• exp -- exposure per cell, in cm^2 Ms units, if exposure_factor==1e-6
• etot -- total exposure
• costh -- mean cos theta per cell if requested
• exp_energy if exp_fract in the exposure DF, set exposure energy

Note: total exposure vs. energy is: t = df.apply(lambda r: np.array(r.exp * np.array(r.exp_fract), float), axis=1).values u = np.vstack(t)

cumsimpson[source]

cumsimpson(y, dx)

Return the cumulative Simpson integral for an odd number>2 of equally-spaced evaluations

• y -- array of integrand values, 2n+1 values for 2n intervals
• dx -- interval size

Returns array of size n, the cumulative integral estimate for pairs of intervals.

x, dx = np.linspace(0,1, 5), 1/4 print(cumsimpson(x, dx)) [0.125 0.5 ]

Demonstrate superiority of integral using integrand $x\ f(x)\ d\log(x)$

Instead of $f(x)\ dx$.

nbins = 16
f = lambda x: x**-2
ee = np.logspace(2,4,nbins+1)
r = 1e-2-1e-4
print(f'Errors:\n\tLinear: {simpson(f(ee), ee)/r-1:+6.1e}'
      f'\n\tLog:    {simpson( f(ee)*ee, np.log(ee) )/r-1:+6.1e}')

Errors:
	Linear: -4.7e-03
	Log:    +3.8e-05

class BaseAeff[source]

BaseAeff(config, source, elim:tuple=(100.0, 10000.0))

Base class for calculating the weighted effective area

It defines a functor of costh, returning the effective area weighted by a spectrum. The exposure is then a sum over intervals of this weighted effective area times the livetime per interval.

A subclass must define a setup method to define the spectral function and minimum energy to include Back events.

class KerrAeff[source]

KerrAeff(config, source, elim:tuple=(100.0, 10000.0)) :: BaseAeff

Kerr implementation from godot, which uses a $E^{-2.1}$ powerlaw spectrum

class SourceAeff[source]

SourceAeff(config, source, elim:tuple=(100.0, 10000.0)) :: BaseAeff

BaseAeff subclass that uses the actual source spectrum applied only to used bands

class WeightedAeff[source]

WeightedAeff(config=None, spectrum=None, nebins:int=16, imin:tuple=(0, 0), **kwargs)

Manage the weighted effective area calculation

class NewSourceAeff[source]

NewSourceAeff(config, source) :: WeightedAeff

Manage the weighted effective area calculation

Weighted exposure

The computation of the likelihood for a given cell needs a measurement of the exposure to determine the estimated number of source counts expected for the average source rate.

weighted_aeff[source]

weighted_aeff(config, source)

Use config.use_kerr to return appropriate weighted effective area function of cos theta

Compare weighted effective area: Kerr-mode vs actual with Geminga

from wtlike.sources import PointSource
config=Config()
# source_name='LS 5039' #'Geminga'

def aeffplots(ke,se):
    ct = np.linspace(0.4, 1.0, 25)

    fig, (ax1,ax2) = plt.subplots(1,2, figsize=(8,4), sharey=True)
    plt.subplots_adjust(wspace=0.05)
    ax1.semilogx(ke.edom, ke.wts/ke.wts[0], '.-', label='power law')
    ax1.semilogx(se.edom, se.wts/se.wts[0], '.-', label=f'{source_name} spectrum')
    ax1.legend();ax1.grid(alpha=0.5)
    ax1.set(title='energy dependence', ylabel='relative weight')

    ax2.plot(ct, se(ct)/se(1), '.-', label='power law')
    ax2.plot(ct, ke(ct)/ke(1), '.-', label=f'{source_name} spectrum')
    ax2.set( title='angular dependence',xlim=(0.4,1.0), ylim=(0, None),)
    ax2.grid(alpha=0.5); ax2.legend();

if config.valid:
    source_name='Geminga'
    source=PointSource(source_name)

    ke = KerrAeff(config, source) 
    se = SourceAeff(config, source)
    aeffplots(ke,se)

sc_data_selection[source]

sc_data_selection(config, source, sc_data)

Return a DataFrame with the S/C data for the source direction, wtih cos theta and zenith cuts

columns:

• start, stop, livetime -- from the FT2 info
• cos_theta -- angle between bore and direction
• exp -- the exposure: effective area at angle weighted by a default spectral function, times livetime

binned_exposure[source]

binned_exposure(config, exposure, time_edges)

Bin the exposure into cells

• exposure -- A DataFrame derived from FT2
• time_bins: list of edges, an interleaved start/stop array

Returns:

An array of exposure integrated over each time bin. Assumes that the time bins are contained within the exposure.

it is interleaved, client must apply [0::2] selection.