#!/usr/bin/env python
# This software is covered by the GNU General Public License 3.
# Copyright (C) 2024-2025 Killian Bouzoud, Jacopo Ghiglieri, Greg Jackson# Add Authors, LICENSE, etc. here
# This script controls the execution of the analytical part of AutoTherm
import os
import subprocess
import sys
import configparser
import re
import warnings
import json
[docs]
def check_particle(thepart):
"""
Checks if a particle name of the form X[n] is correct.
It considers the following particles incorrect:
* Weyl fields
* Rarita-Schwinger fields (FeynRules can generate the associated Feynman rules but they can't be inserted into FeynArts topologies AFAIK)
* Ghost fields
The function only checks the "X" part, not the number n
:param str thepart: Name of the particle to check
:return: Whether the particle name is valid or not
:rtype: bool
"""
supported = r"(S|F|V|T|-S|-F|-V|-T)\[.*\]"
unsupported = r"(W|R|RW|U)\[.*\]"
if re.match(supported,thepart) is not None:
return True
elif re.match(unsupported,thepart) is not None:
warnings.warn("Particle "+thepart+" is recognized but not supported.",RuntimeWarning)
return False
else:
return False
[docs]
class LagrangianError(Exception):
"""
This is a class to have a custom error message when we check the Lagragian and find a Lagrangian term without a ``:=``
:param str message: Apparently required for an ``Exception`` object.
:param str lterm: The problematic Lagrangian term, or rather the first line on which it is
:param int lineno: Line number of the problematic term
:param str other: Additional text
"""
def __init__(self,message,lterm,lineno,other):
super().__init__(message)
self.lterm = lterm
self.lineno = lineno
# I can't seem to be able to access self.message in __str__ so I need to do this
self.other = other
def __str__(self):
return f"Problematic Lagrangian term -> Line {self.lineno}: {self.lterm}\n{self.other}"
[docs]
def check_for_equal(filename):
"""
Read the model file and check the Lagrangian terms.
For now, it is assumed that all terms in the Lagrangian begin with a "L".
:param str filename: The path to the model file, as defined by the user in the configuration
:rtype: NoneType
:raises LagrangianError: if a Lagrangian term without a ``:=`` is found
"""
# Count the lines to help the user locate the error
counter = 0
with open(os.path.expanduser(filename),'r') as themodel:
for line in themodel:
counter+=1
# HACK: For the time being, only lines beginning with a "L" (like "Lterm := aaaaaaaaaaa") are checked
if re.match(r"^L",line):
# If we can find a :=, everything is fine
if re.match(r"(.*):=(.*)",line):
continue
# Raise an error with the correct arguments
else:
raise LagrangianError("",line.strip(),counter,"This Lagrangian term is defined without a := operator.\nThis might cause problems when contracting indices.")
[docs]
def do_the_thermal_bath(thebath):
"""
Construct the thermal bath.
The outputs are either:
#. The SM thermal bath if ``thebath`` is empty
#. The SM thermal bath plus and/or minus a few particles if the correct syntax is used
#. Just ``thebath`` as a list if the previous 2 cases are not met
:param str thebath: The thermal bath as configured by the user in the configuration file
:return: The thermal bath, as list, correctly parsed by AutoTherm
:rtype: list[str]
"""
# The SM thermal bath
thermal_bath_SM = "V[1] V[2] V[3] F[1] F[2] F[3] F[4] F[5] -F[1] -F[2] -F[3] -F[4] -F[5] S[1] -S[1]".split(' ')
# If the bath is empty, assume SM
if not thebath:
return thermal_bath_SM
# If the bath is not empty, we need to determine if we make the bath equal to the argument in the function
# Or if we apply a modification
else:
# List of particles to add to the SM thermal bath
to_add = []
# Regex to see that
# The "+" needs to be escaped as it's a regex operator
to_add_re = r"(S|F|V|T|-S|-F|-V|-T)\[.*\]\+$"
# List of particles to remove from the SM thermal bath
to_remove = []
# Regex to see that
to_remove_re = r"(S|F|V|T|-S|-F|-V|-T)\[.*\]-$"
for part in thebath.split(' '):
# If there is a match, add the particle without the ending "+" or "-" to the adequate list
if re.match(to_add_re,part):
to_add.append(part[:-1])
if re.match(to_remove_re,part):
to_remove.append(part[:-1])
# If there was no modification to do, return the bath as stated
if (len(to_add)==0 and len(to_remove)==0):
print("There was no particle to add to or remove from the SM thermal bath")
return thebath.split(' ')
else:
# In the event that there is a list such as F[1]- V[1]+ S[3]
# The S[3] won't be in to_add or in to_remove and will be ignored
# So we need to warn the user
if(len(to_add)+len(to_remove) != len(thebath.split(' '))):
warnings.warn("It seems a particle was ignored. \n Execution will continue regardless.",RuntimeWarning)
# Add the particles that need to be added
for part in to_add:
# Security if the particle is already in the thermal bath
if part in thermal_bath_SM:
print("The particle "+part+" is already in the SM thermal bath. Ignoring.")
continue
else:
thermal_bath_SM.append(part)
# Remove the particles that need to be removed
for part in to_remove:
# Security if there is no particle to remove
if (not part in thermal_bath_SM):
print("The particle "+part+" is already absent from the SM thermal bath. Ignoring.")
continue
else:
thermal_bath_SM.remove(part)
return thermal_bath_SM
[docs]
def do_the_dict(thefile):
"""
Compute the dictionary that the class ``NumRate`` requires.
:param str thefile: ``.json`` file containing the result. Generated by ``processes.wl``
:return: Parsed matrix element squared, couplings and thermal masses (in this order)
:rtype: dict, dict, tuple
:raises RuntimeError: if the result file can't be found for whatever reason
:raises json.JSONDecodeError: if the ``.json`` file is not formatted correctly
:raises AssertionError: if the nature of the particles are incorrect
"""
# FIXME add explanation of what happens when fermionflag is on
# Check if the file does exist
if not os.path.isfile(thefile):
raise RuntimeError
matdict = {}
coupldict = {}
masses = []
flagsdict = {}
with open(thefile) as data_file:
data = json.load(data_file)
flagsdict = data['flags']
fermionflag=flagsdict['fermionflag']
# Format the dictionary for the matrix element squared
for key,value in data['mat2'].items():
# Convert string to tuple
# Taken from https://www.geeksforgeeks.org/python-convert-string-to-tuple/
stats = eval(key)
# Check if all elements in stats are either +1 or -1
assert [i**2 for i in stats] == [1,1,1]
the_mat2 = value
# Update the dictionary
# Change the order from s1,s2,t1 to t1,s1,s2 to match NumRate
matdict.update({(stats[2],stats[0],stats[1]):the_mat2})
# Format the list for the masses
for key,value in data['masses'].items():
# The key is the particle name which isn't used later
masses.append(value)
# NumRate uses a tuple, not a list
masses = tuple(masses)
# Format the dictionary for the parameters
for key,value in data['couplings'].items():
# The keys remain the same
# The values need to be a tuple of strings
if value == '':
mytuple = ()
else:
mytuple = tuple(value.split(','))
coupldict.update({key:mytuple})
if fermionflag == False:
#in this case we have a single return dict
if len(masses)==1:
masses=masses[0]
return matdict, coupldict, masses
else:
nmasses = len(data['masses'])
matdictlist=[]
masseslist=[]
def genreplace(mystring,index):
mystring = mystring.replace(f"AUTgencontrib[{index}]","1")
pattern = r"AUTgencontrib\[\d+\]"
mystring = re.sub(pattern, "0", mystring)
return mystring
for genindex in range(1,nmasses+1):
matdictgen = {key: genreplace(value,genindex) for key,value in matdict.items()}
matdictlist.append(matdictgen)
massesgen=masses[genindex-1]
masseslist.append(massesgen)
return matdictlist, coupldict, masseslist
[docs]
def analytical_pipeline(confname):
"""
Function corresponding to the analytical pipeline.
The function goes through the following steps:
#. Validate the configuration
#. Generate two ``.m`` files containing the configuration information for later steps *(if chosen)*
#. Generate the Feynman rules by running ``frules.wl`` *(if chosen)*
#. Generate all processes, compute the total matrix element squared and the thermal masses of mediators by running ``processes.wl`` *(if chosen)*
#. Parse the result into a form that can be used as parameters for :class:`numerical.manipulate.NumRate`
The function can be used as follows:
.. code-block:: python
mat2,masses,couplings = analytical_pipeline('thefile.cfg')
:param str confname: The name of the configuration file. It is expected to be a ``.cfg`` file
:return: Same as ``do_the_dict``
:rtype: dict, dict, tuple
"""
cfg = configparser.ConfigParser()
cfg.read(confname)
# Validate the configuration
# - Check for the existence of each tool
for (tool,path) in cfg.items('Tools'):
if not (os.path.isfile(os.path.expanduser(path)) or os.path.isdir(os.path.expanduser(path))):
print("Required tool \'"+tool+"\' not found. Exiting.")
return
# - Check for the model file
if not os.path.isfile(os.path.expanduser(cfg['Model'].get('modelpath'))):
print("Model file not found. Exiting.")
return
# - Check that the produced particle is correct
if not check_particle(cfg['Model'].get('produced')):
print("Produced particle was set as "+cfg['Model'].get('produced')+" but this option is incorrect. Exiting.")
return
# Empty strings are considered "false"
# (see Andrew Clark and Mateen Ulhaq's answer in https://stackoverflow.com/questions/9573244/how-to-check-if-the-string-is-empty-in-python)
if not cfg['Model'].get('inbath'):
print("Assuming SM thermal bath")
else:
# Check that the particles in the thermal bath are correct
for particle in cfg['Model'].get('inbath').split(' '):
if not check_particle(particle):
print("Particle "+particle+" in the thermal bath is incorrect. Exiting.")
return
# - Check that the AutoTherm directory is correct (using the package configuration file)
if not os.path.isfile(os.path.expanduser(cfg['Run'].get('autothermdir')+'pyproject.toml')):
print("AutoTherm directory is incorrect. Exiting.")
return
# Check the model file for incorrect Lagrangian definitions
try:
check_for_equal(cfg['Model'].get('modelpath'))
except LagrangianError as le:
print(le)
return
# In NumRate only one non-equilibrium coupling is allowed
# Throw a warning for now
if "," in cfg['Model'].get('noneq'):
warnings.warn("Multiple non-equilibrium couplings detected. Execution will continue regardless.",RuntimeWarning)
print("Configuration check completed successfully.")
cdirname= cfg['Run'].get('autothermdir')+"analytical"
# The name of the generated configuration is the same as the configuration file but with a .m extension
conf_name = re.sub("(.*)\\.cfg$",r"\1.m",confname)
# We have to generate a second configuration to avoid problems with shadowing definitions
# The name of that second file is the same as the configuration file but with a "_2" at the end and a .m extension
conf_name_2 = re.sub("(.*)\\.cfg$",r"\1_2.m",confname)
if cfg['Run'].getboolean('conf'):
print("Generating the configuration.")
# Exporting the configuration in a format compatible with Mathematica
# First configuration file
assoc = "<|"
for (tool,path) in cfg.items('Tools'):
if tool == 'math':
continue
assoc+="\""+tool+"\"->"
assoc+="\""+path+"\","
assoc+="\"modelpath\"->\""+cfg['Model'].get('modelpath')+"\","
assoc+="\"lagrangian\"->"+cfg['Model'].get('lagrangian')+","
assoc+="\"autothermdir\"->\""+cfg['Run'].get('autothermdir')+"\","
assoc+="\"verbose\"->"+str(cfg['Run'].getboolean('verbose'))+","
assoc+="\"includeSM\"->"+str(cfg['Model'].getboolean('includeSM'))
assoc+="|>"
# Second configuration file
assoc2 = "<|"
thermal_bath = do_the_thermal_bath(cfg['Model'].get('inbath'))
assoc2+="\"inbath\"->{"
for part in thermal_bath:
assoc2+=part
# If we did the naive solution of adding a comma after each element, then the list will be such as {a,b,c,...,z,}
# This poses a problem to Mathematica which adds a Null element after the last comma
# This is the way I found to avoid this problem for now
# This works because all the elements in the thermal bath are unique
if part == thermal_bath[-1]:
continue
else:
assoc2+=","
assoc2+="},"
assoc2+="\"produced\"->"+cfg['Model'].get('produced')+","
assoc2+="\"assumptions\"->"+"{"+cfg['Model'].get('assumptions')+"},"
assoc2+="\"replacements\"->"+"{"+cfg['Model'].get('replacements')+"},"
assoc2+="\"noneq\"->"+"{"+cfg['Model'].get('noneq')+"},"
assoc2+="\"flavorexpand\"->"+"{"+cfg['Model'].get('flavorexpand')+"}"
assoc2+="|>"
conf_file = open(conf_name,"w")
conf_file.write(assoc)
conf_file.close()
conf_file_2 = open(conf_name_2,"w")
conf_file_2.write(assoc2)
conf_file_2.close()
else:
print("Skipping configuration.")
if (cfg['Run'].getboolean('rules') or cfg['Run'].getboolean('proc')):
warnings.warn("Running the code without generating the configuration first. \n Unexpected effects may arise. \n Execution will continue regardless.",RuntimeWarning)
if cfg['Run'].getboolean('rules'):
print("Generating Feynman rules.")
# According to https://stackoverflow.com/questions/12060863/python-subprocess-call-a-bash-alias (see Jonathan Hartley's answer),
# The run function tries to find an executable, hence why the full path to "math" must be provided
# If I put shell=True, the options passed to "math" are ignored, which is problematic
subprocess.run([cfg['Tools'].get('math'),"-script",cdirname+"/frules.wl",conf_name])
else:
print("Skipping generation of Feynman rules.")
if cfg['Run'].getboolean('proc'):
print("Generating the processes.")
subprocess.run([cfg['Tools'].get('math'),"-script",cdirname+"/processes.wl",conf_name])
else:
print("Skipping generation of all processes.")
# Run the function for parsing the result into the form used by the numerical part
# Move it here for better testing, don't know if it should be run only if the processes have been generated
result_filename = confname[:-4]+'_result.json'
try:
the_mat,the_coupling,the_mass = do_the_dict(result_filename)
# Handle the file not existing / not having been created properly
except RuntimeError:
print("The file wasn't found")
# Handle the file being formatted incorrectly
except json.JSONDecodeError:
print("The JSON formatting is incorrect")
# Handle if the particle natures are incorrect
except AssertionError:
print("Something went wrong.")
print("Analytical pipeline completed")
return the_mat, the_coupling, the_mass
def run():
try:
filename = sys.argv[1]
except IndexError:
print("No configuration file provided. Exiting.")
sys.exit(1)
if not os.path.isfile(filename):
print("Configuration file not found. Exiting.")
sys.exit(1)
a,b,c = analytical_pipeline(filename)
print("The matrix elements squared read:")
if type(a) == list:
print("[",end="")
listlen=len(a)
alist=a
else:
listlen=1
alist=[a]
j=0
endl=","
for a in alist:
if j == listlen-1:
if listlen>1:
endl="]"
else:
endl=""
print("{",end="")
alen=len(a)
i=0
endc=","
for key, value in a.items():
if i ==alen-1:
endc="}"+endl
print(f"{key}: '{value}'{endc}")
i+=1
j += 1
print("The couplings are:")
print(b)
print("The thermal masses are:")
print(c)
# Because we want to use the code as a library (from a notebook for example), it's better if the analytical pipeline is a function
if __name__ == "__main__":
run()
