-----------------
Running AutoTherm
-----------------

.. analytical:

The code can be run entirely within a ``jupyter`` notebook or through 
shell commands and python scripts. Here we discuss extensively the latter option.
An :doc:`example notebook<../examples/quickstart>` is also provided. To run it,
you need to install jupyter within the virtual environment wherein AutoTherm has 
been installed, or vice-versa.


The analytical, Mathematica-based part
--------------------------------------

We start by activating the virtual environment, if not already active

.. code-block:: bash
    
    cd $AutothermDir
    source .venv/bin/activate

Here ``$AutothermDir`` should be replaced by its actual value, as discussed in the :doc:`installation section<install>`.

We then move to the ``$ModelDir`` directory where the model file and config file have been created in the :doc:`previous section<model>` and 
we launch the ``python`` executable handling the creation of configuration files and the execution of the Wolfram language 
scripts. This is done through

.. code-block:: bash
    
    cd $ModelDir
    autotherm-run axionSM.cfg

Barring any bugs with the model file or with the code, this should go through all steps. The final lines of the output should then read

.. code-block::

    Analytical pipeline completed
    The matrix elements squared read:
    {(1, -1, -1): '(24*ct^2*ht^2*s)/fPQ^2 + ((5*c1^2*g1^6 + 9*c2^2*g2^6 + 24*c3^2*g3^6)*(t^2 + u^2))/(32*fPQ^2*pi^4*s)', (-1, -1, 1): '-1/32*(768*ct^2*ht^2*pi^4*t^2 + 5*c1^2*g1^6*(t^2 + 2*t*u + 2*u^2) + 9*c2^2*g2^6*(t^2 + 2*t*u + 2*u^2) + 24*c3^2*g3^6*(t^2 + 2*t*u + 2*u^2))/(fPQ^2*pi^4*t)', (-1, 1, -1): '-1/32*(768*ct^2*ht^2*pi^4*u^2 + 5*c1^2*g1^6*(2*t^2 + 2*t*u + u^2) + 9*c2^2*g2^6*(2*t^2 + 2*t*u + u^2) + 24*c3^2*g3^6*(2*t^2 + 2*t*u + u^2))/(fPQ^2*pi^4*u)', (1, 1, 1): '((c1^2*g1^6 + 15*c2^2*g2^6 + 48*c3^2*g3^6)*(t^2 + t*u + u^2)^2)/(64*fPQ^2*pi^4*s*t*u)'}
    The couplings are:
    {'gauge': ('g1', 'g2', 'g3'), 'noneq': ('fPQ',), 'others': ('c1', 'c2', 'c3', 'ct', 'ht')}
    The thermal masses are:
    ('(11*g1^2)/6', '(11*g2^2)/6', '2*g3^2')

These are the main results of this step: the matrix elements squared and the thermal masses that are needed. 
They have been written in the ``axionSM_result.json`` file within ``$ModelDir``. On an M3 MacBook Pro with Python 3.12 
and Mathematica 14 this takes about 20 seconds.

.. numerical:

The numerical part
------------------

The final step is to perform the phase-space integration over these matrix elements. This is carried out by the ``numerical``
module. To this end, we can create this simple python script within ``$ModelDir``, which we call ``get_rate.py``.

.. code-block:: python

    #start by importing the analytical and numerical modules
    from analytical.controller import *
    from numerical.manipulate import *


    #this function calls the analytical pipeline. To avoid re-running all the steps,
    #set the conf, rules and proc flags in axionSM.cfg to no
    axion_analytical=analytical_pipeline("axionSM.cfg")

    #we now create an instance of the NumRate class using the output of the 
    #analytical pipeline. 1 is the degeneracy of the axion (no spin or any other quantum number)
    axion_rate=NumRate(*axion_analytical,1)

    #now all methods of the class are available. For illustration purposes we output 
    #the full rate at 10 values of k/T for an axion-gluon coupling only. We rescale 
    #fpq by setting its value to 1 and we use a strong coupling g_3 corresponding to
    #alpha_s=0.06
    k=numpy.linspace(1,10,10)
    #the tuple for the couplings is the union of the gauge and other couplings
    #from the output of the analytical pipeline we have 
    # {'gauge': ('g1', 'g2', 'g3'), 'noneq': ('fPQ',), 'others': ('c1', 'c2', 'c3', 'ct', 'ht')}
    rate = axion_rate.rate(k,1.,(0.,0.,numpy.sqrt(4.*numpy.pi*0.06),0,0,1.,0,0),0)
    print(rate)

Before executing this file, it is recommended to set the ``conf``, ``rules`` and ``proc`` flags in ``axionSM.cfg`` to ``no``, so 
as to avoid re-running the Wolfram language scripts. 
The script is executed with

.. code-block:: shell

    python get_rate.py

Barring bugs, it should output

.. code-block::

    Assuming SM thermal bath
    Configuration check completed successfully.
    Skipping configuration.
    Skipping generation of Feynman rules.
    Skipping generation of all processes.
    Analytical pipeline completed
    /Users/jacopo/Nextcloud/AUTOTHERM/autotherm/MyModels/axion/get_rate.py:12: AutothermWarning: The ratio between the gauge boson mass 11*g1**2/6 and the coefficient of the IR divergence -11*c1**2*g1**6/(384*pi**4) depends on the coupling constants entering the mass.
    This is most likely not problematic, continuing with evaluation.
    axion_rate=NumRate(*axion_analytical,1)
    /Users/jacopo/Nextcloud/AUTOTHERM/autotherm/MyModels/axion/get_rate.py:12: AutothermWarning: The ratio between the gauge boson mass 11*g2**2/6 and the coefficient of the IR divergence -11*c2**2*g2**6/(128*pi**4) depends on the coupling constants entering the mass.
    This is most likely not problematic, continuing with evaluation.
    axion_rate=NumRate(*axion_analytical,1)
    /Users/jacopo/Nextcloud/AUTOTHERM/autotherm/MyModels/axion/get_rate.py:12: AutothermWarning: The ratio between the gauge boson mass 2*g3**2 and the coefficient of the IR divergence -c3**2*g3**6/(4*pi**4) depends on the coupling constants entering the mass.
    This is most likely not problematic, continuing with evaluation.
    axion_rate=NumRate(*axion_analytical,1)
    [array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.]), array([6.24468720482499e-6, 2.48806377023693e-5, 3.35216329472938e-5,
        3.89866431064388e-5, 4.30820774272845e-5, 4.64391827454089e-5,
        4.93197028763330e-5, 5.18530080251627e-5, 5.41151703842284e-5,
        5.61573931901464e-5], dtype=object), array([1.32458357276370e-5, 2.68517915080078e-5, 3.44196852859558e-5,
        3.94963238927934e-5, 4.34096375330550e-5, 4.66671741757908e-5,
        4.94874378831309e-5, 5.19815453260401e-5, 5.42167932436927e-5,
        5.62397439824078e-5], dtype=object), array([2.13516715648551e-5, 3.50716746267359e-5, 4.25672158519041e-5,
        4.74800571298319e-5, 5.12012593583113e-5, 5.42733060334197e-5,
        5.69285281412678e-5, 5.92808393892183e-5, 6.13957319484079e-5,
        6.33165674566682e-5], dtype=object)]

The warnings about proportionality are harmless. The numerical integration in the three schemes has been successfully carried out. On 
an M3 MacBook Pro with Python 3.12 this takes about 3 seconds.