ppllnnjj_ew.cpp

#define GNU_NAME_MANGLING
#include <recola.h>

#include <array>
#include <cassert>
#include <iostream>
#include <iomanip>
#include <vector>

extern "C"
{

    void ol_setparameter_int(const char* param, int val);
    void ol_setparameter_double(const char* param, double val);
    int ol_register_process(const char* process, int amptype);
    void ol_start();
    void ol_finish();
    void ol_evaluate_tree(int id, double* pp, double* m2tree);
    void ol_evaluate_cc(int id, double* pp, double* m2tree, double* m2cc, double* m2ew);
    void ol_evaluate_loop(int id, double* pp, double* m2tree, double* m2loop, double* acc);
}

enum conv_t {
    blha = 0,
    coli = 1
};

double openloops(int alphas, conv_t conv) {
    // TODO: additional parameters not set for Recola
    ol_setparameter_int("ew_scheme", 1);
    ol_setparameter_int("ew_renorm_scheme", 1);
    ol_setparameter_int("ew_renorm", 1);
    ol_setparameter_int("verbose", 2);

    ol_setparameter_double("mass(23)", 91.153480619183);
    ol_setparameter_double("width(23)", 2.4942663787728);
    ol_setparameter_double("mass(24)", 80.351984549666);
    ol_setparameter_double("width(24)", 2.0837993977753);
    ol_setparameter_double("mass(25)", 125.0);
    ol_setparameter_double("width(25)", 4.07e-3);
    ol_setparameter_double("mass(6)", 173.0);
    ol_setparameter_double("width(6)", 0.0);
    ol_setparameter_double("gmu", 1.1663787e-5);
    ol_setparameter_int("complex_mass_scheme", 2);

    ol_setparameter_double("alphas", 1.0);
    ol_setparameter_double("mureg", 100.0);
    ol_setparameter_double("muren", 100.0);

    switch (conv)
    {
    case blha:
        ol_setparameter_int("polenorm", 0);
        break;

    case coli:
        //ol_setparameter_int("polenorm", 0);
        ol_setparameter_int("polenorm", 1);
        break;

    default:
        assert(false);
    }

    // color-correlated matrix elements are in the library containing the virtuals
    ol_setparameter_int("loop_order_qcd", alphas);
    int id = ol_register_process("2 2 -> -11 12 -13 14 1 1", 11);

    //ol_setparameter_int("order_qcd", alphas);
    //int id = ol_register_process("2 2 -> -11 12 -13 14 1 1", 1);

    //ol_setparameter_int("write_parameters", 1);
    //ol_setparameter_int("verbose", 1);

    ol_start();

    double mom[] = {
        2474.392899476,    0.000000000,    0.000000000, -2474.392899476, 0.0,
        2474.392899476,    0.000000000,    0.000000000,  2474.392899476, 0.0,
         256.283676529,  202.559375939,   -5.247813058,   156.918713137, 0.0,
         372.294941417,  254.501327544,  -95.429348283,   254.412729968, 0.0,
         962.071502755, -592.225093683, -719.648513817,  -238.656722905, 0.0,
         147.178745133,  -63.691489650, -127.366643437,   -37.184880055, 0.0,
        1508.209608516,  -89.033013685,  -21.715829078, -1505.422787280, 0.0,
        1702.747324602,  287.888893536,  969.408147672,  1369.932947134, 0.0
    };

    double result;
    //ol_evaluate_tree(id, mom, &result);

    double m2tree;
    //std::vector<double> m2cc(8 * (8 - 1) / 2);
    //double m2ew;
    //ol_evaluate_cc(id, mom, &m2tree, m2cc.data(), &m2ew);
    //double result = m2cc.at(i + j * (j - 1) / 2);
    double m2loop[3];
    double acc;
    ol_evaluate_loop(id, mom, &m2tree, m2loop, &acc);

    std::cout << "m[0]: " << m2loop[0] << '\n';
    std::cout << "m[1]: " << m2loop[1] << '\n';
    std::cout << "m[2]: " << m2loop[2] << '\n';
    std::cout << "mtree:" << m2tree << '\n';

    result = m2loop[0];

    //if (conv == coli)
    //{
    //    double pi = std::acos(-1.0);
    //    result -= pi * pi / 6.0 * m2loop[2];
    //}

    ol_finish();

    return result;
}

double recola(int alphas, conv_t conv) {
    Recola::set_complex_mass_scheme_rcl();
    Recola::set_pole_mass_z_rcl(91.153480619183, 2.4942663787728);
    Recola::set_pole_mass_w_rcl(80.351984549666, 2.0837993977753);
    Recola::set_pole_mass_h_rcl(125.0, 4.07e-3);
    Recola::set_pole_mass_top_rcl(173.0, 0.0);
    Recola::use_gfermi_scheme_and_set_gfermi_rcl(1.1663787e-5);
    Recola::set_print_level_squared_amplitude_rcl(1);

    switch (conv)
    {
    case blha:
        {
            double pi = std::acos(-1.0);
            Recola::set_delta_ir_rcl(0.0, pi * pi / 6.0);
        }
        break;

    case coli:
        break;

    default:
        assert(false);
    }

    // this isn't Q=100 GeV, but we don't vary the scale so it shouldn't matter
    Recola::set_alphas_rcl(1.0, 100.0, 5);

    Recola::define_process_rcl(1, "u u -> e+ nu_e mu+ nu_mu d d", "NLO");
    Recola::unselect_all_gs_powers_LoopAmpl_rcl(1);
    Recola::unselect_all_gs_powers_BornAmpl_rcl(1);

    std::array<int, 2> born_gs = { 0, 2 };
    std::array<int, 3> loop_gs = { 0, 2, 4 };

    for (int born : born_gs) {
        for (int loop : loop_gs) {
            if ((born + loop) == 2 * alphas) {
                Recola::select_gs_power_LoopAmpl_rcl(1, loop);
                Recola::select_gs_power_BornAmpl_rcl(1, born);
            }
        }
    }

    Recola::generate_processes_rcl();

    double mom[][4] = {
        { 2474.392899476,    0.000000000,    0.000000000, -2474.392899476 },
        { 2474.392899476,    0.000000000,    0.000000000,  2474.392899476 },
        {  256.283676529,  202.559375939,   -5.247813058,   156.918713137 },
        {  372.294941417,  254.501327544,  -95.429348283,   254.412729968 },
        {  962.071502755, -592.225093683, -719.648513817,  -238.656722905 },
        {  147.178745133,  -63.691489650, -127.366643437,   -37.184880055 },
        { 1508.209608516,  -89.033013685,  -21.715829078, -1505.422787280 },
        { 1702.747324602,  287.888893536,  969.408147672,  1369.932947134 },
    };

    //double result = 0.0;
    //Recola::set_print_level_squared_amplitude_rcl(1);
    //Recola::compute_process_rcl(1, mom, "LO");
    //Recola::get_squared_amplitude_rcl(1, alphas, "LO", result);

    //double result = 0.0;
    //Recola::set_print_level_correlations_rcl(1);
    //Recola::compute_colour_correlation_rcl(1, mom, i + 1, j + 1);
    //Recola::get_colour_correlation_rcl(1, alphas - 1, i + 1, j + 1, result);
    //result *= 4.0 / 3.0;
    double result[2] = {};
    Recola::compute_process_rcl(1, mom, "NLO", result);
    Recola::get_squared_amplitude_rcl(1, alphas, "NLO", result[1]);
    Recola::get_squared_amplitude_rcl(1, 0, "LO", result[0]);

    //std::cout << "a(0): " << result[0] << '\n';
    //std::cout << "a(1): " << result[1] << '\n';
    //std::cout << "a   : " << result[2] << '\n';

    return result[1];
}

int main(int argc, char* argv[]) {
    int alphas = std::stoi(argv[1]);
    conv_t conv = static_cast <conv_t> (std::stoi(argv[2]));

    auto recola = ::recola(alphas, conv);
    auto openloops = ::openloops(alphas, conv);

    std::cout << "O(as^" << alphas << ") ";
    switch (conv) {
    case blha:
        std::cout << "(BLHA): \n";
        break;

    case coli:
        std::cout << "(COLI): \n";
        break;

    default:
        assert(false);
    }

    std::cout << "Recola:    " << std::setprecision(13) << recola << '\n';
    std::cout << "OpenLoops: " << std::setprecision(13) << openloops << '\n';
}