diff --git a/raremodel-tf.py b/raremodel-tf.py index 9d75d5d..42bda47 100644 --- a/raremodel-tf.py +++ b/raremodel-tf.py @@ -78,83 +78,64 @@ return prefactor * _sum -class resonance(zfit.pdf.ZPDF): - _N_OBS = 1 # dimension, can be omitted - _PARAMS = ['mass', 'width', 'phase'] # the name of the parameters +def resonance(q, _mass, width, phase, scale): - def _unnormalized_pdf(self, x): + q2 = tf.pow(x, 2) - x = x.unstack_x() # returns a list with the columns: do x, y, z = ztf.unstack_x(x) for 3D - q2 = tf.pow(x, 2) + mmu = ztf.constant(pdg['muon_M']) - _mass = self.params['mass'] - width = self.params['width'] - phase = self.params['phase'] + p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2)) - mmu = ztf.constant(pdg['muon_M']) + p0 = 0.5 * ztf.sqrt(_mass**2 - 4*mmu**2) - p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2)) + gamma_j = tf.divide(p, q2) * _mass * width / p0 - p0 = 0.5 * ztf.sqrt(_mass**2 - 4*mmu**2) + #Calculate the resonance - gamma_j = tf.divide(p, q2) * _mass * width / p0 + _top = tf.complex(_mass * width, ztf.constant(0.0)) - #Calculate the resonance + _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j) - _top = tf.complex(_mass * width, ztf.constant(0.0)) + com = _top/_bottom - _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j) + #Rotate by the phase - com = _top/_bottom + r = tf.abs(com) - #Rotate by the phase + _phase = tf.angle(com) - r = tf.abs(com) + _phase += phase - _phase = tf.angle(com) + x = tf.cos(phase)*r + y = tf.sin(phase)*r - _phase += phase + com = tf.complex(x, y) - x = tf.cos(phase)*r - y = tf.sin(phase)*r + return com - com = tf.complex(x, y) +def bifur_gauss(q, mean, sigma_L, sigma_R, scale): - return com + x_left = tf.where(x < mean, x, False) -class bifur_gauss(zfit.pdf.ZPDF): - _N_OBS = 1 # dimension, can be omitted - _PARAMS = ['mean', 'amp', 'sigma_L', 'sigma_R'] # the name of the parameters + x_right = tf.where(q >= mean, x, False) - def _unnormalized_pdf(self, x): - x = ztf.unstack_x(x) # returns a list with the columns: do x, y, z = ztf.unstack_x(x) for 3D + #Calculate the exponential part of the cusp - mean = self.params['mean'] - amp = self.params['amp'] - sigma_L = self.params['sigma_L'] - sigma_R = self.params['sigma_R'] - - x_left = tf.where(x < mean, x, False) - - x_right = tf.where(q >= mean, x, False) - - #Calculate the exponential part of the cusp - - _exp_left = ztf.exp(- tf.pow((x_left-mean),2) / (2 * sigma_L**2)) - _exp_right = ztf.exp(- tf.pow((x_right-mean),2) / (2 * sigma_R**2)) + _exp_left = ztf.exp(- tf.pow((x_left-mean),2) / (2 * sigma_L**2)) + _exp_right = ztf.exp(- tf.pow((x_right-mean),2) / (2 * sigma_R**2)) - #Scale so the total area under curve is 1 and the top of the cusp is continuous + #Scale so the total area under curve is 1 and the top of the cusp is continuous - _exp = _exp_left + _exp_right + _exp = _exp_left + _exp_right - dgamma = scale*_exp/(ztf.sqrt(2*np.pi))*2*(sigma_L*sigma_R)/(sigma_L+sigma_R) + dgamma = scale*_exp/(ztf.sqrt(2*np.pi))*2*(sigma_L*sigma_R)/(sigma_L+sigma_R) - com = ztf.complex(dgamma, 0) + com = ztf.complex(dgamma, 0) - return com + return com -def axiv_nonres(q2, scale_axiv): +def axiv_nonres(q, scale_axiv): GF = ztf.constant(pdg['GF']) alpha_ew = ztf.constant(pdg['alpha_ew']) @@ -194,9 +175,9 @@ #Note sqrt(q2) comes from derivation as we use q2 and plot q - return ztf.to_complex(scale_axiv * prefactor1 * (bracket_left + bracket_middle) * 2 *ztf.sqrt(q2)) + return scale_axiv * prefactor1 * (bracket_left + bracket_middle) * 2 *ztf.sqrt(q2) -def vec(q2, scale_vec): +def vec(q, scale_vec, funcs): scale = self.params['scale'] @@ -207,7 +188,6 @@ Vtb = ztf.constant(pdg['Vtb']) Vts = ztf.constant(pdg['Vts']) C7eff = ztf.constant(pdg['C7eff']) - C9eff = ztf.constant(pdg['C9eff']) mmu = ztf.constant(pdg['muon_M']) mb = ztf.constant(pdg['bquark_M']) @@ -229,15 +209,26 @@ prefactor2 = kabs**2 * (1. - 1./3. * beta**2) - abs_bracket = tf.abs(C9eff * formfactor(q2, "+") + 2 * C7eff * (mb + ms)/(mB + mK) * formfactor(q2, "T"))**2 + abs_bracket = tf.abs(c9eff(x, funcs) * formfactor(q2, "+") + 2 * C7eff * (mb + ms)/(mB + mK) * formfactor(q2, "T"))**2 bracket_right = prefactor2 * abs_bracket #Note sqrt(q2) comes from derivation as we use q2 and plot q - return ztf.to_complex(scale_vec * prefactor1 * bracket_right * 2 * ztf.sqrt(q2)) + return scale_vec * prefactor1 * bracket_right * 2 * ztf.sqrt(q2) -class total_func(zfit.func.ZFunc): +def c9eff(q, funcs): + + C9eff_nr = ztf.constant(pdg['C9eff']) + + c9 = C9eff_nr + + for func in funcs: + c9 += func(q) + + return c9 + +class total_pdf(zfit.pdf.ZPDF): _N_OBS = 1 # dimension, can be omitted _PARAMS = ['jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width', 'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width', @@ -245,10 +236,26 @@ 'cusp_mass', 'sigma_L', 'sigma_R', 'cusp_scale' ] # the name of the parameters - def _func(self, x): - tot = tf.abs(vec_nonres + axiv_nonres) + def _unnormalized_pdf(self, x): - return a + def jpsi_res(q): + return resonace(q, _mass = self.params['jpsi_mass'], scale = self.params['jpsi_scale'], phase = self.params['jpsi_phase'], width = self.params['jpsi_width']) + + def psi2s_res(q): + return resonace(q, _mass = self.params['psi2s_mass'], scale = self.params['psi2s_scale'], phase = self.params['psi2s_phase'], width = self.params['psi2s_width']) + + def cusp(q): + return bifur_gauss(q, mean = self.params['cusp_mass'], sigma_L = self.params['sigma_L'], sigma_R = self.params['sigma_R'], scale = self.params['cusp_scale']) + + funcs = [jpsi_res(x), psi2s_res(x), cusp(x)] + + vec_f = vec(x, self.params['scale_vec'], funcs) + + axiv_nr = axiv(x, self.params['scale_axiv']) + + tot = vec_f + axiv_nr + + return tot #Load data @@ -268,18 +275,10 @@ #Build pdf -#Old false nonres +#vec/axiv scale -nr_of_pdfs = 2 - -scale_test = zfit.Parameter("scale_test", ztf.constant(1.0)) - -frac1 = zfit.Parameter("frac1", 1/nr_of_pdfs) -frac2 = zfit.Parameter("frac2", 1/nr_of_pdfs) - -axiv_nr = axiv_nonres(obs = obs) - -vec_nr = vec_nonres(obs = obs, scale = scale_test) +scale_v = zfit.Parameter("scale_vec", ztf.constant(1)) +scale_a = zfit.Parameter("scale_axiv", ztf.constant(1)) #jpsi @@ -288,21 +287,38 @@ jpsi_m = zfit.Parameter("jpsi_m", ztf.constant(jpsi_mass)) jpsi_w = zfit.Parameter("jpsi_w", ztf.constant(jpsi_width)) jpsi_p = zfit.Parameter("jpsi_p", ztf.constant(jpsi_phase)) - -# jpsi_res = resonance(obs = obs, mass = jpsi_m, width = jpsi_w, phase = jpsi_p) +jpsi_s = zfit.Parameter("jpsi_s", ztf.constant(jpsi_scale)) #psi2s psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"] +psi2s_m = zfit.Parameter("psi2s_m", ztf.constant(psi2s_mass)) +psi2s_w = zfit.Parameter("psi2s_w", ztf.constant(psi2s_width)) +psi2s_p = zfit.Parameter("psi2s_p", ztf.constant(psi2s_phase)) +psi2s_s = zfit.Parameter("psi2s_s", ztf.constant(psi2s_scale)) -total_pdf = total_func.as_pdf() +#cusp + +cusp_mass, sigma_R, sigma_L, cusp_scale = 3550, 3e-7, 200, 7 + +cusp_m = zfit.Parameter("cusp_m", ztf.constant(cusp_mass)) +sig_L = zfit.Parameter("sig_L", ztf.constant(sigma_L)) +sig_R = zfit.Parameter("sig_R", ztf.constant(sigma_R)) +cusp_s = zfit.Parameter("cusp_s", ztf.constant(cusp_scale)) + +total_f = total_pdf(obs=obs, jpsi_mass = jpsi_m, jpsi_scale = jpsi_s, jpsi_phase = jpsi_p, jpsi_width = jpsi_w, + psi2s_mass = psi2s_m, psi2s_scale = psi2s_s, psi2s_phase = psi2s_p, psi2s_width = psi2s_w, + scale_vec = scale_v, scale_axiv = scale_a, + cusp_mass = cusp_m, sigma_L = sig_L, sigma_R = sig_R, cusp_scale = cusp_s) + +# total_pdf = total_func.as_pdf() print(total_pdf.obs) -nll = zfit.loss.UnbinnedNLL(model=total_pdf, data=data) +nll = zfit.loss.UnbinnedNLL(model=total_pdf, data=data, fit_range = (x_min, x_max)) minimizer = zfit.minimize.MinuitMinimizer() -result = minimizer.minimize(nll, ) +result = minimizer.minimize(nll) param_errors = result.error()