Skip to content
Snippets Groups Projects
Commit 22e97e90 authored by Tobias Quadfasel's avatar Tobias Quadfasel
Browse files

Some minor bug fixes and re-formatting to PEP8 style. Double-checked that... 

Some minor bug fixes and re-formatting to PEP8 style. Double-checked that reformatting does not influence results by comparing output of
parent 8d666a3f
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
AdditionalPDFs.py
+ 47
− 48
View file @ 22e97e90
import numpy as np
from scipy.stats import rv_discrete, rv_continuous, uniform
from scipy.stats import uniform
import scipy.special as sc
import matplotlib.pyplot as plt
import scipy
from scipy.stats._distn_infrastructure import (
rv_discrete, _ncx2_pdf, _ncx2_cdf, get_distribution_names)
rv_discrete, get_distribution_names)
class gpd_gen(rv_discrete):
......@@ -34,7 +32,8 @@ class gpd_gen(rv_discrete):
def _logpmf(self, k, mu, lbda):
mu_pls_klmb = mu + lbda*k
return np.log(mu) + sc.xlogy(k-1, mu_pls_klmb) - mu_pls_klmb - sc.gammaln(k+1)
return (np.log(mu) + sc.xlogy(k-1, mu_pls_klmb)
- mu_pls_klmb - sc.gammaln(k+1))
def _munp(self, n, mu, lbda):
if n == 1:
......@@ -83,44 +82,44 @@ class borel_gen(rv_discrete):
return _rnd
def _stats(self, mu):
_mu = 1/(1-mu)
_var = mu/(1-mu)**3
tmp = np.asarray(mu)
mu_nonzero = ((tmp > 0) & (tmp < 1))
#g1 and g2: Lagrangian Probability Distributions, 978-0-8176-4365-2, page 159
g1 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,), lambda x: (1+2*x)/scipy.sqrt(x*(1-x)), np.inf)
g2 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,), lambda x: 3 + (1 + 8*x+6*x**2)/(x*(1-x)), np.inf)
return _mu, _var, g1, g2
# g1 and g2: Lagrangian Probability Distributions,
# 978-0-8176-4365-2, page 159
g1 = scipy._lib._util._lazywhere(mu_nonzero, (tmp,),
lambda x: (1+2*x)/scipy.sqrt(x*(1-x)),
np.inf)
borel= borel_gen(name='borel')
g2 = scipy._lib._util._lazywhere(
mu_nonzero, (tmp,),
lambda x: 3 + (1 + 8*x+6*x**2)/(x*(1-x)), np.inf)
return _mu, _var, g1, g2
class erlang_gen(rv_discrete):
borel = borel_gen(name='borel')
class erlang_gen(rv_discrete):
def _pdf(self, x, a):
# gamma.pdf(x, a) = x**(a-1) * exp(-x) / gamma(a)
return np.exp(self._logpdf(x, a))
def _logpdf(self, k, mu, nu):
# TODO: x and a not defined!
return sc.xlogy(a - 1.0, x) - x - sc.gammaln(a)
# def _rvs(self, mu, nu, size=None, random_state=None):
# u = scipy.stats.uniform.rvs(loc=0, scale = 1, size=size)
# cum = np.cumsum([self._pmf(_k, mu, nu) for _k in range(0, 100)])
# rnd = [ np.argmax( cum>=_u ) for _u in u ]
# return rnd
pairs = list(globals().items())
_distn_names, _distn_gen_names = get_distribution_names(pairs, rv_discrete)
......
HelperFunctions.py
+ 38
− 41
View file @ 22e97e90
from iminuit import Minuit
from iminuit.util import describe
from scipy.interpolate import interp1d
from astropy.stats import bootstrap
import numpy as np
import matplotlib.pyplot as plt
from scipy import special as sc
#HELPER FUNCTIONS
class FakeFuncCode:
def __init__(self, f, prmt=None, dock=0, append=None):
# f can either be tuple or function object
......@@ -20,7 +17,8 @@ class FakeFuncCode:
for i, p in enumerate(prmt):
self.co_varnames[i] = p
if isinstance(append, str): append = [append]
if isinstance(append, str):
append = [append]
if append is not None:
old_count = self.co_argcount
......@@ -31,8 +29,6 @@ class FakeFuncCode:
list(self.co_varnames[old_count:]))
def LatexFormat(value, scirange=[0.01, 1000]):
if(np.abs(value) > scirange[0] and np.abs(value) < scirange[1]):
float_str = r"${:3.3f}$".format(value)
......@@ -42,6 +38,7 @@ def LatexFormat(value, scirange=[0.01,1000]):
base, exponent = float_str.split("E")
float_str = r"${0} \times 10^{{{1}}}$".format(base, int(exponent))
except:
# TODO: Bare except! Is this a ValueError?
float_str = str(value)
return float_str
......@@ -67,11 +64,11 @@ def FormatExponent(ax, axis='y'):
# Run plt.tight_layout() because otherwise the offset text doesn't update
plt.tight_layout()
##### THIS IS A BUG
##### Well, at least it's sub-optimal because you might not
##### want to use tight_layout(). If anyone has a better way of
##### ensuring the offset text is updated appropriately
##### please comment!
# THIS IS A BUG
# Well, at least it's sub-optimal because you might not
# want to use tight_layout(). If anyone has a better way of
# ensuring the offset text is updated appropriately
# please comment!
# Get the offset value
offset = ax_axis.get_offset_text().get_text()
......@@ -96,11 +93,9 @@ def Bootstrap(data, statistic, n_bootstrap, alpha=0.95):
if not (0 < alpha < 1):
raise ValueError("confidence level must be in (0, 1)")
if len(data) < 1:
raise ValueError("data must contain at least one measurement.")
boot_stat = bootstrap(data, n_bootstrap, bootfunc=statistic)
stat_data = statistic(data)
......@@ -110,21 +105,23 @@ def Bootstrap(data, statistic, n_bootstrap, alpha=0.95):
z_score = np.sqrt(2.0)*sc.erfinv(alpha)
conf_interval = est_stat + z_score*np.array((-std_err, std_err))
return est_stat, std_err, conf_interval
def Linear(x, m, c):
return m*x + c
def HistMoment(counts, bin_centres, i, mu):
return np.sum(counts*(bin_centres-mu)**i)
def HistMean(counts, bin_centres):
N = np.sum(counts)
mu = HistMoment(counts, bin_centres, 1, 0)/N
return mu
def HistStd(counts, bin_centres):
N = np.sum(counts)
mu = HistMean(counts, bin_centres)
......@@ -132,34 +129,34 @@ def HistStd(counts, bin_centres):
sigma = np.sqrt(var)
return sigma
def HistSkew(counts, bin_centres):
N = np.sum(counts)
mu = HistMean(counts, bin_centres)
gamma = (N*np.sqrt((N-1))/(N-2))*HistMoment(counts, bin_centres, 3, mu)/(HistMoment(counts, bin_centres, 2, mu)**(1.5))
gamma = ((N*np.sqrt((N-1))/(N-2))*HistMoment(counts, bin_centres, 3, mu)
/ (HistMoment(counts, bin_centres, 2, mu)**(1.5)))
return gamma
def f_tau(v, v_breakdown, v_0):
c_tau = (v - v_breakdown)/v_0
result = -1/np.log((1-np.exp(c_tau*np.exp(-1)))/(1 - np.exp(c_tau)))
return result
def GP_lbda(mu, sigma, gamma):
lbda = 0.5*(((mu*gamma)/(sigma)) - 1)
return lbda
def GP_gain(mu, sigma, gamma):
lbda = GP_lbda(mu, sigma, gamma)
gain = (sigma**2/mu)*((1 - lbda)**2)
return gain
def GP_mu(mu, sigma, gamma):
lbda = GP_lbda(mu, sigma, gamma)
mu_gp = (1/(1-lbda))*(mu**2/sigma**2)
return mu_gp
LossFunctions.py
+ 29
− 51
View file @ 22e97e90
import numpy as np
from scipy.interpolate import interp1d
from iminuit.util import describe
from Model import epsilon
from HelperFunctions import FakeFuncCode
import matplotlib.pyplot as plt
class BinnedLH:
def __init__(self, f, bin_centres, counts, bw):
self.f = f
self.x = bin_centres
......@@ -25,38 +19,28 @@ class BinnedLH:
self.eps = epsilon()
self.mask = (self.counts > 0)
def __call__(self, *arg):
self.last_arg = arg
y_hat = self.f(self.x, *arg)
E = y_hat*self.dx
h = self.counts
mask = self.mask
nlogL = np.zeros(self.len_x)
nlogL[mask] = h[mask]*(np.log(E[mask]) - np.log(h[mask])) + (h[mask]-E[mask])
nlogL[~mask] = -E[~mask]
nlogL[mask] = (h[mask]*(np.log(E[mask]) - np.log(h[mask]))
+ (h[mask]-E[mask]))
nlogL[~mask] = -E[~mask]
nlogL = -np.sum(nlogL)
self.n_calls += 1
return nlogL
class Chi2Regression:
def __init__(self, f, x, y, y_err, epsilon=1.35):
self.f = f
self.x = x
......@@ -68,23 +52,17 @@ class Chi2Regression:
self.func_code = FakeFuncCode(f, dock=True)
self.ndof = len(self.y) - (self.func_code.co_argcount - 1)
def __call__(self, *arg):
self.last_arg = arg
y_hat = self.f(self.x, *arg)
loss = ((y_hat - self.y)/(self.y_err))**2
return np.sum(loss)
class HuberRegression:
def __init__(self, f, x, y, y_err, delta=1.345):
self.f = f
self.x = x
......@@ -97,7 +75,6 @@ class HuberRegression:
self.func_code = FakeFuncCode(f, dock=True)
self.ndof = len(self.y) - (self.func_code.co_argcount - 1)
def __call__(self, *arg):
self.last_arg = arg
......@@ -105,6 +82,7 @@ class HuberRegression:
a = abs((self.y - self.f(self.x, *arg))/self.y_err)
cond_flag = (a > self.delta)
loss = np.sum((~cond_flag) * (0.5 * a ** 2) - (cond_flag) * self.delta * (0.5 * self.delta - a), -1)
loss = np.sum(((~cond_flag) * (0.5 * a ** 2) - (cond_flag)
* self.delta * (0.5 * self.delta - a)), -1)
return np.sum(loss)
Model.py
+ 660
− 595
View file @ 22e97e90
This diff is collapsed.
PeakOTron.py
+ 1448
− 1464
View file @ 22e97e90
This diff is collapsed.
example.py
+ 37
− 59
View file @ 22e97e90
import os
import sys
import numpy as np
import pandas as pd
from PeakOTron import PeakOTron
from joblib import dump
import time
from HelperFunctions import f_tau
C_tau = lambda V, V_bd, V_0: (V - V_bd)/V_0
f_tau = lambda V, V_bd, V_0: -1/np.log((1-np.exp(C_tau(V, V_bd, V_0)*np.exp(-1)))/(1 - np.exp(C_tau(V, V_bd, V_0))))
V_bd_hmt = 51.570574 + 0.307
V_bd_hmt = 51.877574
V_0_hmt = 2.906
tau = 21.953 ##SLOW COMPONENT OF SIPM PULSE
t_0 = 100.0 ## PRE-INTEGRATION TIME
t_gate = 104.0 ## GATE LENGTH
bin_0=-100.0 ## SELECT FIRST BIN OF SPECTRUM (CAN BE AUTOMATIC)
truncate_nsigma0_up = 2.0 ## SCAN SPECTRUM FROM Q < Q_0 - 4 sigma_0
truncate_nsigma0_do = 2.0 ## EVALUATE SPECTRUM CHI2 IN Q_0 - x*sigma_0 < Q < Q_0 + 2*sigma_0
prefit_only=False ## FIT THE WHOLE SPECTRUM
tau = 21.953 # SLOW COMPONENT OF SIPM PULSE
t_0 = 100.0 # PRE-INTEGRATION TIME
t_gate = 104.0 # GATE LENGTH
bin_0 = -100.0 # SELECT FIRST BIN OF SPECTRUM (CAN BE AUTOMATIC)
truncate_nsigma0_up = 2.0 # SCAN SPECTRUM FROM Q < Q_0 - 4 sigma_0
# EVALUATE SPECTRUM CHI2 IN Q_0 - x*sigma_0 < Q < Q_0 + 2*sigma_0
truncate_nsigma0_do = 2.0
prefit_only = False # FIT THE WHOLE SPECTRUM
print("--------------------")
print('EXAMPLE SIPM CALIBRATION RUN')
print("--------------------")
out_dict = {}
files_to_fit = []
## Find all histograms in directory
# Find all histograms in directory
for root, dirs, files in os.walk("./data/hamamatsu_pcb6/Light"):
for file in files:
if file.endswith(".txt"):
files_to_fit.append([file, os.path.join(root, file)])
## Print files.
# Print files.
print("Files to fit:")
for i, (file, _) in enumerate(files_to_fit):
print('File {0}: {1}'.format(i, file))
SiPM = "PM1125NS_SBO"
## Loop thorough files
# Loop thorough files
for i, (file, path) in enumerate(files_to_fit):
items = file.split('_')
......@@ -59,21 +46,21 @@ for i, (file, path) in enumerate(files_to_fit):
f_tau_hmt = f_tau(V, V_bd_hmt, V_0_hmt)
print("\n\n")
print("===============================================================")
print("FIT {:d} - {:s}".format(i, file))
print("===============================================================")
print("\n\n")
## Load files.
# Load files.
data = np.loadtxt(path, skiprows=8)
## Create a PeakOTron Fit Object.
# Create a PeakOTron Fit Object.
f_data = PeakOTron(verbose=False)
## Perform fit.
# Perform fit.
# NOTE: BINNING RULE "knuth", "freedman", "scott"
# - use bw= #### to override. It is still required for DCR calculation.
f_data.Fit(data,
tau=tau, # SLOW PULSE COMPONENT TIME CONSTANT (ns)
t_gate=t_gate, # GATE LENGTH (ns)
......@@ -81,19 +68,13 @@ for i, (file, path) in enumerate(files_to_fit):
tau_R=f_tau_hmt*tau,
bin_0=bin_0,
truncate_nsigma0_up=truncate_nsigma0_up,
truncate_nsigma0_do = truncate_nsigma0_do
) #BINNING RULE "knuth", "freedman", "scott" - use bw= #### to override. it is still required for DCR calculation.
f_data.PlotFit(plot_in_bins=True, display=False, save_directory="./Results/{0}_fit.png".format(os.path.splitext(file)[0]))
dump(f_data, "./Results/{0}.joblib".format(os.path.splitext(file)[0]))
truncate_nsigma0_do=truncate_nsigma0_do)
f_data.PlotFit(
plot_in_bins=True, display=False,
save_directory=f"./Results/{os.path.splitext(file)[0]}_fit.png")
dump(f_data, f"./Results/{os.path.splitext(file)[0]}.joblib")
fit_out = {}
fit_val, fit_err = f_data.GetFitResults()
......@@ -118,8 +99,5 @@ for i, (file, path) in enumerate(files_to_fit):
print("===============================================================")
print("\n\n")
df = pd.DataFrame.from_dict(out_dict)
df.to_csv("./fit_results_{:s}.csv".format(SiPM))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment