def distr_func_lppl_pileup(gamma, gamma0_log_parab, gamma_inj, r, s, gamma_eq, ratio_pile_up , alpha, gamma_cut_acc):
b = np.zeros(gamma.shape)
a = np.zeros(gamma.shape)
m = gamma < gamma_inj
s1=s+0.5
b[m] = np.power(gamma[m]/gamma0_log_parab, s1)
f1 = np.zeros(gamma.shape)
m1 = np.logical_and(gamma < gamma0_log_parab,gamma>gamma_inj)
f1[m1] = np.power(gamma[m1]/gamma0_log_parab, -s)
f1[~m1] = np.power(gamma[~m1]/gamma0_log_parab, (-s -r*np.log10(gamma[~m1] / gamma0_log_parab )))
b[~m] = np.power(gamma_inj/gamma0_log_parab,s1+s)*f1[~m]*np.exp(-np.power((gamma[~m]/gamma_cut_acc),alpha)/alpha)
pile_up=gamma*gamma*np.exp(-np.power((gamma/gamma_eq),alpha)/alpha)
pile_up=pile_up/(gamma_eq*gamma_eq*np.exp(-np.power((gamma_eq/gamma_eq),alpha)/alpha))
g_ID=np.argmin(np.fabs(gamma-gamma_eq))
c=b[g_ID]*ratio_pile_up
a=pile_up*c
return a+b
The exponential cut-off factor in LPPL_pileup dist (
distr_func_lppl_pileup) is only applied to the LPPL part, which is usually fine in cases wheregamma_eqis large, but ifgamma_eq<1E7 andalphais also small (<0.5), this can produce a jump discontinuity atgamma_injin n(gamma), even worse for n(gamma)*gamma^3