Source code for smefit.rge

# -*- coding: utf-8 -*-
import pathlib
import warnings
from copy import deepcopy
from functools import partial

import ckmutil.ckm
import jax.numpy as jnp
import numpy as np
import pandas as pd
import wilson
from numpy import ComplexWarning

from smefit import log
from smefit.loader import Loader
from smefit.wcxf import inverse_wcxf_translate, wcxf_translate

### Patch of a CKM function, so that the CP violating
### phase is set to gamma and not computed explicitly
### See https://github.com/wilson-eft/wilson/issues/113#issuecomment-2179273979
### This needs to be done before the import of wilson

# copying so we keep the original function
ckm_tree = deepcopy(ckmutil.ckm.ckm_tree)
ckmutil.ckm.ckm_tree = partial(ckm_tree, delta_expansion_order=0)
### End of patch

# switch off the SM - EFT mixing, since SMEFiT assumes that the
# RGE solution is linearised
# Keep a reference to the original beta function
original_beta = wilson.run.smeft.beta.beta


[docs] def beta_wrapper(C, HIGHSCALE=np.inf, *args, **kwargs): return original_beta(C, HIGHSCALE, *args, **kwargs)
wilson.run.smeft.beta.beta = beta_wrapper wilson.run.smeft.beta.beta_array = partial( wilson.run.smeft.beta.beta_array, HIGHSCALE=np.inf ) # Suppress the ComplexWarning warnings.filterwarnings("ignore", category=ComplexWarning) _logger = log.logging.getLogger(__name__) ########################### # Input parameter options # ########################### # copying so we could use the default parameters later default_params = wilson.run.smeft.smpar.p.copy() top_yukawa = { "Vus": 0.0, "Vub": 0.0, "Vcb": 0.0, "gamma": 0.0, "m_b": 0.0, "m_s": 0.0, "m_c": 0.0, "m_u": 0.0, "m_d": 0.0, "m_e": 0.0, "m_mu": 0.0, "m_tau": 0.0, } no_yukawa = { "Vus": 0.0, "Vub": 0.0, "Vcb": 0.0, "gamma": 0.0, "m_b": 0.0, "m_s": 0.0, "m_c": 0.0, "m_u": 0.0, "m_d": 0.0, "m_e": 0.0, "m_mu": 0.0, "m_tau": 0.0, "m_t": 0.0, } QCD_only = { "alpha_e": 0.0, "m_W": 1e-20, "m_h": 1e-20, } # gs at MZ alpha_s = 0.118 gs = np.sqrt(4 * np.pi * alpha_s)
[docs] def evolve_gs(scale): # evolve gs from MZ to scale beta0 = 11 - 2 / 3 * 6 return gs / np.sqrt( 1 + 2 * beta0 * gs**2 / (4 * np.pi) ** 2 * np.log(scale / 91.1876) )
[docs] class RGE: """ Class to compute the RGE matrix for the SMEFT Wilson coefficients. The RGE matrix is computed at the initial scale `init_scale` and evolved to the scale of interest. Parameters ---------- wc_names: list list of Wilson coefficient names to be included in the RGE matrix init_scale: float initial scale of the Wilson coefficients accuracy: str accuracy of the RGE integration. Options: "leadinglog" or "integrate". Default is 'integrate'. Inherited behaviour from wilson package. adm_QCD: bool if True, only the QCD anomalous dimension is used. Default is False. yukawa: str Yukawa parameterization to be used. Options: "top", "none" or "full". Default is "top". """ def __init__( self, wc_names, init_scale, accuracy="integrate", adm_QCD=False, yukawa="top", ): # order the Wilson coefficients alphabetically self.wc_names = sorted(wc_names) self.init_scale = init_scale self.accuracy = accuracy # set the anomalous dimension matrix parameters if yukawa == "top": wilson.run.smeft.smpar.p.update(**top_yukawa) elif yukawa == "none": wilson.run.smeft.smpar.p.update(**no_yukawa) elif yukawa == "full": wilson.run.smeft.smpar.p.update(**default_params) else: raise ValueError(f"Yukawa parameter not supported: {yukawa}") _logger.info(f"Using Yukawa parameterization: {yukawa}.") if adm_QCD: wilson.run.smeft.smpar.p.update(**QCD_only) _logger.info("Using anomalous dimension order: QCD.") else: _logger.info("Using anomalous dimension order: full.") _logger.info( f"Initializing RGE runner with initial scale {init_scale} GeV and accuracy {accuracy}." )
[docs] def RGEmatrix_dict(self, scale): """ Compute the RGE solution at the scale `scale` and return it as a dictionary. """ # compute the RGE matrix at the scale `scale` rge_matrix_dict = {} for wc_name, wc_vals in self.RGEbasis.items(): _logger.info(f"Computing RGE for {wc_name} at {scale} GeV.") wc_init = wilson.Wilson( wc_vals, scale=self.init_scale, eft="SMEFT", basis="Warsaw" ) wc_init.set_option("smeft_accuracy", self.accuracy) wc_final = wc_init.match_run(scale=scale, eft="SMEFT", basis="Warsaw") # Remove small values wc_final_vals = { key: value for key, value in wc_final.dict.items() if abs(value) > 1e-10 } # check that imaginary values are small if any(abs(val.imag) > 1e-10 for val in wc_final_vals.values()): raise ValueError( f"Imaginary values in Wilson coefficient for operator {wc_name}." ) rge_matrix_dict[wc_name] = self.map_to_smefit(wc_final_vals, scale) return rge_matrix_dict
[docs] def RGEmatrix(self, scale): """ Compute the RGE solution at the scale `scale` and return it as a pandas DataFrame. """ # compute the RGE matrix dict at the scale `scale` rge_matrix_dict = self.RGEmatrix_dict(scale) # create the RGE matrix as pandas dataframe rge_matrix = pd.DataFrame( columns=self.wc_names, index=self.all_ops, dtype=float ) for wc_name, wc_dict in rge_matrix_dict.items(): for op in self.all_ops: rge_matrix.loc[op, wc_name] = wc_dict.get(op, 0.0) # if there are rows with all zeros, remove them rge_matrix = rge_matrix.loc[(rge_matrix != 0).any(axis=1)] return rge_matrix
@property def RGEbasis(self): """ Returns the RGE basis translated from smefit to Warsaw. """ # computes the translation from the smefit basis to the Warsaw basis # as expected by the Wilson package wc_basis = {} for wc_name in self.wc_names: try: wcxf_dict = wcxf_translate[wc_name] except KeyError: _logger.warning( f"Wilson coefficient {wc_name} not present in the WCxf translation dictionary." ) _logger.warning( "Assuming it is a UV coupling and associating it to the null vector." ) wc_basis[wc_name] = {} continue wc_warsaw_name = wcxf_dict["wc"] if "value" not in wcxf_dict: wc_warsaw_value = [1] * len(wcxf_dict["wc"]) else: # check if value is gs # (this is a special case for OtG) if wcxf_dict["value"] == ["gs"]: wc_warsaw_value = [evolve_gs(self.init_scale)] else: wc_warsaw_value = wcxf_dict["value"] # 1e-6 is because the Warsaw basis is in GeV^-2 wc_value = { wc: val * 1e-6 for wc, val in zip(wc_warsaw_name, wc_warsaw_value) } wc_basis[wc_name] = wc_value return wc_basis
[docs] def map_to_smefit(self, wc_final_vals, scale): """ Map the Wilson coefficients from the Warsaw basis to the SMEFiT basis. """ wc_dict = {} for wc_basis, wc_inv_dict in inverse_wcxf_translate.items(): wc_warsaw_name = wc_inv_dict["wc"] if "coeff" not in wc_inv_dict: wc_warsaw_coeff = [1] * len(wc_warsaw_name) else: # check if coeff is 1/gs # (this is a special case for OtG) if wc_inv_dict["coeff"] == ["1/gs"]: wc_warsaw_coeff = [1 / evolve_gs(scale)] else: wc_warsaw_coeff = wc_inv_dict["coeff"] value = 0.0 for wc, coeff in zip(wc_warsaw_name, wc_warsaw_coeff): if wc in wc_final_vals: # 1e6 is to transform from GeV^-2 to TeV^2 value += 1e6 * wc_final_vals[wc].real * coeff wc_dict[wc_basis] = value return wc_dict
@property def all_ops(self): return sorted(wcxf_translate.keys())
[docs] def RGEevolve(self, wcs, scale): """ Evolve the Wilson coefficients from the initial scale to the scale of interest. """ wc_wilson = {} for op, values in self.RGEbasis.items(): for key in values: if key not in wc_wilson: wc_wilson[key] = values[key] * wcs[op] else: wc_wilson[key] += values[key] * wcs[op] wc_init = wilson.Wilson( wc_wilson, scale=self.init_scale, eft="SMEFT", basis="Warsaw" ) wc_init.set_option("smeft_accuracy", self.accuracy) wc_final = wc_init.match_run(scale=scale, eft="SMEFT", basis="Warsaw").dict # remove small values wc_final = {key: value for key, value in wc_final.items() if abs(value) > 1e-10} return self.map_to_smefit(wc_final, scale)
[docs] def load_scales(datasets, theory_path, default_scale=1e3): """ Load the energy scales for the datasets. Parameters ---------- datasets: list list of datasets theory_path: str path to the theory files default_scale: float default scale to use if the dataset does not have a scale Returns ------- scales: list list of scales for the datasets """ scales = [] for dataset in np.unique(datasets): Loader.theory_path = pathlib.Path(theory_path) # dummy call just to get the scales _, _, _, _, dataset_scales = Loader.load_theory( dataset, operators_to_keep={}, order="LO", use_quad=False, use_theory_covmat=False, use_multiplicative_prescription=False, ) # check that dataset_scales is not a list filled with None # otherwise, assume the initial scale if not all(scale is None for scale in dataset_scales): scales.extend(dataset_scales) else: scales.extend([default_scale] * len(dataset_scales)) _logger.info(f"Loaded scales for dataset {dataset}: {dataset_scales}") return scales
[docs] def load_rge_matrix(rge_dict, coeff_list, datasets=None, theory_path=None): """ Load the RGE matrix for the SMEFT Wilson coefficients. Parameters ---------- rge_dict: dict dictionary with the RGE input parameter options coeff_list: list list of Wilson coefficients to be included in the RGE matrix datasets: list list of datasets theory_path: str path to the theory files Returns ------- rgemat: numpy.ndarray RGE matrix """ init_scale = rge_dict.get("init_scale", 1e3) obs_scale = rge_dict.get("obs_scale", 91.1876) smeft_accuracy = rge_dict.get("smeft_accuracy", "integrate") adm_QCD = rge_dict.get("adm_QCD", "full") yukawa = rge_dict.get("yukawa", "top") rge_runner = RGE(coeff_list, init_scale, smeft_accuracy, adm_QCD, yukawa) # if it is a float, it is a static scale if type(obs_scale) is float or type(obs_scale) is int: rgemat = rge_runner.RGEmatrix(obs_scale) gen_operators = list(rgemat.index) operators_to_keep = {k: {} for k in gen_operators} return rgemat.values, operators_to_keep elif obs_scale == "dynamic": scales = load_scales(datasets, theory_path, default_scale=init_scale) operators_to_keep = {} rgemat = [] rge_cache = {} for scale in scales: if scale == init_scale: # produce an identity matrix with row and columns coeff_list rgemat_scale = pd.DataFrame( np.eye(len(coeff_list), len(coeff_list)), columns=sorted(coeff_list), index=sorted(coeff_list), ) # Check if the RGE matrix has already been computed elif scale in rge_cache: rgemat_scale = rge_cache[scale] else: rgemat_scale = rge_runner.RGEmatrix(scale) gen_operators = list(rgemat_scale.index) # Fill with the operators if not already present in the dictionary for op in gen_operators: if op not in operators_to_keep: operators_to_keep[op] = {} rge_cache[scale] = rgemat_scale rgemat.append(rgemat_scale) # now loop through the rgemat and if there are operators # that are not present in the matrix, fill them with zeros for mat in rgemat: for op in operators_to_keep: if op not in mat.index: mat.loc[op] = np.zeros(len(mat.columns)) # order the rows alphabetically in the index mat.sort_index(inplace=True) # now stack the matrices in a 3D array stacked_mats = jnp.stack([mat.values for mat in rgemat]) return stacked_mats, operators_to_keep else: raise ValueError( f"obs_scale must be either a float/int or 'dynamic'. Passed: {obs_scale}" )