Initial commit

.gitignore

@@ -99,3 +99,8 @@ ENV/
 
 # mypy
 .mypy_cache/
+
+#LaTeX
+*.aux
+*.log
+*.synctex.gz

TeX/fits.tex

@@ -0,0 +1,43 @@
+\documentclass[11pt, a4paper]{article}
+\usepackage{graphicx}
+\usepackage{amsmath, amsfonts, amssymb}
+\usepackage{float}
+\usepackage{pgf}
+\usepackage[margin=2cm]{geometry}
+
+\begin{document}
+\section{Imaginary Wilson coefficients and $R_K$ observables}
+\begin{figure}[H]
+\centering
+\includegraphics[width = 0.6\textwidth]{Figures/RK_im}
+\caption{Values of $R_K$ and $R_{K^*}$ whith imaginary Wilson coefficients.}
+\label{im:RK_im}
+\end{figure}
+
+ Figure \ref{im:RK_im} shows the values of the ratios $R_K$ and $R_{K^*}$ in their respective $q^2$ ranges, when both Wilson coefficients $C_9$ and $C_{10}$ are imaginary, and also we assume $C_9 = -C_{10}$. In all cases the minimum value for the ratio is attained at the SM point $C_9 = -C_{10} =0$. The addition of non-zero imaginary Wilson coefficients results in larger values of $R_K$ and $R_{K^*}$, at odds with the experimental values $R_{K^{(*)}}^{\mathrm{exp}} < R_{K^{(*)}}^{\mathrm{SM}}$.
+
+\begin{figure}[H]
+\centering
+\resizebox{0.6\textwidth}{!}{\input{Figures/fitim.pgf}}
+\caption{Allowed regions for imaginary Wilson coefficients.}
+\label{im:fitim}
+\end{figure}
+
+In Figure \ref{im:fitim} we plot the allowed regions for imaginary values of $C_9$ and $C_{10}$ when fitting to measurements of a series of $b\to s \mu^+\mu^-$ observables. In the fit we have included the ratios $R_K$ and $R_{K^*}$, the angular observables $P_4'$ and $P_5'$ and the branching ratios $\mathrm{BR}(B_s\to \mu^+ \mu^-)$ and $\mathrm{BR}(B^0 \to \mu^+ \mu^-)$. The fit was carried by the software \texttt{flavio}, wich computes the $\chi^2$ function with each $(C_9, C_{10})$ pair. The global minimum of the $\chi^2$ was found at $C_9 = 0.31\,i,\ C_{10} = 0.29\;i$. The pull of the SM, defined as the probability that the SM scenario can describe the best fit assuming that $\Delta \chi^2 = \chi^2_{\mathrm{SM}} - \chi^2_{\mathrm{min}}$ follows a $\chi^2$ distribution with 5 degrees of freedom, is of just $2.1\times 10^{-5} \sigma$. That is, the SM result and the best fit are indistinguishable from an statistical point of view. The shaded regions in the plot are $1\sigma$ (darker pink) and $2\sigma$ (lighter pink) away from the best fit.
+\section{Z' fit}
+\begin{figure}[H]
+\centering
+\resizebox{0.6\textwidth}{!}{\input{Figures/WCZ.pgf}}
+\caption{Bounds on $Z'$ parameter space imposed by $b\to s \mu^+ \mu^-$ decays and $B_s$ mixing.}\label{im:WCZ}
+\end{figure}
+
+Figure \ref{im:WCZ} shows the bounds on the $Z'$ mass $M_{Z'}$ and the imaginary coupling coupling $\lambda_{23}^Q$ (setting $\lambda_{22}^L=1$) imposed by $b\to s \mu^+ \mu^-$ decays and $B_s$ mixing. Blue lines correspond to the fit to $B_s$-mixing observables $\Delta M_s$ and $A_{\mathrm{CP}}^{\mathrm{mix}}$ (solid line: $\Delta \chi^2 = 1$, dashed lines $\Delta \chi^2 = 4$), red lines to $b\to s \mu^+ \mu^-$ as in Figure \ref{im:fitim} (solid line: $\Delta \chi^2 = 1$, dashed lines $\Delta \chi^2 = 4$), and green regions to the combined fit (dark green: $\Delta \chi^2 = 1$, medium green $\Delta \chi^2 = 4$, light green: $\Delta \chi^2 = 9$). For the $B_s$-mixing fit, the SM offers the best fit, while the best fits for the $b\to s \mu^+ \mu^-$ and global fits is ($M_{Z'} = 11.0\ \mathrm{TeV},\ \lambda_{23}^Q = 0.001\.i$), although the improvement with respect to the SM is negligible. In the low mass region ($M_{Z'} \lesssim 7\ \mathrm{TeV}$) the global fit is dominated by $b\to s \mu^+ \mu^-$ observables, and $B_s$-mixing becomes more important at larger masses. The allowed regions correspond to 'small' Wilson coefficients ($\mathrm{Im} C_{10} < 0.7$, $C_{bs}^{LL} > -3\times 10^{-4}$), except for a narrow region in the $B_s$ fit with $C_{bs}^{LL} \sim -2.7 \times 10^{-3}$.
+\section{Leptoquark fit}
+\begin{figure}[H]
+\centering
+\resizebox{0.6\textwidth}{!}{\input{Figures/WCLQ.pgf}}
+\caption{Bounds on $S_3$ leptoquark parameter space imposed by $b\to s \mu^+ \mu^-$ decays and $B_s$ mixing.}\label{im:WCLQ}
+\end{figure}
+
+Figure \ref{im:WCLQ} shows the bounds on the $S_3$ leptoquark mass $M_{S_3}$ and the imaginary coupling coupling $y^{QL}_{23} y^{QL*}_{23}$  imposed by $b\to s \mu^+ \mu^-$ decays and $B_s$ mixing. The observables used in the respective fits are the same as in Figure \ref{im:WCZ}. The SM scenario was prefered in all the fits. In contrast to $Z'$ models, $b\to s \mu^+ \mu^-$ observables remain the main contributors to the bounds at masses as large as $M_{S_3} \sim 50\ \mathrm{TeV}$.
+\end{document}

src/Wilsonffit.py

@@ -0,0 +1,50 @@
+import texfig # download from https://github.com/nilsleiffischer/texfig
+import flavio
+import flavio.plots
+import flavio.statistics.fits
+import matplotlib.pyplot as plt
+
+def wc(C9, C10):
+	'''
+	Wilson coeffients settings
+	'''
+	return {
+		'C9_bsmumu': C9*1j,
+		'C10_bsmumu': C10*1j,
+	}
+
+fast_fit = flavio.statistics.fits.FastFit(
+                name = "C9-C10 SMEFT fast fit",
+                observables = [ 'BR(Bs->mumu)', 'BR(B0->mumu)', ('<Rmue>(B0->K*ll)', 0.045, 1.1), ('<Rmue>(B0->K*ll)', 1.1, 6.0), ('<Rmue>(B+->Kll)', 1.0, 6.0), ('<P4p>(B0->K*mumu)', 1.1, 6.0) , ('<P5p>(B0->K*mumu)', 1.1, 6.0) ],
+                fit_wc_function = wc,
+                input_scale = 4.8,
+                include_measurements = ['LHCb Bs->mumu 2017', 'CMS Bs->mumu 2013', 'LHCb RK* 2017', 'LHCb B->Kee 2014', 'LHCb B->K*mumu 2015 P 1.1-6', 'ATLAS B->K*mumu 2017 P4p', 'ATLAS B->K*mumu 2017 P5p'],
+            )
+
+fast_fit.make_measurement(threads=2)
+
+
+def plot():
+	'''
+	Plots the allowed regions in the C9-C10 plane for imaginary Wilson coefficients
+	'''
+	fig = texfig.figure()
+	opt = dict(x_min=-2, x_max=2, y_min=-2, y_max=2, n_sigma=(1,2), interpolation_factor=5)
+	flavio.plots.likelihood_contour(fast_fit.log_likelihood, col=0, **opt, threads=2)
+	#flavio.plots.flavio_branding(y=0.07, x=0.05) #crashes LaTeX
+	plt.gca().set_aspect(1)
+	plt.axhline(0, c='k', lw=0.2)
+	plt.axvline(0, c='k', lw=0.2)
+	plt.plot(0.31, 0.29, marker='x') #compute best fit first!
+	plt.xlabel(r'$\mathrm{Im} C_9$')
+	plt.ylabel(r'$\mathrm{Im} C_{10}$')
+	texfig.savefig('fitim')
+
+
+def best_fit(x0=[0.3,0.3]):
+	'''
+	Computes the best fit starting at point x0 = [Im C9, Im C10]
+	'''
+	bf_global = fast_fit.best_fit(x0)
+	print('Global: C9='+str(bf_global['x'][0])+ 'i\tC10='+str(bf_global['x'][1])+'i\nchi2 = ' + str(bf_global['log_likelihood']))
+

src/measCVLL.yml

@@ -0,0 +1,7 @@
+Rule them all:
+ experiment: HFLAV
+ inspire: DiLuzio:2017fdq
+ url: http://arxiv.org/abs/arXiv:1712.06572
+ values:
+  Delta_Ms: 1.1688(14)e-11
+  ACP_mix: -0.021 ± 0.031

src/modelfit.py

@@ -0,0 +1,151 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+'''
+Fit of LFUV and Bs mixing observables to Wilson Coefficients in the Z' and LQ models
+'''
+
+from math import pi, sqrt, sin
+from cmath import phase
+import flavio
+import flavio.plots
+import flavio.statistics.fits
+import matplotlib.pyplot as plt
+from flavio.classes import Observable, Prediction, Parameter
+from flavio.physics.running.running import get_alpha
+import flavio.measurements
+
+
+Parameter('eta_LL')
+flavio.default_parameters.set_constraint('eta_LL', u'0.77 ± 0.02')
+Parameter('DM_s')
+flavio.default_parameters.set_constraint('DM_s', '1.1688 ± 0.111103 e-11')
+R = 1.3397e-3
+Parameter('beta_s')
+flavio.default_parameters.set_constraint('beta_s', u'0.01852 ± 0.00032')
+par = flavio.default_parameters.get_central_all()
+
+GF = par['GF']
+Vtb = sqrt(1-par['Vub']**2-par['Vcb']**2)
+Vts = sqrt(1-Vtb**2)
+sq2 = sqrt(2)
+
+def Delta_Ms(wc_obj, par):
+	CVLL = wc_obj.get_wc('bsbs', scale=4.8, par=par)['CVLL_bsbs']
+	return par['DM_s']*abs(1.0+CVLL/R)
+
+Observable('Delta_Ms')
+Prediction('Delta_Ms', lambda wc_obj, par: Delta_Ms(wc_obj,par))
+
+def ACP_mix(wc_obj, par):
+	CVLL = wc_obj.get_wc('bsbs', scale=4.8, par=par)['CVLL_bsbs']
+	phi = phase(1.0+CVLL/R)
+	return sin(phi-2*par['beta_s'])
+
+Observable('ACP_mix')
+Prediction('ACP_mix', ACP_mix)
+
+
+def wc_Z(lambdaQ, MZp):
+	"Wilson coefficients as functions of Z' couplings"
+	alpha = get_alpha(par, MZp+0.2)['alpha_e']
+	return {
+            'C9_bsmumu': -pi/(sq2* GF * (MZp+0.2)**2 *alpha) * lambdaQ/(Vtb*Vts) *1j,
+            'C10_bsmumu': pi/(sq2* GF * (MZp+0.2)**2 *alpha) * lambdaQ/(Vtb*Vts) *1j,
+            'CVLL_bsbs': par['eta_LL']/(4 * sq2 *GF * (MZp+0.2)**2) * (lambdaQ/(Vtb*Vts))**2 *(-1)
+	}
+
+def wc_LQ(yy, MS):
+	"Wilson coefficients as functions of Leptoquark couplings"
+	alpha = get_alpha(par, MS+0.2)['alpha_e']
+	return {
+		'C9_bsmumu': pi/(sq2* GF * (MS+0.2)**2 *alpha) * yy/(Vtb*Vts) * 1j,
+		'C10_bsmumu': -pi/(sq2* GF * (MS+0.2)**2 *alpha) * yy/(Vtb*Vts) * 1j,
+		'CVLL_bsbs': par['eta_LL']/(4 * sq2 *GF * (MS+0.2)**2) * (yy/(Vtb*Vts))**2 * 5 /(64*pi**2) * (-1)
+	}
+
+flavio.measurements.read_file('measCVLL.yml')
+wc = flavio.WilsonCoefficients()
+
+
+'''
+exp = []
+uncExp = []
+expD = {}
+uncExpD = {}
+for m in measurements:
+	expD.update(flavio.Measurement[m].get_central_all())
+	uncExpD.update(flavio.Measurement[m].get_1d_errors_random())
+
+unc = []
+uncTot = []
+observablesSM = [ ('BR(Bs->mumu)',), ('BR(B0->mumu)',), ('<Rmue>(B0->K*ll)', 0.045, 1.1), ('<Rmue>(B0->K*ll)', 1.1, 6.0), ('<Rmue>(B+->Kll)', 1.0, 6.0), ('<P4p>(B0->K*mumu)', 1.1, 6.0) , ('<P5p>(B0->K*mumu)', 1.1, 6.0), ('Delta_Ms',), ('ACP_mix',) ]
+measurements = ['LHCb Bs->mumu 2017', 'CMS Bs->mumu 2013', 'LHCb RK* 2017', 'LHCb B->Kee 2014', 'LHCb B->K*mumu 2015 P 1.1-6', 'ATLAS B->K*mumu 2017 P4p', 'ATLAS B->K*mumu 2017 P5p', 'Rule them all']
+for obs in observablesSM:
+	unc += [flavio.sm_uncertainty(*obs)]
+	if len(obs)==1:
+		obs = obs[0]
+	exp += [expD[obs]]
+	uncExp += [uncExpD[obs]]
+	uncTot += [sqrt(unc[-1]**2 + uncExpD[obs]**2) ]
+'''
+
+exp = [3.03448275862069e-09, 3.620689655172414e-10, 0.652542372881356, 0.6813559322033899, 0.745, 0.07, 0.01, 1.1688e-11, -0.021]
+uncTot = [1.0224162247111247e-09, 2.1633572944255006e-10, 0.09783105808975297, 0.1049727973615827, 0.09067364365960803, 0.32969017997333544, 0.23984462133464168, 1.1110302680969927e-12, 0.03128684613867722]
+
+
+observablesNP = [ ('BR(Bs->mumu)', wc), ('BR(B0->mumu)', wc), ('<Rmue>(B0->K*ll)', wc, 0.045, 1.1), ('<Rmue>(B0->K*ll)', wc, 1.1, 6.0), ('<Rmue>(B+->Kll)', wc, 1.0, 6.0), ('<P4p>(B0->K*mumu)', wc, 1.1, 6.0) , ('<P5p>(B0->K*mumu)', wc, 1.1, 6.0), ('Delta_Ms', wc), ('ACP_mix', wc) ]
+
+
+def chicalc(l , M, wc):
+	"Chi-squared statistic for the fit"
+	wc.set_initial(wc(l, M), scale=4.8)
+	chi2 = 0
+	for o in range(0, len(observablesNP)):
+		chi2 += ((flavio.np_prediction(*observablesNP[o]) - exp[o] ))**2/(uncTot[o]**2 )
+	return chi2
+
+def chicalcRK(l , M, wc):
+	"Chi-squared statistic for the fit to RK(*)-related observables"
+	wc.set_initial(wc(l, M), scale=4.8)
+	chi2 = 0
+	for o in range(0, len(observablesNP)-2):
+		chi2 += ((flavio.np_prediction(*observablesNP[o]) - exp[o] ))**2/(uncTot[o]**2 )
+	return chi2
+
+def chicalcBs(l , M, wc):
+	"Chi-squared statistic for the fit to Bs-mixing-related observables"
+	wc.set_initial(wc(l, M), scale=4.8)
+	chi2 = 0
+	for o in (-2,-1):
+		chi2 += ((flavio.np_prediction(*observablesNP[o]) - exp[o] ))**2/(uncTot[o]**2 )
+	return chi2
+
+
+def makefit(wc, stepM, stepL, maxM, maxL, filename):
+	'''
+	Fit and print results to file
+	wc: Function that computes Wilson Coefficients from the model parameters (e.g. wc_Z, wc_LQ)
+	stepM, maxM in TeV
+	'''
+	print('There we go!')
+	chiRK0 = chicalcRK(0, 5000, wc)
+	chiBs0 = chicalcBs(0, 5000, wc)
+	chitot0 = chiRK0 + chiBs0
+	numM = int(maxM/stepM)
+	numL = int(maxL/stepL)
+	for M in range(1, numM+1):
+		chi = []
+		for l in range(1, numL+1):
+			chiRK = chicalcRK(l*stepL, M*stepM*1000, wc)
+			chiBs = chicalcBs(l*stepL, M*stepM*1000, wc)
+			chitot = chiRK + chiBs
+			chi += [(chiRK, chiBs, chitot)]
+		f = open(filename, 'at')
+		f.write(str(M*stepM) + '\t0\t0\t0\t' + str(chiRK0) + '\t' + str(chiBs0) + '\t' + str(chitot0) + '\n' )
+		for l in range(1, numL+1):
+			WCs = wc(l*stepL, M*stepM*1000)
+			C9 = WCs['C9_bsmumu']
+			CVLL = WCs['CVLL_bsbs']
+			f.write(str(M*stepM) + '\t' + str(l*stepL) + '\t'  + str(C9.imag) + '\t' + str(CVLL) + '\t' + str(chi[l-1][0]) + '\t' + str(chi[l-1][1]) + '\t' + str(chi[l-1][2]) + '\n' )
+		f.close() 

src/plotter.py

@@ -0,0 +1,58 @@
+import texfig #download from https://github.com/nilsleiffischer/texfig
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+
+def readfile(filename, col, threshold):
+	'''
+	Reads the data in file and saves it in three lists (first column, second column, and column number col) if the z-value is under the threshold
+	'''
+	f = open(filename, 'rt')
+	xtot = []
+	ytot = []
+	ztot = []
+	for l in f.readlines():
+		znew = float(l.split('\t')[col])
+		if znew < threshold:
+			xtot += [float(l.split('\t')[0])]
+			ytot += [float(l.split('\t')[1])]
+			ztot += [znew]
+	f.close()
+	return [xtot, ytot, ztot]
+
+def triang(X):
+	'''
+	Interpolation by triangulation used to plot contour levels
+	'''
+	trg = tri.Triangulation(X[0], X[1])
+	refiner = tri.UniformTriRefiner(trg)
+	trg_ref, z_ref = refiner.refine_field(X[2], subdiv=3)
+
+def drawplot(filein, fileout):
+	'''
+	Draw the plot using shaded areas for the global fit and contour lines for Bs-only and RK-only fits
+	'''
+	fig = texfig.figure()
+
+	trgtot_ref, ztot_ref = triang(readfile(filein, -1, 12))
+	plt.tricontourf(trgtot_ref, ztot_ref, levels = [0.0, 1.0, 4.0, 9.0], colors = ('#008000', '#00FF00', '#BFFF80'))
+	trgtot = tri.Triangulation(xtot, ytot)
+	trgtot_Bs, ztot_Bs = triang(readfile(filein, -2, 12))
+	plt.tricontour(trgBs_ref, zBs_ref, levels = [1.0, 4.0], colors = 'b', linestyles = ('solid', 'dashed'))
+	trgRK_ref, zRK_ref = triang(readfile(filein, -3, 12))
+	plt.tricontour(trgRK_ref, zRK_ref, levels = [1.0, 4.0], colors = 'r', linestyles = ('solid', 'dashed'))
+
+	plt.xlabel(r"$M_{S_3} [\mathrm{TeV}]$")
+	plt.ylabel(r'$\mathrm{Im}\ y_{32}^{QL} y_{32}^{QL*}$')
+	#axes = plt.gca()
+	#axes.set_ylim([0, 1.5])
+	texfig.savefig(fileout)
+
+
+
+
+#Including the plot in LaTex doc:
+#% in the preamble
+#\usepackage{pgf}
+#% somewhere in your document
+#\input{WCLQ.pgf}
+