diff --git a/_PeakOTronAP2.py b/_PeakOTronAP2.py
index 2ac9e512eebdd3bfcfa2b6777dbf45b9f430e7ec..dd2c6601caf8df7acb06ca5f8ae4c121b792d8a7 100644
--- a/_PeakOTronAP2.py
+++ b/_PeakOTronAP2.py
@@ -3,7 +3,7 @@
import os
import numpy as np
import pandas as pd
-from itertools import chain
+from itertools import chain, zip_longest
from iminuit import Minuit
from iminuit.util import describe
import matplotlib.pyplot as plt
@@ -532,67 +532,22 @@ class PeakOTron:
return x_kde, y_kde, y_kde_err, y_kde_conf_interval
-
- ##MODEL DEFINITIONS
-
-
-
-# def DRM_basic(self, x, mu, x_0, G, sigma_0, sigma_1, alpha, r_beta, lbda, k_low, k_hi):
-
-# k = np.arange(k_low, k_hi)
-# k = sorted(list(k[k!=0]))
-# k_is = [k[0:int(_k)] for _k in k]
-
-# f0 = gpoisson.pmf(0, mu, lbda)*norm.pdf(x, x_0, sigma_0)
-# f1s = np.asarray([
-# binom.pmf(0, _k, alpha)*gpoisson.pmf(_k, mu, lbda)*norm.pdf(x, x_0 + _k*G, self.SigmaK(_k, sigma_0, sigma_1))
-# for _k in k
-# ])
-
-# f2s = np.vstack([
-# np.asarray([
-# binom.pmf(_i, _k, alpha)*gpoisson.pmf(_k, mu, lbda)*self.XTM(x, _k, _i, x_0, G, sigma_0, sigma_1, r_beta)
-# for _i in _is
-# ])
-
-# for _k, _is in zip(k, k_is)
-# ])
-
-
-# f = np.sum(np.vstack([f0, f1s, f2s]), axis=0)
-
-# return f
-# def XTM(self, x, k, i, x_0, G, sigma_0, sigma_1, r_beta):
-
-# x_k = x - (x_0 + k*G)
-# _sigma_k = self.SigmaK(k, sigma_0, sigma_1)
-# x_ct_pad = np.zeros_like(x_k)
-
-# if(i==1):
-# cond = (x_k > -5*_sigma_k)
-# error_f = 0.5*(1.0 + sc.erf(x_k/(_sigma_k*np.sqrt(2.0))))
-# log_f = -x_k/(r_beta*G) - np.log(r_beta*G)
-# x_ct = np.where(cond, np.exp(log_f)*error_f, x_ct_pad)
-
-# elif(i>1):
-# cond = (x_k > 0)
-# log_f = sc.xlogy(i-1, x_k) - sc.gammaln(i) - sc.xlogy(i, r_beta*G) - x_k/(r_beta*G)
-# x_ct = np.where(cond, np.exp(log_f), x_ct_pad)
-
-# else:
-# x_ct = x_k
-
-# return x_ct
-
-
-
- def DRM_basic(self, x, mu, x_0, G, sigma_0, sigma_1, lbda, tau, r_Ap, tgate, pAp, k_low, k_hi):
-
+ def DRM_basic(self,
+ x,
+ mu,
+ x_0,
+ G,
+ lbda,
+ sigma_0,
+ sigma_1,
+ k_low,
+ k_hi):
+
k = np.arange(k_low, k_hi)
k = sorted(list(k[k!=0]))
k_is = [k[0:int(_k)] for _k in k]
@@ -603,67 +558,10 @@ class PeakOTron:
for _k in k
])
- f2s = np.asarray([
- gpoisson.pmf(_k, mu, lbda)*_k*self.dPApdPH(x,
- _k,
- x_0,
- G,
- tau,
- r_Ap,
- tgate,
- pAp)
- for _k in k
- ])
- f = np.sum(np.vstack([f0, f1s, f2s]), axis=0)
+ f = np.sum(np.vstack([f0, f1s]), axis=0)
return f
-
-
-
-
-
-# def dPApdPH(
-# self,
-# x,
-# k,
-# x_0,
-# G,
-# tau,
-# r_Ap,
-# tgate,
-# pAp):
-
-# PH = (x - x_0)/G - k
-
-# PH_max = 1 + np.exp(-tgate/tau) - 2*np.exp(-0.5*tgate/tau)
-
-# cond_PH_max = (PH<PH_max)
-
-# PH_temp = np.where(cond_PH_max, PH, PH_max)
-
-# t_b0 = 0.5*tgate - tau*np.arccosh(0.5*((1 - PH_temp)*np.exp(tgate/(2*tau)) + np.exp(-tgate/(2*tau))))
-# t_b1 = tgate-t_b0
-
-# dpApdt_b0 = pAp*(1-np.exp(-t_b0/tau))*np.exp(-t_b0/r_Ap)/r_Ap
-# dpApdt_b1 = pAp*(1-np.exp(-t_b1/tau))*np.exp(-t_b1/r_Ap)/r_Ap
-
-# dPHdt_b0 = -2*np.sinh((t_b0-0.5*tgate)/tau)*np.exp(-0.5*tgate/tau)/tau
-# dPHdt_b1 = -2*np.sinh((t_b1-0.5*tgate)/tau)*np.exp(-0.5*tgate/tau)/tau
-
-# dpApdPH_b0 = PH_max*dpApdt_b0/((np.where(np.abs(dPHdt_b0)>1e-5, np.abs(dPHdt_b0), 1e-5)))
-# dpApdPH_b1 = PH_max*dpApdt_b1/((np.where(np.abs(dPHdt_b1)>1e-5, np.abs(dPHdt_b1), 1e-5)))
-
-
-
-# dpApdPH_b0 = np.where(cond_PH_max, dpApdPH_b0, 0)
-# dpApdPH_b1 = np.where(cond_PH_max, dpApdPH_b1, 0)
-
-# dpApdPH_b0 = np.where(dpApdPH_b0>0, dpApdPH_b0, 0)
-# dpApdPH_b1 = np.where(dpApdPH_b0>0, dpApdPH_b1, 0)
-
-
-# return dpApdPH_b0/G
def dPApdPH(self,
@@ -707,11 +605,56 @@ class PeakOTron:
dpApdPH_b1 = np.where(dpApdPH_b0>0, dpApdPH_b1, 0)
- return dpApdPH_b0/G
+ return PH_max*(dpApdPH_b0)/G
+
+
+
+ def ApplyAp(self,
+ x,
+ y_light,
+ mu,
+ x_0,
+ G,
+ lbda,
+ DCR,
+ tgate,
+ t_0,
+ tau,
+ r_Ap,
+ pAp,
+ k_low,
+ k_hi,
+ k_dcr_low,
+ k_dcr_hi):
+
+ k = np.arange(k_low, k_hi)
+ k = sorted(list(k[k!=0]))
+
+ mu_dcr = DCR*(tgate+abs(t_0))
+
+ k_dcr = np.arange(k_dcr_low, k_dcr_hi)
+ k_dcr = sorted(list(k_dcr[k_dcr!=0]))
+
+
+ ps_light = [gpoisson.pmf(_k, mu, lbda) for _k in k]
+ ps_dark = [gpoisson.pmf(_k, mu_dcr, lbda) for _k in k_dcr]
+
+ f0s = np.asarray([_p_l + _p_d for _p_l, _p_d in zip_longest(ps_light, ps_dark, fillvalue=0)]).reshape(-1,1)
+
+
+ f1s = np.asarray([
+ _k*self.dPApdPH(x, _k, x_0, G, tau, r_Ap, tgate, pAp)
+ for _k in k
+ ])
+
+ f = np.sum(f0s*f1s, axis=0)
+
+ return y_light+f
+
def PH_range(self, t, t_0, r_fast, tau, tgate):
@@ -748,7 +691,20 @@ class PeakOTron:
return PHs
- def ApplyDCR(self, x, y_light, x_0, G, DCR, tgate, t_0, tau, r_fast, lbda, dx, k_dcr_low, k_dcr_hi):
+ def ApplyDCR(self,
+ x,
+ y_light,
+ x_0,
+ G,
+ lbda,
+ DCR,
+ tgate,
+ t_0,
+ tau,
+ r_fast,
+ dx,
+ k_dcr_low,
+ k_dcr_hi):
mu_dcr = DCR*(abs(t_0) + tgate)
@@ -836,36 +792,7 @@ class PeakOTron:
- def DRM_nodcr(
- self,
- x,
- mu,
- x_0,
- G,
- sigma_0,
- sigma_1,
- alpha,
- r_beta,
- lbda,
- k_low,
- k_hi):
-
- y_model = self.DRM_basic(x,
- mu,
- x_0,
- G,
- sigma_0,
- sigma_1,
- alpha,
- r_beta,
- lbda,
- k_low,
- k_hi,
- dx
- )
-
- return y_model
-
+
def DRM(
@@ -894,33 +821,49 @@ class PeakOTron:
y_light = self.DRM_basic(x,
- mu,
+ mu,
x_0,
G,
+ lbda,
sigma_0,
sigma_1,
- lbda,
- tau,
- r_Ap,
- tgate,
- pAp,
k_low,
k_hi)
-
- y_model = self.ApplyDCR(x,
- y_light,
- x_0,
- G,
- DCR,
- tgate,
- t_0,
- tau,
- r_fast,
- lbda,
- dx,
- k_dcr_low,
- k_dcr_hi)
+ y_lightap = self.ApplyAp(x,
+ y_light,
+ mu,
+ x_0,
+ G,
+ lbda,
+ DCR,
+ tgate,
+ t_0,
+ tau,
+ r_Ap,
+ pAp,
+ k_low,
+ k_hi,
+ k_dcr_low,
+ k_dcr_hi)
+
+
+
+ y_model = self.ApplyDCR(x,
+ y_lightap,
+ x_0,
+ G,
+ lbda,
+ DCR,
+ tgate,
+ t_0,
+ tau,
+ r_fast,
+ dx,
+ k_dcr_low,
+ k_dcr_hi)
+
+
return y_model