API

The user interface is given in terms of the ParamPoint class. For each parameter point that one wishes to analyze, an instance of this class has to be created, giving the values of the free parameters as input.

ParamPoint class

The ParamPoint class contains the functions to perform the possible computations, such as the RGE running, the calculation of the branching ratios, or to confront the parameter point with the various theoretical and experimental constraints. The functions are documented in detail below.

s2hdmTools.paramPoint.ParamPoint

The base class for an S2HDM parameter point.

From the input parameters the required theoretical predictions are computed. Various functions exist in order to confront a parameter point with theoretical or experimental constraints.

Source code in s2hdmTools/paramPoint.py

class ParamPoint:
    """
    The base class for an S2HDM parameter point.

    From the input parameters the required
    theoretical predictions are computed.
    Various functions exist in order to confront
    a parameter point with theoretical or
    experimental constraints.
    """

    def __init__(self, dc, debug=True):
        """
        Initialize an instance of a parameter
        point given a set of input parameters.
        The parameters have to be given as
        the dictionary dc. There are two
        possible input formats:

        Angle input, for instance:
        ```
        dc = {
            'type': 2,
            'tb': 1.26,
            'al1': 1.27,
            'al2': -1.08,
            'al3': -1.24,
            'mH1': 96.5,
            'mH2': 125.09,
            'mH3': 535.86,
            'mXi': 266.0,
            'mA': 712.578,
            'mHp': 737.8,
            'vs': 272.0,
            'm12sq': 80644.0}
        ```

        Lambda input, for instance:
        ```
        dc = {
            'type': 2,
            'tb': 0.91655,
            'mXi': 134.03,
            'vs': 64.987,
            'm12sq': 1.7696e5,
            'lam1': 1.5297,
            'lam2': 1.2074,
            'lam3': 1.5741,
            'lam4': 5.3967,
            'lam5': -7.8556,
            'lam6': 6.0689,
            'lam7': 0.80378,
            'lam8': -0.83745}
        ```

        Args:
            dc (dict): Dictionary with input values
                for free parameters, either in angle-
                or in lambda-input format.
            debug (bool): If set to `True` then the
                debugging functions are called.
        """

        self.dc = dc
        self.read_constants()
        self.read_input()
        if debug:
            self.check_masses()
        self.debug = debug

    def read_constants(self):

        self.v = constants.vSM
        self.mt = constants.mt
        self.mb = constants.mb
        self.ml = constants.ml
        self.mW = constants.mW
        self.mZ = constants.mZ
        self.alphas = constants.als

        self.g1 = 2. * np.sqrt(self.mZ**2 - self.mW**2) / self.v
        self.g2 = 2. * self.mW / self.v
        self.g3 = np.sqrt(4. * np.pi * self.alphas)

    def read_input(self):

        # Find input format: OS or DR
        if 'al1' in self.dc:
            self.read_angle_input()
            self.get_phi_vevs()
            self.get_lambdas()
            self.get_bilinears()
        elif 'lam1' in self.dc:
            self.read_lambda_input()
            self.get_phi_vevs()
            self.get_angle_input()
        else:
            raise RuntimeError(
                'Invalid input format.')
        if not self.yuktype in [1, 2, 3, 4]:
            raise ValueError(
                'Wrong input for Yukawa type.')

    def read_angle_input(self):

        dc = self.dc

        # Read values
        self.yuktype = dc['type']
        self.tb = dc['tb']
        self.al1 = dc['al1']
        self.al2 = dc['al2']
        self.al3 = dc['al3']
        self.mH1 = dc['mH1']
        self.mH2 = dc['mH2']
        self.mH3 = dc['mH3']
        self.mXi = dc['mXi']
        self.mA = dc['mA']
        self.mHp = dc['mHp']
        self.vs = dc['vs']
        self.m12sq = dc['m12sq']

        if not (self.mH1 <= self.mH2 <= self.mH3):
            raise ValueError(
                'Only normal mass ordering supported. ' +
                'Change input accordingly.')

        self.check_angle_ranges()

        self.calc_rot_matrix()

    def check_angle_ranges(self):

        a1 = self.al1
        a2 = self.al2
        a3 = self.al3

        pi2 = np.pi / 2.0

        # for a in [a1, a2, a3]:

        #     if not (-pi2 <= a <= pi2):

        #         raise ValueError(
        #             'Choose angle within the range ' +
        #             '-pi / 2 <= alpha_i <= pi / 2.')

    def read_lambda_input(self):

        dc = self.dc

        # Read values
        self.yuktype = dc['type']
        self.tb = dc['tb']
        self.vs = dc['vs']
        self.mXi = dc['mXi']
        self.m12sq = dc['m12sq']
        self.lam1 = dc['lam1']
        self.lam2 = dc['lam2']
        self.lam3 = dc['lam3']
        self.lam4 = dc['lam4']
        self.lam5 = dc['lam5']
        self.lam6 = dc['lam6']
        self.lam7 = dc['lam7']
        self.lam8 = dc['lam8']

        if self.mXi < 0.:
            raise ValueError(
                'mXi has to be positive.')

        if self.m12sq < 0.:
            raise ValueError(
                'm12sq has to be positive.')

    def get_phi_vevs(self):

        beta = np.arctan(self.tb)
        self.cb = np.cos(beta)
        self.sb = np.sin(beta)
        self.v1 = self.v * self.cb
        self.v2 = self.v * self.sb

    def calc_rot_matrix(self):

        c = np.cos
        s = np.sin

        a1 = self.al1
        a2 = self.al2
        a3 = self.al3

        R = np.zeros(shape=(3, 3))

        R[0, 0] = c(a1) * c(a2)
        R[0, 1] = c(a2) * s(a1)
        R[0, 2] = s(a2)
        R[1, 0] = -(c(a3) * s(a1)) - c(a1) * s(a2) * s(a3)
        R[1, 1] = c(a1) * c(a3) - s(a1) * s(a2) * s(a3)
        R[1, 2] = c(a2) * s(a3)
        R[2, 0] = -(c(a1) * c(a3) * s(a2)) + s(a1) * s(a3)
        R[2, 1] = -(c(a3) * s(a1) * s(a2)) - c(a1) * s(a3)
        R[2, 2] = c(a2) * c(a3)

        self.R = R

    def get_lambdas(self):

        mH1 = self.mH1
        mH2 = self.mH2
        mH3 = self.mH3
        m12sq = self.m12sq
        tb = self.tb
        R = self.R
        v = self.v
        mHp = self.mHp
        mA = self.mA
        vs = self.vs

        self.lam1 = (
            (1. + tb**2) * (
                -(m12sq * tb) +
                mH1**2 * R[0, 0]**2 +
                mH2**2 * R[1, 0]**2 +
                mH3**2 * R[2, 0]**2)) / v**2

        self.lam2 = ((1. + tb**2) * (
            -m12sq +
            mH1**2 * tb * R[0, 1]**2 +
            mH2**2 * tb * R[1, 1]**2 +
            mH3**2 * tb * R[2, 1]**2)) / (tb**3 * v**2)

        self.lam3 = (
            2. * mHp**2 * tb -
            m12sq * (1. + tb**2) +
            (1. + tb**2) * (
                mH1**2 * R[0, 0] * R[0, 1] +
                mH2**2 * R[1, 0] * R[1, 1] +
                mH3**2 * R[2, 0] * R[2, 1])) / (tb * v**2)

        self.lam4 = (
            m12sq + mA**2 * tb -
            2. * mHp**2 * tb +
            m12sq * tb**2) / (tb * v**2)

        self.lam5 = (
            m12sq - mA**2 * tb + m12sq * tb**2) / (tb * v**2)

        self.lam6 = (
            mH1**2 * R[0, 2]**2 +
            mH2**2 * R[1, 2]**2 +
            mH3**2 * R[2, 2]**2) / vs**2

        self.lam7 = (np.sqrt(1. + tb**2) * (
            mH1**2 * R[0, 0] * R[0, 2] +
            mH2**2 * R[1, 0] * R[1, 2] +
            mH3**2 * R[2, 0] * R[2, 2])) / (v * vs)

        self.lam8 = (np.sqrt(1. + tb**2) * (
            mH1**2 * R[0, 1] * R[0, 2] +
            mH2**2 * R[1, 1] * R[1, 2] +
            mH3**2 * R[2, 1] * R[2, 2])) / (tb * v * vs)

    def get_bilinears(self):

        L1 = self.lam1
        L2 = self.lam2
        L3 = self.lam3
        L4 = self.lam4
        L5 = self.lam5
        L6 = self.lam6
        L7 = self.lam7
        L8 = self.lam8
        v1 = self.v1
        v2 = self.v2
        vs = self.vs
        m12sq = self.m12sq

        self.m11sq = -(L1 * v1**2) / 2. + (m12sq * v2) / v1 - \
            (L3 * v2**2) / 2. - (L4 * v2**2) / 2. - \
                (L5 * v2**2) / 2. - (L7 * vs**2) / 2.

        self.m22sq = -(L3 * v1**2) / 2. - (L4 * v1**2) / 2. - \
            (L5 * v1**2) / 2. + (m12sq * v1) / v2 - \
                (L2 * v2**2) / 2. - (L8 * vs**2) / 2.

        self.mXsq = self.mXi**2

        self.mSsq = self.mXsq - L7 * v1**2 - L8 * v2**2 - L6 * vs**2

    def get_angle_input(self):

        self.get_bilinears()

        L1 = self.lam1
        L2 = self.lam2
        L3 = self.lam3
        L4 = self.lam4
        L5 = self.lam5
        L6 = self.lam6
        L7 = self.lam7
        L8 = self.lam8
        tb = self.tb
        sb = self.sb
        cb = self.cb
        v = self.v
        vs = self.vs
        m12sq = self.m12sq

        # CP evens
        M = np.zeros(shape=(3, 3))
        M[0, 0] = L1 * v**2 * cb**2 + m12sq * tb
        M[0, 1] = -m12sq + L3 * v**2 * cb * sb + \
            L4 * v**2 * cb * sb + L5 * v**2 * cb * sb
        M[1, 0] = M[0, 1]
        M[0, 2] = L7 * v * vs * cb
        M[2, 0] = M[0, 2]
        M[1, 1] = m12sq / tb + L2 * v**2 * sb**2
        M[1, 2] = L8 * v * vs * sb
        M[2, 1] = M[1, 2]
        M[2, 2] = L6 * vs**2
        diag, R = np.linalg.eig(M)

        if np.any(diag < 0.):
            raise ValueError(
                'Tachyonic scalar state.')

        # Get mass ordering mh1 < mh2 < mh3
        switched = True

        while switched:

            switched = False

            for i in range(0, 2):

                if diag[i] > diag[i + 1]:

                    temp = diag[i]
                    diag[i] = diag[i + 1]
                    diag[i + 1] = temp

                    R[:, [i, i + 1]] = R[:, [i + 1, i]]

                    switched = True

        self.mH1 = np.sqrt(diag[0])
        self.mH2 = np.sqrt(diag[1])
        self.mH3 = np.sqrt(diag[2])

        R = R.T
        self.R = R

        self.al2 = np.arcsin(R[0, 2])
        self.al1 = np.arcsin(R[0, 1] / np.cos(self.al2))
        self.al3 = np.arcsin(R[1, 2] / np.cos(self.al2))

        self.determine_signs_alphas()

        # CP odds
        M = np.zeros(shape=(2, 2))
        M[0, 0] = -L5 * v**2 * sb**2 + m12sq * tb
        M[0, 1] = -m12sq + L5 * v**2 * cb * sb
        M[1, 0] = M[0, 1]
        M[1, 1] = -L5 * v**2 * cb**2 + m12sq / tb
        Rp = np.array([[cb, sb], [-sb, cb]])
        diag = np.dot(np.dot(Rp, M), Rp.T)
        if diag[1, 1] > 0.:
            self.mA = np.sqrt(diag[1, 1])
        else:
            raise ValueError(
                'Tachyonic pseudoscalar state.')

        # Chargeds
        M = np.zeros(shape=(2, 2))
        M[0, 0] = -L4 * v**2 * sb**2 / 2 - \
            L5 * v**2 * sb**2 / 2 + m12sq * tb
        M[0, 1] = -m12sq + L4 * v**2 * sb * cb / 2 + \
            L5 * v**2 * sb * cb / 2
        M[1, 0] = M[0, 1]
        M[1, 1] = -L4 * v**2 * cb**2 / 2 - \
            L5 * v**2 * cb**2 / 2 + m12sq / tb
        Rp = np.array([[cb, sb], [-sb, cb]])
        diag = np.dot(np.dot(Rp, M), Rp.T)
        if diag[1, 1] > 0.:
            self.mHp = np.sqrt(diag[1, 1])
        else:
            raise ValueError(
                'Tachyonic charged scalar state.')

    def determine_signs_alphas(self):

        R = self.R

        Rcheck = np.zeros(shape=(3, 3))

        c = np.cos
        s = np.sin

        L1 = self.lam1
        L2 = self.lam2
        L3 = self.lam3
        L4 = self.lam4
        L5 = self.lam5
        L6 = self.lam6
        L7 = self.lam7
        L8 = self.lam8
        tb = self.tb
        sb = self.sb
        cb = self.cb
        v = self.v
        vs = self.vs
        m12sq = self.m12sq

        M = np.zeros(shape=(3, 3))
        M[0, 0] = L1 * v**2 * cb**2 + m12sq * tb
        M[0, 1] = -m12sq + L3 * v**2 * cb * sb + \
            L4 * v**2 * cb * sb + L5 * v**2 * cb * sb
        M[1, 0] = M[0, 1]
        M[0, 2] = L7 * v * vs * cb
        M[2, 0] = M[0, 2]
        M[1, 1] = m12sq / tb + L2 * v**2 * sb**2
        M[1, 2] = L8 * v * vs * sb
        M[2, 1] = M[1, 2]
        M[2, 2] = L6 * vs**2

        signpos = [
            [1, 1, 1],
            [-1, 1, 1],
            [1, -1, 1],
            [1, 1, -1],
            [-1, -1, 1],
            [-1, 1, -1],
            [1, -1, -1],
            [-1, -1, -1]]

        i = 0
        diag = np.array([
            [self.mH1**2, 0., 0.],
            [0., self.mH2**2, 0.],
            [0., 0., self.mH3**2]])

        diagcheck = np.zeros(shape=(3, 3))

        try:

            while not np.allclose(diag, diagcheck, atol=1e-6):

                a1 = signpos[i][0] * self.al1
                a2 = signpos[i][1] * self.al2
                a3 = signpos[i][2] * self.al3

                Rcheck[0, 0] = c(a1) * c(a2)
                Rcheck[0, 1] = c(a2) * s(a1)
                Rcheck[0, 2] = s(a2)
                Rcheck[1, 0] = -(c(a3) * s(a1)) - c(a1) * s(a2) * s(a3)
                Rcheck[1, 1] = c(a1) * c(a3) - s(a1) * s(a2) * s(a3)
                Rcheck[1, 2] = c(a2) * s(a3)
                Rcheck[2, 0] = -(c(a1) * c(a3) * s(a2)) + s(a1) * s(a3)
                Rcheck[2, 1] = -(c(a3) * s(a1) * s(a2)) - c(a1) * s(a3)
                Rcheck[2, 2] = c(a2) * c(a3)

                diagcheck = np.dot(np.dot(Rcheck, M), Rcheck.T)

                i += 1

        except IndexError:

            raise RuntimeError(
                'Could not parametrize mixing matrix with al1, al2 and al3.')

        self.al1 = a1
        self.al2 = a2
        self.al3 = a3

        self.R = Rcheck

    def check_masses(self):

        L1 = self.lam1
        L2 = self.lam2
        L3 = self.lam3
        L4 = self.lam4
        L5 = self.lam5
        L6 = self.lam6
        L7 = self.lam7
        L8 = self.lam8

        tb = self.tb
        sb = self.sb
        cb = self.cb

        v = self.v
        vs = self.vs
        m12sq = self.m12sq

        R = self.R

        # CP evens
        M = np.zeros(shape=(3, 3))
        M[0, 0] = L1 * v**2 * cb**2 + m12sq * tb
        M[0, 1] = -m12sq + L3 * v**2 * cb * sb + \
            L4 * v**2 * cb * sb + L5 * v**2 * cb * sb
        M[1, 0] = M[0, 1]
        M[0, 2] = L7 * v * vs * cb
        M[2, 0] = M[0, 2]
        M[1, 1] = m12sq / tb + L2 * v**2 * sb**2
        M[1, 2] = L8 * v * vs * sb
        M[2, 1] = M[1, 2]
        M[2, 2] = L6 * vs**2
        diag = np.dot(np.dot(R, M), R.T)
        M1 = np.sqrt(diag[0, 0])
        M2 = np.sqrt(diag[1, 1])
        M3 = np.sqrt(diag[2, 2])
        if not np.allclose(
            np.array([M1, M2, M3]),
            np.array([self.mH1, self.mH2, self.mH3])):
            raise RuntimeError(
                'Problem with input: Scalar masses do not fit.')

        # CP odds
        M = np.zeros(shape=(2, 2))
        M[0, 0] = -L5 * v**2 * sb**2 + m12sq * tb
        M[0, 1] = -m12sq + L5 * v**2 * cb * sb
        M[1, 0] = M[0, 1]
        M[1, 1] = -L5 * v**2 * cb**2 + m12sq / tb
        Rp = np.array([[cb, sb], [-sb, cb]])
        diag = np.dot(np.dot(Rp, M), Rp.T)
        M1 = np.sqrt(abs(diag[0, 0]))
        M2 = np.sqrt(diag[1, 1])
        if not np.allclose(
            np.array([M1, M2]),
            np.array([0.0, self.mA]),
            atol=1e-3): # Use atol here because comparing to zero
            raise RuntimeError(
                'Problem with input: Pseudoscalar masses do not fit.')

        # Chargeds
        M = np.zeros(shape=(2, 2))
        M[0, 0] = -L4 * v**2 * sb**2 / 2 - \
            L5 * v**2 * sb**2 / 2 + m12sq * tb
        M[0, 1] = -m12sq + L4 * v**2 * sb * cb / 2 + \
            L5 * v**2 * sb * cb / 2
        M[1, 0] = M[0, 1]
        M[1, 1] = -L4 * v**2 * cb**2 / 2 - \
            L5 * v**2 * cb**2 / 2 + m12sq / tb
        Rp = np.array([[cb, sb], [-sb, cb]])
        diag = np.dot(np.dot(Rp, M), Rp.T)
        M1 = np.sqrt(abs(diag[0, 0]))
        M2 = np.sqrt(diag[1, 1])
        if not np.allclose(
            np.array([M1, M2]),
            np.array([0.0, self.mHp]),
            atol=1e-3): # Use atol here because comparing to zero
            raise RuntimeError(
                'Problem with input: Charged scalar masses do not fit.')

    def calculate_branching_ratios(self):
        """
        Calculate the branching ratios of
        the Higgs bosons.

        Branching ratios are stored in:

        ```
        self.b_H1
        self.b_H2
        self.b_H3
        self.b_A
        self.b_Hp
        ```

        Total widths are stored in:

        ```
        self.wTot
        ```
        """

        n2hdecay = CallN2HDecay(self)

        invs = calc_gammas_inv(self)

        self.b_A = n2hdecay.b_A
        self.b_Hp = n2hdecay.b_Hp
        self.b_t = n2hdecay.b_t

        gamH1 = {}
        wH1 = n2hdecay.w['H1']
        for key in n2hdecay.b_H1:
            gamH1[key] = n2hdecay.b_H1[key] * wH1
        gamH1['XX'] = invs[0]
        wH1 += invs[0]
        self.b_H1 = {}
        for key in gamH1:
            self.b_H1[key] = gamH1[key] / wH1
        if self.debug:
            if abs(sum(self.b_H1.values()) - 1.) > 1e-6:
                raise RuntimeError(
                    'Branching ratios of H1 do not add up to one.')

        gamH2 = {}
        wH2 = n2hdecay.w['H2']
        for key in n2hdecay.b_H2:
            gamH2[key] = n2hdecay.b_H2[key] * wH2
        gamH2['XX'] = invs[1]
        wH2 += invs[1]
        self.b_H2 = {}
        for key in gamH2:
            self.b_H2[key] = gamH2[key] / wH2
        if self.debug:
            if abs(sum(self.b_H2.values()) - 1.) > 1e-6:
                raise RuntimeError(
                    'Branching ratios of H2 do not add up to one.')

        gamH3 = {}
        wH3 = n2hdecay.w['H3']
        for key in n2hdecay.b_H3:
            gamH3[key] = n2hdecay.b_H3[key] * wH3
        gamH3['XX'] = invs[2]
        wH3 += invs[2]
        self.b_H3 = {}
        for key in gamH3:
            self.b_H3[key] = gamH3[key] / wH3
        if self.debug:
            if abs(sum(self.b_H3.values()) - 1.) > 1e-6:
                raise RuntimeError(
                    'Branching ratios of H3 do not add up to one.')

        self.wTot = {
            'A': n2hdecay.w['A'],
            'H1': wH1,
            'H2': wH2,
            'H3': wH3,
            'Hp': n2hdecay.w['Hp'],
            't': n2hdecay.w['t']}

    def check_tree_pert_uni(
            self, scale=None, num=100000,
            loop_level=2,
            cutoff=EightPi):
        """
        Check against tree-level perturbative
        unitarity constraints.

        Args:
            scale (float): Energy scale up to
                which the constraints should be
                applied. If not given, the check
                is only performed at the electroweak scale
                $(\mu = v = 246\ \mathrm{GeV})$.
                Value must be larger than $v$.
            num (int): Number of points that are calculated
                for the RGE running when `scale` is given.
            loop_level (int): Loop-level of RGE running.
                Can be set to 1 or 2.
            cutoff (float): Defines
                the upper limit for the eigenvalues
                of the scattering matrix.

        Returns:
            bool/dict: If `scale=None` then return `True`
                if constraints are fulfilled and `False`
                otherwise. If `scale` is given then return
                dictionary with information about validity
                depending on the energy scale.
        """

        trPrU.CUTOFF = cutoff

        if scale is None:

            raw = trPrU.pertuni(
                self.lam1, self.lam2, self.lam3, self.lam4,
                self.lam5, self.lam6, self.lam7, self.lam8)

            self.tree_pert_uni_vals = raw[0]

            return raw[1]

        else:

            if scale < self.v:

                raise ValueError(
                    'scale has to be larger than initial ' +
                    'scale (=vSM).')

            dcrge, solrge = self.run_to_scale(
                scale, num=num,
                loop_level=loop_level)

            lamsAtScale = solrge[:, 5:13]

            dcL = self._check_lams_for_landau(lamsAtScale, dcrge, num)
            try:
                upper_ind = dcL['Landau_index']
            except KeyError:
                upper_ind = num

            dcP = trPrU.pertuni_to_scale(
                self,
                lamsAtScale[0:upper_ind, ...],
                dcrge['scale'][0:upper_ind])

            return {**dcL, **dcP}

    def _check_lams_for_landau(self, lamsAtScale, dcrge, num):

        for i in range(0, num):

            if (np.abs(lamsAtScale[i, ...]) > 1.0e3).any():

                cutoffL = dcrge['scale'][i]
                IsLandau = True
                iL = i
                break

        else:

                IsLandau = False

        dc = {
            'IsLandau': IsLandau,
            'num': num,
            'Scale_max': dcrge['scale'][-1]}

        if IsLandau:

            dc['Landau_scale'] = cutoffL
            dc['Landau_index'] = iL

        return dc

    def check_boundedness(
            self, scale=None, num=100000,
            loop_level=2):
        """
        Check against tree-level bounded-from-below
        constraints.

        Args:
            scale (float): Energy scale up to
                which the constraints should be
                applied. If not given, the check
                is only performed at the electroweak scale
                $(\mu = v = 246\ \mathrm{GeV})$.
                Value must be larger than $v$.
            num (int): Number of points that are calculated
                for the RGE running when `scale` is given.
            loop_level (int): Loop-level of RGE running.
                Can be set to 1 or 2.

        Returns:
            bool/dict: If `scale=None` then return `True`
                if constraints are fulfilled and `False`
                otherwise. If `scale` is given then return
                dictionary with information about validity
                depending on the energy scale.
        """

        if scale is None:

            return boundedness(
                self.lam1, self.lam2, self.lam3, self.lam4,
                self.lam5, self.lam6, self.lam7, self.lam8)

        else:

            if scale < self.v:

                raise ValueError(
                    'scale has to be larger than initial ' +
                    'scale (=vSM).')

            dcrge, solrge = self.run_to_scale(
                scale, num=num,
                loop_level=loop_level)

            lamsAtScale = solrge[:, 5:13]

            dcL = self._check_lams_for_landau(lamsAtScale, dcrge, num)
            try:
                upper_ind = dcL['Landau_index']
            except KeyError:
                upper_ind = num

            dcB = boundedness_to_scale(
                self,
                lamsAtScale[0:upper_ind, ...],
                dcrge['scale'][0:upper_ind])

            return {**dcL, **dcB}

    def check_metastability(self):
        """
        Check whether the electroweak minimum
        defined by the values of $v_1$, $v_2$
        and $v_S$ is the global minimum of
        the scalar potential.

        Returns:
            bool: `True` when the electroweak
                minimum is the global minimum
                and `False` otherwise.
        """

        self.xew = np.array([
            self.v1, 0., self.v2, 0., self.vs, 0])

        meta = Hom4ps2(self)

        meta.execute()

        self.real_roots = meta.real_roots
        self.local_minima = meta.local_minima
        self.dangerous_minima = meta.dangerous_minima

        if len(self.dangerous_minima) > 0:
            y = False
        else:
            y = True

        return y

    def check_ewpo(self, mode='ST'):
        """
        Check whether the electroweak
        precision observables in terms
        of the oblique parameters are
        predicted to be in agreement
        with the experimental values.

        Args:
            mode (str): When set to `'ST'` then
                a two-dimensional $\chi^2$ fit is
                performed in terms of the $S$ and
                the $T$ parameters.
                When set to `'STU'` then a
                three-dimensional $\chi^2$ fit is
                performed in terms of the $S$, the
                $T$ and the $U$ parameter.

        Returns:
            bool: `False` when the parameter point
                is excluded at the 95% confidence
                level, and `True` otherwise.
        """

        return check_stu(self, mode=mode)

    def eval_darkmatter(self, mode='micro', gauge='unitary'):
        """
        Evaluates the predictions for the
        dark-matter sector.

        Depending on the chosen mode, calls
        external programmes to compute
        predictions for the relic abundace,
        the direct detection of dark matter
        or the indirect detection of dark matter.
        Results are stored in:
        ```
        self.Micromegas['relic']
        self.MadDM['indirect']
        self.directDetection
        ```

        Args:
            mode (str): Selects the computation
                of relic-abundance if set
                to 'micro', indirect-detection
                cross sections if set to
                'maddm' or direct-detection
                cross sections if set to 'direct'.
            gauge (str): Gauge-fixing in the
                CalcHEP files used by MicrOmegas.
                Can be set to `unitary` or `feynman`.
        """
        if mode == 'micro':
            micro = Micromegas(self, gauge=gauge)
            micro.execute()
        if mode == 'maddm':
            maddm = MadDM(self, gauge=gauge)
            maddm.execute()
        if mode == 'direct':
            calc_nucleon_scattering(self)

    def check_darkmatter(self, gauge='unitary'):
        """
        Check whether the parameter point
        is excluded because of a too large
        thermal relic abundance of dark
        matter or due to constraints
        from the indirect-detection constraints
        from dSph observations by the Fermi sattelite.
        Also calculates dark-matter nucleon scattering
        cross section for direct-detection constraints.

        Detailed information about the calculations
        of relic abundance, indirect-detection
        cross sections and direct-detection cross
        sections are stored in:
        ```
        self.Micromegas['relic']
        self.MadDM['indirect']
        self.directDetection
        ```

        Args:
            gauge (str): Gauge-fixing in the
                CalcHEP files used by MicrOmegas.
                Can be set to `unitary` or `feynman`.

        Returns:
            bool: `True` when the constraints (except
                direct detection) are
                respected and `False` otherwise.
        """

        self.eval_darkmatter(mode='micro', gauge=gauge)
        self.eval_darkmatter(mode='maddm', gauge='feynman')
        self.eval_darkmatter(mode='direct')

        # Check if overclosed and indirect dwarf bounds

        yrelic = True
        yindir = True

        if self.Micromegas['relic']['Omegahsq'] > \
            Omegah2_Planck + Omegah2_Planck_Err:

            yrelic = False

        for k, v in self.MadDM['indirect'].items():

            if not k == 'xi':
                if v[1] > v[2]:
                    if v[2] > 0.:
                        yindir = False

        y = yrelic and yindir
        return y

    def eval_center_of_galaxy(self):

        maddm = MadDM(self, gauge='feynman', vrel=1e-3)
        maddm.execute()

    def run_to_scale(
            self, scale, num=1000,
            loop_level=2):
        """
        Computes the parameters of the model
        as a function of the energy scale
        by making use of the running group
        equations.

        Args:
            scale (float): Energy scale up to
                which the constraints should be
                applied. Must be larger than
                the initial scale
                $\mu = v = 246\ \mathrm{GeV}$.
            num (int): Number of points that are calculated
                for the RGE running.
            loop_level (int): Loop-level of RGE running.
                Can be set to 1 or 2.

        Returns:
            (dict, array): First element is a
                dictionary containing
                the parameters as a function of the
                energy scale.The second element
                is a numpy array with the raw
                output of the `odeint` call of SciPy.
        """

        if not loop_level in [1, 2]:
            raise ValueError(
                'Only one-loop (1) and two-loop (2) ' +
                'available. Please choose different ' +
                'value for loop_level.')

        alphas = self.alphas
        MW = self.mW
        MZ = self.mZ
        v = self.v
        MT = self.mt
        MB = self.mb

        g1 = self.g1
        g2 = self.g2
        g3 = self.g3

        if self.yuktype == 2:

            RGE = RGEII

            Yt = np.sqrt(2.) * MT / self.v2
            Yb = np.sqrt(2.) * MB / self.v1

        elif self.yuktype == 1:

            RGE = RGEI

            Yt = np.sqrt(2.) * MT / self.v2
            Yb = np.sqrt(2.) * MB / self.v2

        elif self.yuktype == 3:

            RGE = RGEIII

            Yt = np.sqrt(2.) * MT / self.v2
            Yb = np.sqrt(2.) * MB / self.v2

        elif self.yuktype == 4:

            RGE = RGEIV

            Yt = np.sqrt(2.) * MT / self.v2
            Yb = np.sqrt(2.) * MB / self.v1

        scale_ini = np.log(v)
        scale_end = np.log(scale)

        L1 = self.lam1
        L2 = self.lam2
        L3 = self.lam3
        L4 = self.lam4
        L5 = self.lam5
        L6 = self.lam6
        L7 = self.lam7
        L8 = self.lam8
        m11sq = self.m11sq
        m22sq = self.m22sq
        m12sq = self.m12sq
        mSsq = self.mSsq
        mXsq = self.mXsq

        y = np.array([
            Yt, Yb, g1, g2, g3,
            L1, L2, L3, L4,
            L5, L6, L7, L8,
            m11sq, m22sq, m12sq,
            mSsq, mXsq])

        t = np.linspace(
            scale_ini,
            scale_end,
            num=num)

        if loop_level == 2:
            f = RGE.betafunctions.calc_betas
        else:
            f = RGE.betafunctions.calc_betas_oneloop

        sol = odeint(f, y, t)

        dc = {
            'scale': np.exp(t),
            'Yt': sol[:, 0],
            'Yb': sol[:, 1],
            'g1': sol[:, 2],
            'g2': sol[:, 3],
            'g3': sol[:, 4],
            'L1': sol[:, 5],
            'L2': sol[:, 6],
            'L3': sol[:, 7],
            'L4': sol[:, 8],
            'L5': sol[:, 9],
            'L6': sol[:, 10],
            'L7': sol[:, 11],
            'L8': sol[:, 12],
            'm11sq': sol[:, 13],
            'm22sq': sol[:, 14],
            'm12sq': sol[:, 15],
            'mSsq': sol[:, 16],
            'mSsx': sol[:, 17]}

        self.RunningParas = dc

        return dc, sol

    def check_theory_constraints(
            self, bounded=True,
            meta=False, pert=True,
            pert_cut=EightPi,
            scale=None, num=10000,
            loop_level=2):
        """
        Wrapper function that checks
        gainst all implemented theoretical
        constraints.

        Results regarding each constraint
        are stored in the attributes:
        ```
        self.boundedness
        self.pertUni
        self.EWminIsGlobal
        ```

        Args:
            bounded (bool): Whether the
                bounded-from-below constraints
                should be checked against
            meta (bool): Whether it should be
                checked if the EW minimum is the
                global minium.
            pert (bool): Whether the tree-level
                perturbative unitarity constraints
                should be checked against.
            pert_cut (float): Argument `cutoff` given
                to `check_tree_pert_uni`.
            scale (float): Argument `scale` given to
                `check_tree_pert_uni` and
                `check_boundedness`.
            num (int): Argument `num` given to
                `check_tree_pert_uni` and
                `check_boundedness`.
            loop_level (int): Argument `loop_level`
                given to
                `check_tree_pert_uni` and
                `check_boundedness`.
        """

        # Wrapper function for theory constraints
        # -> Fast checks first

        if not scale is None:

            # If not bounded or pert at mu = v
            #  -> reduce num to save time
            ch1 = self.check_boundedness()
            ch2 = self.check_tree_pert_uni()
            if not (ch1 and ch2):
                num = 40

            if bounded and not pert:

                self.boundedness = self.check_boundedness(
                    scale=scale,
                    num=num,
                    loop_level=loop_level)

            if pert and not bounded:

                self.pertUni = self.check_tree_pert_uni(
                    scale=scale,
                    num=num,
                    loop_level=loop_level,
                    cutoff=pert_cut)

            if pert and bounded:

                # If both then this function to not make running twice

                y = self.check_bounded_pertuni(
                    scale=scale,
                    num=num,
                    loop_level=loop_level,
                    cutoff=pert_cut)

                self.boundedness = y[0]
                self.pertUni = y[1]

        else:

            if bounded:

                self.boundedness = self.check_boundedness()

            if pert:

                self.pertUni = self.check_tree_pert_uni()

        if meta:

            self.EWminIsGlobal = self.check_metastability()


    def check_bounded_pertuni(
            self, scale, num=100000,
            loop_level=2, cutoff=EightPi):

        trPrU.CUTOFF = cutoff

        if scale < self.v:

            raise ValueError(
                'scale has to be larger than initial ' +
                'scale (=vSM).')

        dcrge, solrge = self.run_to_scale(
            scale, num=num,
            loop_level=loop_level)

        lamsAtScale = solrge[:, 5:13]

        dcL = self._check_lams_for_landau(lamsAtScale, dcrge, num)
        try:
            upper_ind = dcL['Landau_index']
        except KeyError:
            upper_ind = num

        dcB = boundedness_to_scale(
            self,
            lamsAtScale[0:upper_ind, ...],
            dcrge['scale'][0:upper_ind])

        dcP = trPrU.pertuni_to_scale(
            self,
            lamsAtScale[0:upper_ind, ...],
            dcrge['scale'][0:upper_ind])

        return {**dcL, **dcB}, {**dcL, **dcP}

    def check_collider_constraints(self):
        """
        Calls HiggsTools to test the parameter
        point against the cross-section limits
        from searches for additional Higgs bosons
        (HiggsBounds) and to perform a $\chi^2$-fit
        to the cross-section measurements of the
        discovered Higgs boson at $125~\mathrm{GeV}$
        (HiggsSignals).

        Results are stored in:

        ```
        self.hb_result
        self.HT['HB']
        self.HT['HS']
        ```

        The object hb_result is the object returned
        by HiggsTools when calling the HiggsBounds test.
        The dictionary HT contains the most relevant
        information from the HiggsBounds
        and the HiggsSignals analysis under the keys
        'HB' and 'HS', respectively.
        """
        HiggsTools.check_point(self)

__init__(dc, debug=True)

Initialize an instance of a parameter point given a set of input parameters. The parameters have to be given as the dictionary dc. There are two possible input formats:

Angle input, for instance:

dc = {
    'type': 2,
    'tb': 1.26,
    'al1': 1.27,
    'al2': -1.08,
    'al3': -1.24,
    'mH1': 96.5,
    'mH2': 125.09,
    'mH3': 535.86,
    'mXi': 266.0,
    'mA': 712.578,
    'mHp': 737.8,
    'vs': 272.0,
    'm12sq': 80644.0}

Lambda input, for instance:

dc = {
    'type': 2,
    'tb': 0.91655,
    'mXi': 134.03,
    'vs': 64.987,
    'm12sq': 1.7696e5,
    'lam1': 1.5297,
    'lam2': 1.2074,
    'lam3': 1.5741,
    'lam4': 5.3967,
    'lam5': -7.8556,
    'lam6': 6.0689,
    'lam7': 0.80378,
    'lam8': -0.83745}

Parameters:

Name	Type	Description	Default
dc	dict	Dictionary with input values for free parameters, either in angle- or in lambda-input format.	required
debug	bool	If set to True then the debugging functions are called.	True

Source code in s2hdmTools/paramPoint.py

def __init__(self, dc, debug=True):
    """
    Initialize an instance of a parameter
    point given a set of input parameters.
    The parameters have to be given as
    the dictionary dc. There are two
    possible input formats:

    Angle input, for instance:
    ```
    dc = {
        'type': 2,
        'tb': 1.26,
        'al1': 1.27,
        'al2': -1.08,
        'al3': -1.24,
        'mH1': 96.5,
        'mH2': 125.09,
        'mH3': 535.86,
        'mXi': 266.0,
        'mA': 712.578,
        'mHp': 737.8,
        'vs': 272.0,
        'm12sq': 80644.0}
    ```

    Lambda input, for instance:
    ```
    dc = {
        'type': 2,
        'tb': 0.91655,
        'mXi': 134.03,
        'vs': 64.987,
        'm12sq': 1.7696e5,
        'lam1': 1.5297,
        'lam2': 1.2074,
        'lam3': 1.5741,
        'lam4': 5.3967,
        'lam5': -7.8556,
        'lam6': 6.0689,
        'lam7': 0.80378,
        'lam8': -0.83745}
    ```

    Args:
        dc (dict): Dictionary with input values
            for free parameters, either in angle-
            or in lambda-input format.
        debug (bool): If set to `True` then the
            debugging functions are called.
    """

    self.dc = dc
    self.read_constants()
    self.read_input()
    if debug:
        self.check_masses()
    self.debug = debug

calculate_branching_ratios()

Calculate the branching ratios of the Higgs bosons.

Branching ratios are stored in:

self.b_H1
self.b_H2
self.b_H3
self.b_A
self.b_Hp

Total widths are stored in:

self.wTot

Source code in s2hdmTools/paramPoint.py

def calculate_branching_ratios(self):
    """
    Calculate the branching ratios of
    the Higgs bosons.

    Branching ratios are stored in:

    ```
    self.b_H1
    self.b_H2
    self.b_H3
    self.b_A
    self.b_Hp
    ```

    Total widths are stored in:

    ```
    self.wTot
    ```
    """

    n2hdecay = CallN2HDecay(self)

    invs = calc_gammas_inv(self)

    self.b_A = n2hdecay.b_A
    self.b_Hp = n2hdecay.b_Hp
    self.b_t = n2hdecay.b_t

    gamH1 = {}
    wH1 = n2hdecay.w['H1']
    for key in n2hdecay.b_H1:
        gamH1[key] = n2hdecay.b_H1[key] * wH1
    gamH1['XX'] = invs[0]
    wH1 += invs[0]
    self.b_H1 = {}
    for key in gamH1:
        self.b_H1[key] = gamH1[key] / wH1
    if self.debug:
        if abs(sum(self.b_H1.values()) - 1.) > 1e-6:
            raise RuntimeError(
                'Branching ratios of H1 do not add up to one.')

    gamH2 = {}
    wH2 = n2hdecay.w['H2']
    for key in n2hdecay.b_H2:
        gamH2[key] = n2hdecay.b_H2[key] * wH2
    gamH2['XX'] = invs[1]
    wH2 += invs[1]
    self.b_H2 = {}
    for key in gamH2:
        self.b_H2[key] = gamH2[key] / wH2
    if self.debug:
        if abs(sum(self.b_H2.values()) - 1.) > 1e-6:
            raise RuntimeError(
                'Branching ratios of H2 do not add up to one.')

    gamH3 = {}
    wH3 = n2hdecay.w['H3']
    for key in n2hdecay.b_H3:
        gamH3[key] = n2hdecay.b_H3[key] * wH3
    gamH3['XX'] = invs[2]
    wH3 += invs[2]
    self.b_H3 = {}
    for key in gamH3:
        self.b_H3[key] = gamH3[key] / wH3
    if self.debug:
        if abs(sum(self.b_H3.values()) - 1.) > 1e-6:
            raise RuntimeError(
                'Branching ratios of H3 do not add up to one.')

    self.wTot = {
        'A': n2hdecay.w['A'],
        'H1': wH1,
        'H2': wH2,
        'H3': wH3,
        'Hp': n2hdecay.w['Hp'],
        't': n2hdecay.w['t']}

check_boundedness(scale=None, num=100000, loop_level=2)

Check against tree-level bounded-from-below constraints.

Parameters:

Name	Type	Description	Default
scale	float	Energy scale up to which the constraints should be applied. If not given, the check is only performed at the electroweak scale \((\mu = v = 246\ \mathrm{GeV})\). Value must be larger than \(v\).	None
num	int	Number of points that are calculated for the RGE running when scale is given.	100000
loop_level	int	Loop-level of RGE running. Can be set to 1 or 2.	2

Returns:

Type	Description
	bool/dict: If scale=None then return True if constraints are fulfilled and False otherwise. If scale is given then return dictionary with information about validity depending on the energy scale.

Source code in s2hdmTools/paramPoint.py

def check_boundedness(
        self, scale=None, num=100000,
        loop_level=2):
    """
    Check against tree-level bounded-from-below
    constraints.

    Args:
        scale (float): Energy scale up to
            which the constraints should be
            applied. If not given, the check
            is only performed at the electroweak scale
            $(\mu = v = 246\ \mathrm{GeV})$.
            Value must be larger than $v$.
        num (int): Number of points that are calculated
            for the RGE running when `scale` is given.
        loop_level (int): Loop-level of RGE running.
            Can be set to 1 or 2.

    Returns:
        bool/dict: If `scale=None` then return `True`
            if constraints are fulfilled and `False`
            otherwise. If `scale` is given then return
            dictionary with information about validity
            depending on the energy scale.
    """

    if scale is None:

        return boundedness(
            self.lam1, self.lam2, self.lam3, self.lam4,
            self.lam5, self.lam6, self.lam7, self.lam8)

    else:

        if scale < self.v:

            raise ValueError(
                'scale has to be larger than initial ' +
                'scale (=vSM).')

        dcrge, solrge = self.run_to_scale(
            scale, num=num,
            loop_level=loop_level)

        lamsAtScale = solrge[:, 5:13]

        dcL = self._check_lams_for_landau(lamsAtScale, dcrge, num)
        try:
            upper_ind = dcL['Landau_index']
        except KeyError:
            upper_ind = num

        dcB = boundedness_to_scale(
            self,
            lamsAtScale[0:upper_ind, ...],
            dcrge['scale'][0:upper_ind])

        return {**dcL, **dcB}

check_collider_constraints()

Calls HiggsTools to test the parameter point against the cross-section limits from searches for additional Higgs bosons (HiggsBounds) and to perform a \(\chi^2\)-fit to the cross-section measurements of the discovered Higgs boson at \(125~\mathrm{GeV}\) (HiggsSignals).

Results are stored in:

self.hb_result
self.HT['HB']
self.HT['HS']

The object hb_result is the object returned by HiggsTools when calling the HiggsBounds test. The dictionary HT contains the most relevant information from the HiggsBounds and the HiggsSignals analysis under the keys 'HB' and 'HS', respectively.

Source code in s2hdmTools/paramPoint.py

def check_collider_constraints(self):
    """
    Calls HiggsTools to test the parameter
    point against the cross-section limits
    from searches for additional Higgs bosons
    (HiggsBounds) and to perform a $\chi^2$-fit
    to the cross-section measurements of the
    discovered Higgs boson at $125~\mathrm{GeV}$
    (HiggsSignals).

    Results are stored in:

    ```
    self.hb_result
    self.HT['HB']
    self.HT['HS']
    ```

    The object hb_result is the object returned
    by HiggsTools when calling the HiggsBounds test.
    The dictionary HT contains the most relevant
    information from the HiggsBounds
    and the HiggsSignals analysis under the keys
    'HB' and 'HS', respectively.
    """
    HiggsTools.check_point(self)

check_darkmatter(gauge='unitary')

Check whether the parameter point is excluded because of a too large thermal relic abundance of dark matter or due to constraints from the indirect-detection constraints from dSph observations by the Fermi sattelite. Also calculates dark-matter nucleon scattering cross section for direct-detection constraints.

Detailed information about the calculations of relic abundance, indirect-detection cross sections and direct-detection cross sections are stored in:

self.Micromegas['relic']
self.MadDM['indirect']
self.directDetection

Parameters:

Name	Type	Description	Default
gauge	str	Gauge-fixing in the CalcHEP files used by MicrOmegas. Can be set to unitary or feynman.	'unitary'

Returns:

Name	Type	Description
bool		True when the constraints (except direct detection) are respected and False otherwise.

Source code in s2hdmTools/paramPoint.py

def check_darkmatter(self, gauge='unitary'):
    """
    Check whether the parameter point
    is excluded because of a too large
    thermal relic abundance of dark
    matter or due to constraints
    from the indirect-detection constraints
    from dSph observations by the Fermi sattelite.
    Also calculates dark-matter nucleon scattering
    cross section for direct-detection constraints.

    Detailed information about the calculations
    of relic abundance, indirect-detection
    cross sections and direct-detection cross
    sections are stored in:
    ```
    self.Micromegas['relic']
    self.MadDM['indirect']
    self.directDetection
    ```

    Args:
        gauge (str): Gauge-fixing in the
            CalcHEP files used by MicrOmegas.
            Can be set to `unitary` or `feynman`.

    Returns:
        bool: `True` when the constraints (except
            direct detection) are
            respected and `False` otherwise.
    """

    self.eval_darkmatter(mode='micro', gauge=gauge)
    self.eval_darkmatter(mode='maddm', gauge='feynman')
    self.eval_darkmatter(mode='direct')

    # Check if overclosed and indirect dwarf bounds

    yrelic = True
    yindir = True

    if self.Micromegas['relic']['Omegahsq'] > \
        Omegah2_Planck + Omegah2_Planck_Err:

        yrelic = False

    for k, v in self.MadDM['indirect'].items():

        if not k == 'xi':
            if v[1] > v[2]:
                if v[2] > 0.:
                    yindir = False

    y = yrelic and yindir
    return y

check_ewpo(mode='ST')

Check whether the electroweak precision observables in terms of the oblique parameters are predicted to be in agreement with the experimental values.

Parameters:

Name	Type	Description	Default
mode	str	When set to 'ST' then a two-dimensional \(\chi^2\) fit is performed in terms of the \(S\) and the \(T\) parameters. When set to 'STU' then a three-dimensional \(\chi^2\) fit is performed in terms of the \(S\), the \(T\) and the \(U\) parameter.	'ST'

Returns:

Name	Type	Description
bool		False when the parameter point is excluded at the 95% confidence level, and True otherwise.

Source code in s2hdmTools/paramPoint.py

def check_ewpo(self, mode='ST'):
    """
    Check whether the electroweak
    precision observables in terms
    of the oblique parameters are
    predicted to be in agreement
    with the experimental values.

    Args:
        mode (str): When set to `'ST'` then
            a two-dimensional $\chi^2$ fit is
            performed in terms of the $S$ and
            the $T$ parameters.
            When set to `'STU'` then a
            three-dimensional $\chi^2$ fit is
            performed in terms of the $S$, the
            $T$ and the $U$ parameter.

    Returns:
        bool: `False` when the parameter point
            is excluded at the 95% confidence
            level, and `True` otherwise.
    """

    return check_stu(self, mode=mode)

check_metastability()

Check whether the electroweak minimum defined by the values of \(v_1\), \(v_2\) and \(v_S\) is the global minimum of the scalar potential.

Returns:

Name	Type	Description
bool		True when the electroweak minimum is the global minimum and False otherwise.

Source code in s2hdmTools/paramPoint.py

def check_metastability(self):
    """
    Check whether the electroweak minimum
    defined by the values of $v_1$, $v_2$
    and $v_S$ is the global minimum of
    the scalar potential.

    Returns:
        bool: `True` when the electroweak
            minimum is the global minimum
            and `False` otherwise.
    """

    self.xew = np.array([
        self.v1, 0., self.v2, 0., self.vs, 0])

    meta = Hom4ps2(self)

    meta.execute()

    self.real_roots = meta.real_roots
    self.local_minima = meta.local_minima
    self.dangerous_minima = meta.dangerous_minima

    if len(self.dangerous_minima) > 0:
        y = False
    else:
        y = True

    return y

check_theory_constraints(bounded=True, meta=False, pert=True, pert_cut=EightPi, scale=None, num=10000, loop_level=2)

Wrapper function that checks gainst all implemented theoretical constraints.

Results regarding each constraint are stored in the attributes:

self.boundedness
self.pertUni
self.EWminIsGlobal

Parameters:

Name	Type	Description	Default
bounded	bool	Whether the bounded-from-below constraints should be checked against	True
meta	bool	Whether it should be checked if the EW minimum is the global minium.	False
pert	bool	Whether the tree-level perturbative unitarity constraints should be checked against.	True
pert_cut	float	Argument cutoff given to check_tree_pert_uni.	EightPi
scale	float	Argument scale given to check_tree_pert_uni and check_boundedness.	None
num	int	Argument num given to check_tree_pert_uni and check_boundedness.	10000
loop_level	int	Argument loop_level given to check_tree_pert_uni and check_boundedness.	2

Source code in s2hdmTools/paramPoint.py

def check_theory_constraints(
        self, bounded=True,
        meta=False, pert=True,
        pert_cut=EightPi,
        scale=None, num=10000,
        loop_level=2):
    """
    Wrapper function that checks
    gainst all implemented theoretical
    constraints.

    Results regarding each constraint
    are stored in the attributes:
    ```
    self.boundedness
    self.pertUni
    self.EWminIsGlobal
    ```

    Args:
        bounded (bool): Whether the
            bounded-from-below constraints
            should be checked against
        meta (bool): Whether it should be
            checked if the EW minimum is the
            global minium.
        pert (bool): Whether the tree-level
            perturbative unitarity constraints
            should be checked against.
        pert_cut (float): Argument `cutoff` given
            to `check_tree_pert_uni`.
        scale (float): Argument `scale` given to
            `check_tree_pert_uni` and
            `check_boundedness`.
        num (int): Argument `num` given to
            `check_tree_pert_uni` and
            `check_boundedness`.
        loop_level (int): Argument `loop_level`
            given to
            `check_tree_pert_uni` and
            `check_boundedness`.
    """

    # Wrapper function for theory constraints
    # -> Fast checks first

    if not scale is None:

        # If not bounded or pert at mu = v
        #  -> reduce num to save time
        ch1 = self.check_boundedness()
        ch2 = self.check_tree_pert_uni()
        if not (ch1 and ch2):
            num = 40

        if bounded and not pert:

            self.boundedness = self.check_boundedness(
                scale=scale,
                num=num,
                loop_level=loop_level)

        if pert and not bounded:

            self.pertUni = self.check_tree_pert_uni(
                scale=scale,
                num=num,
                loop_level=loop_level,
                cutoff=pert_cut)

        if pert and bounded:

            # If both then this function to not make running twice

            y = self.check_bounded_pertuni(
                scale=scale,
                num=num,
                loop_level=loop_level,
                cutoff=pert_cut)

            self.boundedness = y[0]
            self.pertUni = y[1]

    else:

        if bounded:

            self.boundedness = self.check_boundedness()

        if pert:

            self.pertUni = self.check_tree_pert_uni()

    if meta:

        self.EWminIsGlobal = self.check_metastability()

check_tree_pert_uni(scale=None, num=100000, loop_level=2, cutoff=EightPi)

Check against tree-level perturbative unitarity constraints.

Parameters:

Name	Type	Description	Default
scale	float	Energy scale up to which the constraints should be applied. If not given, the check is only performed at the electroweak scale \((\mu = v = 246\ \mathrm{GeV})\). Value must be larger than \(v\).	None
num	int	Number of points that are calculated for the RGE running when scale is given.	100000
loop_level	int	Loop-level of RGE running. Can be set to 1 or 2.	2
cutoff	float	Defines the upper limit for the eigenvalues of the scattering matrix.	EightPi

Returns:

Type	Description
	bool/dict: If scale=None then return True if constraints are fulfilled and False otherwise. If scale is given then return dictionary with information about validity depending on the energy scale.

Source code in s2hdmTools/paramPoint.py

def check_tree_pert_uni(
        self, scale=None, num=100000,
        loop_level=2,
        cutoff=EightPi):
    """
    Check against tree-level perturbative
    unitarity constraints.

    Args:
        scale (float): Energy scale up to
            which the constraints should be
            applied. If not given, the check
            is only performed at the electroweak scale
            $(\mu = v = 246\ \mathrm{GeV})$.
            Value must be larger than $v$.
        num (int): Number of points that are calculated
            for the RGE running when `scale` is given.
        loop_level (int): Loop-level of RGE running.
            Can be set to 1 or 2.
        cutoff (float): Defines
            the upper limit for the eigenvalues
            of the scattering matrix.

    Returns:
        bool/dict: If `scale=None` then return `True`
            if constraints are fulfilled and `False`
            otherwise. If `scale` is given then return
            dictionary with information about validity
            depending on the energy scale.
    """

    trPrU.CUTOFF = cutoff

    if scale is None:

        raw = trPrU.pertuni(
            self.lam1, self.lam2, self.lam3, self.lam4,
            self.lam5, self.lam6, self.lam7, self.lam8)

        self.tree_pert_uni_vals = raw[0]

        return raw[1]

    else:

        if scale < self.v:

            raise ValueError(
                'scale has to be larger than initial ' +
                'scale (=vSM).')

        dcrge, solrge = self.run_to_scale(
            scale, num=num,
            loop_level=loop_level)

        lamsAtScale = solrge[:, 5:13]

        dcL = self._check_lams_for_landau(lamsAtScale, dcrge, num)
        try:
            upper_ind = dcL['Landau_index']
        except KeyError:
            upper_ind = num

        dcP = trPrU.pertuni_to_scale(
            self,
            lamsAtScale[0:upper_ind, ...],
            dcrge['scale'][0:upper_ind])

        return {**dcL, **dcP}

eval_darkmatter(mode='micro', gauge='unitary')

Evaluates the predictions for the dark-matter sector.

Depending on the chosen mode, calls external programmes to compute predictions for the relic abundace, the direct detection of dark matter or the indirect detection of dark matter. Results are stored in:

self.Micromegas['relic']
self.MadDM['indirect']
self.directDetection

Parameters:

Name	Type	Description	Default
mode	str	Selects the computation of relic-abundance if set to 'micro', indirect-detection cross sections if set to 'maddm' or direct-detection cross sections if set to 'direct'.	'micro'
gauge	str	Gauge-fixing in the CalcHEP files used by MicrOmegas. Can be set to unitary or feynman.	'unitary'

Source code in s2hdmTools/paramPoint.py

def eval_darkmatter(self, mode='micro', gauge='unitary'):
    """
    Evaluates the predictions for the
    dark-matter sector.

    Depending on the chosen mode, calls
    external programmes to compute
    predictions for the relic abundace,
    the direct detection of dark matter
    or the indirect detection of dark matter.
    Results are stored in:
    ```
    self.Micromegas['relic']
    self.MadDM['indirect']
    self.directDetection
    ```

    Args:
        mode (str): Selects the computation
            of relic-abundance if set
            to 'micro', indirect-detection
            cross sections if set to
            'maddm' or direct-detection
            cross sections if set to 'direct'.
        gauge (str): Gauge-fixing in the
            CalcHEP files used by MicrOmegas.
            Can be set to `unitary` or `feynman`.
    """
    if mode == 'micro':
        micro = Micromegas(self, gauge=gauge)
        micro.execute()
    if mode == 'maddm':
        maddm = MadDM(self, gauge=gauge)
        maddm.execute()
    if mode == 'direct':
        calc_nucleon_scattering(self)

run_to_scale(scale, num=1000, loop_level=2)

Computes the parameters of the model as a function of the energy scale by making use of the running group equations.

Parameters:

Name	Type	Description	Default
scale	float	Energy scale up to which the constraints should be applied. Must be larger than the initial scale \(\mu = v = 246\ \mathrm{GeV}\).	required
num	int	Number of points that are calculated for the RGE running.	1000
loop_level	int	Loop-level of RGE running. Can be set to 1 or 2.	2

Returns:

Type	Description
dict, array	First element is a dictionary containing the parameters as a function of the energy scale.The second element is a numpy array with the raw output of the odeint call of SciPy.

Source code in s2hdmTools/paramPoint.py

def run_to_scale(
        self, scale, num=1000,
        loop_level=2):
    """
    Computes the parameters of the model
    as a function of the energy scale
    by making use of the running group
    equations.

    Args:
        scale (float): Energy scale up to
            which the constraints should be
            applied. Must be larger than
            the initial scale
            $\mu = v = 246\ \mathrm{GeV}$.
        num (int): Number of points that are calculated
            for the RGE running.
        loop_level (int): Loop-level of RGE running.
            Can be set to 1 or 2.

    Returns:
        (dict, array): First element is a
            dictionary containing
            the parameters as a function of the
            energy scale.The second element
            is a numpy array with the raw
            output of the `odeint` call of SciPy.
    """

    if not loop_level in [1, 2]:
        raise ValueError(
            'Only one-loop (1) and two-loop (2) ' +
            'available. Please choose different ' +
            'value for loop_level.')

    alphas = self.alphas
    MW = self.mW
    MZ = self.mZ
    v = self.v
    MT = self.mt
    MB = self.mb

    g1 = self.g1
    g2 = self.g2
    g3 = self.g3

    if self.yuktype == 2:

        RGE = RGEII

        Yt = np.sqrt(2.) * MT / self.v2
        Yb = np.sqrt(2.) * MB / self.v1

    elif self.yuktype == 1:

        RGE = RGEI

        Yt = np.sqrt(2.) * MT / self.v2
        Yb = np.sqrt(2.) * MB / self.v2

    elif self.yuktype == 3:

        RGE = RGEIII

        Yt = np.sqrt(2.) * MT / self.v2
        Yb = np.sqrt(2.) * MB / self.v2

    elif self.yuktype == 4:

        RGE = RGEIV

        Yt = np.sqrt(2.) * MT / self.v2
        Yb = np.sqrt(2.) * MB / self.v1

    scale_ini = np.log(v)
    scale_end = np.log(scale)

    L1 = self.lam1
    L2 = self.lam2
    L3 = self.lam3
    L4 = self.lam4
    L5 = self.lam5
    L6 = self.lam6
    L7 = self.lam7
    L8 = self.lam8
    m11sq = self.m11sq
    m22sq = self.m22sq
    m12sq = self.m12sq
    mSsq = self.mSsq
    mXsq = self.mXsq

    y = np.array([
        Yt, Yb, g1, g2, g3,
        L1, L2, L3, L4,
        L5, L6, L7, L8,
        m11sq, m22sq, m12sq,
        mSsq, mXsq])

    t = np.linspace(
        scale_ini,
        scale_end,
        num=num)

    if loop_level == 2:
        f = RGE.betafunctions.calc_betas
    else:
        f = RGE.betafunctions.calc_betas_oneloop

    sol = odeint(f, y, t)

    dc = {
        'scale': np.exp(t),
        'Yt': sol[:, 0],
        'Yb': sol[:, 1],
        'g1': sol[:, 2],
        'g2': sol[:, 3],
        'g3': sol[:, 4],
        'L1': sol[:, 5],
        'L2': sol[:, 6],
        'L3': sol[:, 7],
        'L4': sol[:, 8],
        'L5': sol[:, 9],
        'L6': sol[:, 10],
        'L7': sol[:, 11],
        'L8': sol[:, 12],
        'm11sq': sol[:, 13],
        'm22sq': sol[:, 14],
        'm12sq': sol[:, 15],
        'mSsq': sol[:, 16],
        'mSsx': sol[:, 17]}

    self.RunningParas = dc

    return dc, sol