Gaussian and quadratic example functions

n = 20
sf = _Sampler(lambda x: np.exp(-(x**2)/2), limits=(-4, 4) )

data = sf(10000)
tests = np.array([np.abs(data.mean()), np.abs(data.std()-1) ])
assert np.all(tests<5e-2 ), f'Failed Tests: mean {data.mean()}, std {data.std()}'

func = lambda x: x**2
wfun = _Sampler(func)

test2 = wfun.mean,  np.mean(wfun(1000))
assert np.abs( test2[0]-test2[1] ) < 1e-1, f'Not almost equal: {test2}'

Test generating weights from a source weight histogram

Test with a function peaked at both ends, generate equal signal and background

Simulate times and weights

generate_times[source]

generate_times(start, stop, count, rng=None)

Generate a list of times, distributed randomly

• start, stop: times
• count : expected number to generate with rate=count/(stop-start)
• rng : random state generator

returns : list of times between start and stop. Note that the actual number is Poisson-distributed

rng = np.random.default_rng(42)
nn = np.array( [ len(generate_times(0,1,10, rng=rng )) for i in range(100000) ])
plt.hist(nn, np.linspace(0,25,26), label=f'Mean: {nn.mean():.2f},\n  std: {nn.std():.2f}');
plt.legend()
plt.title(f'Check that generated number is Poisson', fontsize=12);

class WeightFunction[source]

WeightFunction(s=1, b=1, wt_signif=0.1, rng=None)

# test weight function
s, lam = 1,0.1

def wfuntest(s, b=1, lam=0.1, Nb=10000, rng=42):
   
    wfun = WeightFunction(s=s,b=b,wt_signif=lam, rng=rng)
    n = int(Nb*(1+s/b)) 
    samp = wfun.sample(s, n ) 
    wts = wfun.weights(s, n)
    sumwts = np.sum(wts) #wfun.weights(s,b, 10000)

    fig, (ax1, ax3) = plt.subplots(1,2, figsize=(7,3), sharex=True)
    plt.subplots_adjust(top=0.80, wspace=0.3)
    
    r = np.linspace(0,1, 1000)  
    ax1.plot(r, wfun(r)/(s+b), '-',color='orange', lw=2,label='PSF+bkg')
    ax1.set(ylim=(0,None), xlim=(0,1), xlabel='$r^2$', title='position')
    
    ax1.hist(samp, np.linspace(0,1,51), density=True, color='cornflowerblue',
             label=f'generated:\n mean:{np.mean(samp):.2f}\n  cnt: {len(samp)}');
    ax1.legend(fontsize=10)
    #ax1.set(title='sampled', xlabel='$r^2$'); ax1.legend()
    
    
    ax3.hist(wts, bins=np.linspace(0,1,51), label=f'sum:{sumwts:.1f}')
    ax3.set(title='weights', label='weight'); ax3.legend(fontsize=10)
    fig.suptitle(f'WeightFunction test: s={s}')
wfuntest(s)

make_exposure[source]

make_exposure(fexp, start, stop, interval=300, costh=0.8)

• fexp -- exposure in cm^2, a value or a function of time in day units
• start, stop -- range of time in day units
• cos_theta -- just set to costh
• interval [300] -- 5-min interval (fermi data is 30 s)

Returns: a DataFrame with start, stop, exp

def src_flare(t, tzero=5, width=0.25, amp=2):
    return 1e-6*(1 + amp*np.exp(-(t-tzero)**2/2/width))
sim = Simulation('test_sim', src_flux=src_flare, tstart=0, tstop=10, rng=42 , debug=10)
sim.run()
sim.photons.hist(bins=50, color='cornflowerblue');
sim.exposure.describe()

time: 0.000 - 0.003, source 1e-06, exposure/s 3000, expected/generated counts 2/0
	 weights []
time: 0.003 - 0.007, source 1e-06, exposure/s 3000, expected/generated counts 2/0
	 weights []
time: 0.007 - 0.010, source 1e-06, exposure/s 3000, expected/generated counts 2/0
	 weights []
time: 0.010 - 0.014, source 1e-06, exposure/s 3000, expected/generated counts 2/2
	 weights [0.05 0.03]
time: 0.014 - 0.017, source 1e-06, exposure/s 3000, expected/generated counts 2/3
	 weights [0.28 0.02 0.72]
time: 0.017 - 0.021, source 1e-06, exposure/s 3000, expected/generated counts 2/4
	 weights [0.05 0.79 0.   0.  ]
time: 0.021 - 0.024, source 1e-06, exposure/s 3000, expected/generated counts 2/5
	 weights [0.07 0.   0.81 0.78 0.68]
time: 0.024 - 0.028, source 1e-06, exposure/s 3000, expected/generated counts 2/3
	 weights [0.21 0.72 0.01]
time: 0.028 - 0.031, source 1e-06, exposure/s 3000, expected/generated counts 2/2
	 weights [0.02 0.77]
time: 0.031 - 0.035, source 1e-06, exposure/s 3000, expected/generated counts 2/6
	 weights [0.13 0.03 0.7  0.49 0.88 0.89]
generated 5895 photons

	start	stop	cos_theta	exp
count	2880.000000	2880.000000	2880.000000	2880.0
mean	4.998264	5.001736	0.800023	900000.0
std	2.887252	2.887252	0.000023	0.0
min	0.000000	0.003472	0.800000	900000.0
25%	2.499132	2.502604	0.800000	900000.0
50%	4.998264	5.001736	0.800000	900000.0
75%	7.497396	7.500868	0.800000	900000.0
max	9.996528	10.000000	0.800000	900000.0

class Simulation[source]

Simulation(name, src_flux, tstart, tstop, bkg_rate=1e-06, efun=3000, wt_signif=0.1, debug=False, rng=None, config=None)

• src_flux : source flux, scalar or function of days, typically around 1e-7
• tstart, tstop :(days)
• bkg_rate : background flux, scalar or function of day, typicaly 1e-6 for 4-deg cone

• efun : scalar | function (of time in days) of the exposure/s. Typically 3000 cm^2 for fermi

• wt_signif : the width of the PSF in (r/rmax)**2 coordinates

• rng : random generator instance, or integer seed

class WtGen[source]

WtGen(src, aeff=<lambda>) :: dict

Encapsulate the dict of weights used by wtlike for simulation

It is a functor of a list of uniform probabilities, which retuns a DataFrame with columns bandid, weight, pixel index

Test by generating with Geminga

Check SED as measued with sums of weights for each band with input SED

%time self.setup()#alpha = lambda n: 0)

CPU times: user 4.65 ms, sys: 578 µs, total: 5.23 ms
Wall time: 3.16 ms