example_gm2scan.py.in

Go to the documentation of this file.

#!/usr/bin/env python
 
# ====================================================================
# This file is part of GM2Calc.
#
# GM2Calc is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License,
# or (at your option) any later version.
#
# GM2Calc is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with GM2Calc.  If not, see
# <http://www.gnu.org/licenses/>.
# ====================================================================
 
from __future__ import print_function
from gm2_python_interface import *
 
cppyy.include(os.path.join("gm2calc","gm2_1loop.hpp"))
cppyy.include(os.path.join("gm2calc","gm2_2loop.hpp"))
cppyy.include(os.path.join("gm2calc","gm2_uncertainty.hpp"))
cppyy.include(os.path.join("gm2calc","gm2_error.hpp"))
cppyy.include(os.path.join("gm2calc","MSSMNoFV_onshell.hpp"))
 
cppyy.load_library("libgm2calc")
 
# Load data types
from cppyy.gbl import std
from cppyy.gbl import gm2calc
from cppyy.gbl.gm2calc import Error 
 
def setup():
    # load data types
    from cppyy.gbl import Eigen
    # Shorten class call name
    from cppyy.gbl.Eigen import Matrix3d
    
    model = gm2calc.MSSMNoFV_onshell()
    
    Pi = 3.141592653589793
    # Outer Matrix3cd is to convert type from CwiseNullaryOp to Matrix
    UnitMatrix = Matrix3d(Matrix3d().Identity())
    # __mul__ is not defined for Matrix type
    UnitMatrix.__imul__(500*500)
    
    # fill SM parameters
    model.set_alpha_MZ(0.0077552)                # 1L
    model.set_alpha_thompson(0.00729735)         # 2L
    model.set_g3(std.sqrt(4. * Pi * 0.1184))    # 2L
    model.get_physical().MFt   = 173.34          # 2L
    model.get_physical().MFb   = 4.18            # 2L, mb(mb) MS-bar
    model.get_physical().MFm   = 0.1056583715    # 1L
    model.get_physical().MFtau = 1.777           # 2L
    model.get_physical().MVWm  = 80.385          # 1L
    model.get_physical().MVZ   = 91.1876         # 1L
    
    # fill DR-bar parameters
    model.set_TB(10)                       # 1L
    model.set_Ae(1,1,0)                    # 1L
 
    # fill on-shell parameters
    model.set_Mu(350)                      # 1L
    model.set_MassB(150)                   # 1L
    model.set_MassWB(300)                  # 1L
    model.set_MassG(1000)                  # 2L
    model.set_mq2(UnitMatrix)              # 2L
    model.set_ml2(UnitMatrix)              # 1L(smuon)/2L
    model.set_md2(UnitMatrix)              # 2L
    model.set_mu2(UnitMatrix)              # 2L
    model.set_me2(UnitMatrix)              # 1L(smuon)/2L
    model.set_Au(2,2,0)                    # 2L
    model.set_Ad(2,2,0)                    # 2L
    model.set_Ae(2,2,0)                    # 2L
    model.set_MA0(1500)                    # 2L
    model.set_scale(454.7)                 # 2L
    
    return model;
 
tanb_start = 2.
tanb_stop = 100.
nsteps = 100
string = ["tan(beta)","amu","uncertainty","error"]
print("# {0:>14s} {1:>16s} {2:>16s} {3:>16s}".format(*string))
for n in range(0,nsteps):
    amu = 0.
    delta_amu = 0.
    tanb = tanb_start + (tanb_stop - tanb_start) * n / nsteps
    error = std.string()
    
    model = setup()
    model.set_TB(tanb)
    
    try:
        model.calculate_masses()
        amu = gm2calc.calculate_amu_1loop(model) + gm2calc.calculate_amu_2loop(model)
        delta_amu = gm2calc.calculate_uncertainty_amu_2loop(model)
    except gm2calc.Error as e:
        error = "# " + str(e.what())
        amu = std.numeric_limits['double'].signaling_NaN()
        delta_amu = std.numeric_limits['double'].signaling_NaN()
 
    print("{0:>16.8e} {1:>16.8e} {2:>16.8e}".format(tanb,amu,delta_amu)+error.c_str())