file frontends/SUSY_HIT_1_5.cpp

[No description available] More…

Detailed Description

Author: Pat Scott

Date: 2015 Jan-May

Frontend for SUSY-HIT 1.5 backend


Authors (add name and date if you modify):


Source code

//   GAMBIT: Global and Modular BSM Inference Tool
//   *********************************************
///  \file
///
///  Frontend for SUSY-HIT 1.5 backend
///
///  *********************************************
///
///  Authors (add name and date if you modify):
///
/// \author Pat Scott
/// \date 2015 Jan-May
///
///  *********************************************

#include <sstream>

#include "gambit/Backends/frontend_macros.hpp"
#include "gambit/Backends/frontends/SUSY_HIT_1_5.hpp"
#include "gambit/Elements/mssm_slhahelp.hpp"
#include "gambit/Utils/slhaea_helpers.hpp"


// Convenience functions (definitions)
BE_NAMESPACE
{

  /// Shortcut for dealing with SLHA blocks
  void required_block(const str& name, SLHAea::Block& block, const SLHAea::Coll& slha)
  {
    if (slha.find(name) != slha.end()) block = slha.at(name);
    else backend_error().raise(LOCAL_INFO, "Sorry, SUSY-HIT needs SLHA block: " + name + ".\n"
    "If you tried to read in a debug SLHA file with missing entries (e.g. STOPMIX, STAUMIX,\n"
    "etc), then sort out your SLHA file so that it is readable by SUSY-HIT!  If you can't  \n"
    "do that, then depending on what blocks your file has, you may be able to instead use  \n"
    "it to initialise either                                                               \n"
    " 1. a GAMBIT Spectrum object (see SpecBit::get_MSSM_spectrum_from_SLHAfile), or       \n"
    " 2. a GAMBIT DecayTable object (see DecayBit::all_decays_from_SLHA).                  \n");
  }
  /// @}

  /// Runs actual SUSY-HIT decay calculations.
  /// Inputs: m_s_1GeV_msbar    strange mass in GeV, in MSbar scheme at an energy of 1GeV
  ///         W_width, Z_width  EW gauge boson total widths in GeV
  void run_susy_hit(SLHAstruct slha, double W_width, double Z_width)
  {
    using SLHAea::to;

    // Set out all the SLHA1 PDG codes
    int pdg_codes[35];
    pdg_codes[0]  = 24;      // W+
    pdg_codes[1]  = 25;      // h
    pdg_codes[2]  = 35;      // H
    pdg_codes[3]  = 36;      // A
    pdg_codes[4]  = 37;      // H+
    pdg_codes[5]  = 1000001; // ~d_L
    pdg_codes[6]  = 2000001; // ~d_R
    pdg_codes[7]  = 1000002; // ~u_L
    pdg_codes[8]  = 2000002; // ~u_R
    pdg_codes[9]  = 1000003; // ~s_L
    pdg_codes[10] = 2000003; // ~s_R
    pdg_codes[11] = 1000004; // ~c_L
    pdg_codes[12] = 2000004; // ~c_R
    pdg_codes[13] = 1000005; // ~b_1
    pdg_codes[14] = 2000005; // ~b_2
    pdg_codes[15] = 1000006; // ~t_1
    pdg_codes[16] = 2000006; // ~t_2
    pdg_codes[17] = 1000011; // ~e_L
    pdg_codes[18] = 2000011; // ~e_R
    pdg_codes[19] = 1000012; // ~nu_eL
    pdg_codes[20] = 1000013; // ~mu_L
    pdg_codes[21] = 2000013; // ~mu_R
    pdg_codes[22] = 1000014; // ~nu_muL
    pdg_codes[23] = 1000015; // ~tau_1
    pdg_codes[24] = 2000015; // ~tau_2
    pdg_codes[25] = 1000016; // ~nu_tauL
    pdg_codes[26] = 1000021; // ~g
    pdg_codes[27] = 1000022; // ~chi_10
    pdg_codes[28] = 1000023; // ~chi_20
    pdg_codes[29] = 1000025; // ~chi_30
    pdg_codes[30] = 1000035; // ~chi_40
    pdg_codes[31] = 1000024; // ~chi_1+
    pdg_codes[32] = 1000037; // ~chi_2+
    pdg_codes[33] = 5;       // b pole
    pdg_codes[34] = 1000039; // ~G

    // Zero all SLHA2 Q values
    for (int i=1; i<=22; ++i) sd_leshouches2->qvalue(i) = 0.0;

    // Get SLHA2 blocks
    SLHAea::Block sminputs, vckm, msoft, mass, nmix, vmix, umix;
    SLHAea::Block stopmix, sbotmix, staumix, alpha, hmix, gauge;
    SLHAea::Block au, ad, ae, yu, yd, ye, msq2, msd2, msu2, td, tu;
    SLHAea::Block usqmix, dsqmix, selmix, spinfo;
    required_block("SPINFO",   spinfo,   slha);
    required_block("SMINPUTS", sminputs, slha);
    required_block("VCKMIN",   vckm,     slha);
    required_block("MSOFT",    msoft,    slha);
    required_block("MASS",     mass,     slha);
    required_block("NMIX",     nmix,     slha);
    required_block("VMIX",     vmix,     slha);
    required_block("UMIX",     umix,     slha);
    required_block("ALPHA",    alpha,    slha);
    required_block("HMIX",     hmix,     slha);
    required_block("GAUGE",    gauge,    slha);
    required_block("STOPMIX",  stopmix,  slha);
    required_block("SBOTMIX",  sbotmix,  slha);
    required_block("STAUMIX",  staumix,  slha);
    required_block("YU",       yu,       slha);
    required_block("YD",       yd,       slha);
    required_block("YE",       ye,       slha);
    required_block("AU",       au,       slha);
    required_block("AD",       ad,       slha);
    required_block("AE",       ae,       slha);

    // SPINFO
    sd_leshouches1->spinfo1 = (spinfo[1].is_data_line()) ? spinfo[1][1] : ""; // RGE +Spectrum calculator
    sd_leshouches1->spinfo2 = (spinfo[2].is_data_line()) ? spinfo[2][1] : ""; // version number

    // SMINPUTS
    for (int i=1; i<=14; ++i)
    {
      sd_leshouches2->smval(i) = (sminputs[i].is_data_line()) ? to<double>(sminputs[i][1]) : 0.0;
    }
    for (int i=15; i<=20; ++i) sd_leshouches2->smval(i) = 0.0;     // zeroing

    // SUSY-HIT and HDecay non-SLHA inputs
    susyhitin->amsin = to<double>(sminputs.at(23).at(1));          // MSBAR(1): HDECAY claims it wants ms(1GeV)^MSbar, but we don't believe it, and give it m_s(2GeV)^MSBar
    susyhitin->amcin = to<double>(sminputs.at(24).at(1));          // MC: HDECAY claims it wants the c pole mass, but that is not well defined, so we give it mc(mc)^MSBar
    susyhitin->ammuonin = sd_leshouches2->smval(13);               // MMUON: mmu(pole)
    susyhitin->alphin = sd_leshouches2->smval(1);                  // ALPHA: alpha_em^-1(M_Z)^MSbar (scheme and scale not specified in SUSYHIT or HDECAY documentation)
    susyhitin->gamwin = W_width;                                   // GAMW: W width (GeV)
    susyhitin->gamzin = Z_width;                                   // GAMZ: Z width (GeV)
    double lambda = to<double>(vckm.at(1).at(1));                  // Wolfenstein lambda
    double A = to<double>(vckm.at(2).at(1));                       // Wolfenstein A
    double rhobar = to<double>(vckm.at(3).at(1));                  // Wolfenstein rhobar
    double etabar = to<double>(vckm.at(4).at(1));                  // Wolfenstein etabar
    susyhitin->vusin = Spectrum::Wolf2V_us(lambda, A, rhobar, etabar); // VUS: |CKM V_us|
    susyhitin->vcbin = Spectrum::Wolf2V_cb(lambda, A, rhobar, etabar); // VCB: |CKM V_cb|
    susyhitin->rvubin = std::abs(Spectrum::Wolf2V_ub(lambda, A, rhobar, etabar)) / susyhitin->vcbin; // VUB/VCB: Ratio of CKM |V_UB|/|V_CB|

    // MODSEL
    sd_leshouches2->imod(1) = 1;    // 1, x = SUSY model info.
    sd_leshouches2->imod(2) = 0;    // 0 = general MSSM, 1 = SUGRA => do not calculate metastable NLSP decays to gravitinos.
    //sd_leshouches2->imod(2) = 2;  // 2 = GMSB => calculate metastable NLSP decays to gravitinos.  Not yet implemented in GAMBIT.
    checkfavvio->imodfav(1) = 6;    // 6, x =  Flavour violation info.
    checkfavvio->imodfav(2) = 0;    // 0 = no FV, 1 = q FV, 2 = l FV, 3 = both.  For FV decays, must be set non-zero later before calling sdecay() again.  Not yet implemented in GAMBIT.

    // MSOFT
    for (int i=1; i<=100; ++i) sd_leshouches2->msoftval(i) = unlikely(); // indicate undefined
    if (msoft.find_block_def()->size() >= 4) sd_leshouches2->qvalue(3) = to<double>(msoft.find_block_def()->at(3));// Q(GeV)
    int msoft_indices[20] = {1, 2, 3, 21, 22, 31, 32, 33, 34, 35, 36, 41, 42, 43, 44, 45, 46, 47, 48, 49};
    // M_1, M_2, M_3, M^2_Hd, M^2_Hu, M_eL, M_muL, M_tauL, M_eR, M_muR, M_tauR, M_q1L, M_q2L, M_q3L, M_uR, M_cR, M_tR, M_dR, M_sR, M_bR
    for (int i : msoft_indices)
    {
      if (msoft[i].is_data_line()) sd_leshouches2->msoftval(i) = to<double>(msoft.at(i).at(1));
      //std::cout << "msoft[" << i << "] = " << sd_leshouches2->msoftval(i) << std::endl;
    }

    // EXTPAR
    sd_leshouches2->extval(0) = sd_leshouches2->qvalue(3);         // EWSB scale (set to SUSY scale as per MSOFT).  Not used by SUSY-HIT anymore.
    //std::cout << "extval = " << sd_leshouches2->extval(0) << std::endl;

    // MASS
    for (int i=1; i<=35; ++i)
    {
      sd_leshouches2->massval(i) = (mass[pdg_codes[i-1]].is_data_line()) ? to<double>(mass[pdg_codes[i-1]][1]) : unlikely();
      //std::cout << "massval(" << pdg_codes[i-1] << ") = " << sd_leshouches2->massval(i) << std::endl;
    }
    for (int i=36; i<=50; ++i) sd_leshouches2->massval(i) = 0.0; // zeroing

    // NMIX
    for (int i=1; i<=4; ++i) for (int j=1; j<=4; ++j) sd_leshouches2->nmixval(i,j) = (nmix[initVector<int>(i,j)].is_data_line()) ? to<double>(nmix.at(i,j)[2]) : 0.0;
    //for (int i=1; i<=4; ++i) for (int j=1; j<=4; ++j) std::cout << "nmixval(" << i << "," << j << ") = " << sd_leshouches2->nmixval(i,j) << std::endl;


    // VMIX, UMIX, STOPMIX, SBOTMIX, STAUMIX
    for (int i=1; i<=2; ++i)
    {
      for (int j=1; j<=2; ++j)
      {
        std::vector<int> ij = initVector<int>(i,j);
        sd_leshouches2->vmixval(i,j) = (vmix[ij].is_data_line()) ? to<double>(vmix.at(i,j)[2]) : 0.0;
        sd_leshouches2->umixval(i,j) = (umix[ij].is_data_line()) ? to<double>(umix.at(i,j)[2]) : 0.0;
        sd_leshouches2->stopmixval(i,j) = (stopmix[ij].is_data_line()) ? to<double>(stopmix.at(i,j)[2]) : 0.0;
        sd_leshouches2->sbotmixval(i,j) = (sbotmix[ij].is_data_line()) ? to<double>(sbotmix.at(i,j)[2]) : 0.0;
        sd_leshouches2->staumixval(i,j) = (staumix[ij].is_data_line()) ? to<double>(staumix.at(i,j)[2]) : 0.0;
        //std::cout << "vmixval(" << i << "," << j << ") = " << sd_leshouches2->vmixval(i,j) << std::endl;
        //std::cout << "umixval(" << i << "," << j << ") = " << sd_leshouches2->umixval(i,j) << std::endl;
        //std::cout << "stopmixval(" << i << "," << j << ") = " << sd_leshouches2->stopmixval(i,j) << std::endl;
        //std::cout << "sbotmixval(" << i << "," << j << ") = " << sd_leshouches2->sbotmixval(i,j) << std::endl;
        //std::cout << "staumixval(" << i << "," << j << ") = " << sd_leshouches2->staumixval(i,j) << std::endl;
      }
    }

    // ALPHA (value is spectrum generator's "best choice" => can be on-shell, DRbar at a given scale, whatever)
    sd_leshouches2->alphaval = to<double>(alpha.back().at(0));             // Mixing angle in the neutral Higgs boson sector.
    //std::cout << "alphaval = " << sd_leshouches2->alphaval << std::endl;

    // HMIX
    if (hmix.find_block_def()->size()  >= 4) sd_leshouches2->qvalue(1) = to<double>(hmix.find_block_def()->at(3)); // Q(GeV)
    for (int i=1; i<=10; ++i) sd_leshouches2->hmixval(i) = (hmix[i].is_data_line()) ? to<double>(hmix[i][1]) : unlikely();
    //for (int i=1; i<=10; ++i) std::cout << "hmixval(" << i << ") = " << sd_leshouches2->hmixval(i) << std::endl;

    // GAUGE
    if (gauge.find_block_def()->size() >= 4) sd_leshouches2->qvalue(2) = to<double>(gauge.find_block_def()->at(3));// Q(GeV)
    sd_leshouches2->gaugeval(1) = to<double>(gauge.at(1).at(1));                                                   // gprime(Q) DRbar
    sd_leshouches2->gaugeval(2) = to<double>(gauge.at(2).at(1));                                                   // g(Q) DRbar
    sd_leshouches2->gaugeval(3) = to<double>(gauge.at(3).at(1));                                                   // g_3(Q) DRbar
    //for(int i=1;i<=3;i++) std::cout << "gaugeval(" << i << ") = " << sd_leshouches2->gaugeval(i) << std::endl;

    // AU, AD, AE, YU, YD, YE
    if (au.find_block_def()->size() >= 4) sd_leshouches2->qvalue(4)  = to<double>(au.find_block_def()->at(3));     // Q(GeV)
    if (ad.find_block_def()->size() >= 4) sd_leshouches2->qvalue(5)  = to<double>(ad.find_block_def()->at(3));     // Q(GeV)
    if (ae.find_block_def()->size() >= 4) sd_leshouches2->qvalue(6)  = to<double>(ae.find_block_def()->at(3));     // Q(GeV)
    if (yu.find_block_def()->size() >= 4) sd_leshouches2->qvalue(7)  = to<double>(yu.find_block_def()->at(3));     // Q(GeV)
    if (yd.find_block_def()->size() >= 4) sd_leshouches2->qvalue(8)  = to<double>(yd.find_block_def()->at(3));     // Q(GeV)
    if (ye.find_block_def()->size() >= 4) sd_leshouches2->qvalue(9)  = to<double>(ye.find_block_def()->at(3));     // Q(GeV)
    for (int i=1; i<=3; ++i)
    {
      for (int j=1; j<=3; ++j)
      {
        std::vector<int> ij = initVector<int>(i,j);
        sd_leshouches2->auval(i,j) = (au[ij].is_data_line()) ? to<double>(au.at(i,j)[2]) : unlikely();
        sd_leshouches2->adval(i,j) = (ad[ij].is_data_line()) ? to<double>(ad.at(i,j)[2]) : unlikely();
        sd_leshouches2->aeval(i,j) = (ae[ij].is_data_line()) ? to<double>(ae.at(i,j)[2]) : unlikely();
        sd_leshouches2->yuval(i,j) = (yu[ij].is_data_line()) ? to<double>(yu.at(i,j)[2]) : 0.0;
        sd_leshouches2->ydval(i,j) = (yd[ij].is_data_line()) ? to<double>(yd.at(i,j)[2]) : 0.0;
        sd_leshouches2->yeval(i,j) = (ye[ij].is_data_line()) ? to<double>(ye.at(i,j)[2]) : 0.0;
      }
    }

    // SLHA2 blocks zeroed for safety.
    for (int i=1; i<=6; ++i)
    {
      for (int j=1; j<=6; ++j)
      {
        if (i <=3 and j<=3)
        {
          flavviolation->msq2(i,j) = 0.0;
          flavviolation->msd2(i,j) = 0.0;
          flavviolation->msu2(i,j) = 0.0;
          flavviolation->td(i,j) = 0.0;
          flavviolation->tu(i,j) = 0.0;
          flavviolation->vckm(i,j) = 0.0;
        }
        flavviolation->usqmix(i,j) = 0.0;
        flavviolation->dsqmix(i,j) = 0.0;
        sd_selectron->selmix(i,j)  = 0.0;
      }
    }

    // Tell SUSY-HIT not to bother calculating the b pole mass from mb(mb)_MSbar, just use the value we pass it.
    sd_mbmb->i_sd_mbmb = 1;

    // Do calculation without flavour-violating light stop decays.
    flavviolation->ifavvio = 0;

    // The following is only relevant if considering FV stop decays and SLHA2 inputs.  Code below is tested in as much as it sets
    // the right common blocks with the right stuff.  We haven't tested what sorts of results you can get out of it though.
    // In order to turn FV on, some time after the initial call to sdecay() at the end of this function, you
    // need to organise a method to read in the SLHA2 version of the spectrum object, then set
    //  flavviolation->ifavvio = 1;
    //  checkfavvio->imodfav(1) = 6;
    //  checkfavvio->imodfav(2) = 1;
    // and then call sdecay() again.  You might do that somehow in this function, or in a module function.
    // There may be a smart way to order this so that it happens automatically if and when you want it.
    // However, before spending time automating this, it still needs to be tested that running first without
    // flavour violation and then re-running with flavour violation actually works, i.e. it does not break
    // the non-FV results.  The SUSY-HIT authors say they think it should be OK, but they have never done it.
    /*

    // MSQ2, MSD2, MSU2, TD, TU
    if (msq2.find_block_def()->size() >= 4) sd_leshouches2->qvalue(17) = to<double>(msq2.find_block_def()->at(3)); // Q(GeV) corrects minor bug in SUSY-HIT SLHA reader
    if (msd2.find_block_def()->size() >= 4) sd_leshouches2->qvalue(10) = to<double>(msd2.find_block_def()->at(3)); // Q(GeV)
    if (msu2.find_block_def()->size() >= 4) sd_leshouches2->qvalue(11) = to<double>(msu2.find_block_def()->at(3)); // Q(GeV)
    if (td.find_block_def()->size() >= 4) sd_leshouches2->qvalue(12) = to<double>(td.find_block_def()->at(3));     // Q(GeV)
    if (tu.find_block_def()->size() >= 4) sd_leshouches2->qvalue(13) = to<double>(tu.find_block_def()->at(3));     // Q(GeV)
    for (int i=1; i<=3; ++i)
    {
      for (int j=1; j<=3; ++j)
      {
        std::vector<int> ij = initVector<int>(i,j);
        flavviolation->msq2(i,j) = (msq2[ij].is_data_line()) ? to<double>(msq2.at(i,j)[2]) : 0.0;
        flavviolation->msd2(i,j) = (msd2[ij].is_data_line()) ? to<double>(msd2.at(i,j)[2]) : 0.0;
        flavviolation->msu2(i,j) = (msu2[ij].is_data_line()) ? to<double>(msu2.at(i,j)[2]) : 0.0;
        flavviolation->td(i,j) = (td[ij].is_data_line()) ? to<double>(td.at(i,j)[2]) : 0.0;
        flavviolation->tu(i,j) = (tu[ij].is_data_line()) ? to<double>(tu.at(i,j)[2]) : 0.0;
      }
    }

    //VCKM
    if (vckm.find_block_def()->size() >= 4) sd_leshouches2->qvalue(21) = to<double>(vckm.find_block_def()->at(3)); // Q(GeV)
    flavviolation->vckm(1,1) =          Spectrum::Wolf2V_ud(lambda, A, rhobar, etabar);
    flavviolation->vckm(1,2) =          Spectrum::Wolf2V_us(lambda, A, rhobar, etabar);
    flavviolation->vckm(1,3) = std::abs(Spectrum::Wolf2V_ud(lambda, A, rhobar, etabar));
    flavviolation->vckm(2,1) = std::abs(Spectrum::Wolf2V_cb(lambda, A, rhobar, etabar));
    flavviolation->vckm(2,2) = std::abs(Spectrum::Wolf2V_cs(lambda, A, rhobar, etabar));
    flavviolation->vckm(2,3) =          Spectrum::Wolf2V_cb(lambda, A, rhobar, etabar);
    flavviolation->vckm(3,1) = std::abs(Spectrum::Wolf2V_td(lambda, A, rhobar, etabar));
    flavviolation->vckm(3,2) = std::abs(Spectrum::Wolf2V_ts(lambda, A, rhobar, etabar));
    flavviolation->vckm(3,3) =          Spectrum::Wolf2V_tb(lambda, A, rhobar, etabar);

    // USQMIX, DSQMIX, SELMIX
    if (usqmix.find_block_def()->size() >= 4) sd_leshouches2->qvalue(14) = to<double>(usqmix.find_block_def()->at(3)); // Q(GeV)
    if (dsqmix.find_block_def()->size() >= 4) sd_leshouches2->qvalue(15) = to<double>(dsqmix.find_block_def()->at(3)); // Q(GeV) corrects minor bug in SUSY-HIT SLHA reader.
    if (selmix.find_block_def()->size() >= 4) sd_leshouches2->qvalue(16) = to<double>(selmix.find_block_def()->at(3)); // Q(GeV) corrects minor bug in SUSY-HIT SLHA reader.
    for (int i=1; i<=6; ++i)
    {
      for (int j=1; j<=6; ++j)
      {
        std::vector<int> ij = initVector<int>(i,j);
        flavviolation->usqmix(i,j) = (usqmix[ij].is_data_line()) ? to<double>(usqmix.at(i,j)[2]) : 0.0;
        flavviolation->dsqmix(i,j) = (dsqmix[ij].is_data_line()) ? to<double>(dsqmix.at(i,j)[2]) : 0.0;
        sd_selectron->selmix(i,j)  = (selmix[ij].is_data_line()) ? to<double>(selmix.at(i,j)[2]) : 0.0;
      }
    }
    */

    // Set equivalent SLHA common blocks for HDecay.  Only differences are dimension of qvalues and zero vs unlikely for au, ad & ae.
    *slha_leshouches1_hdec = *sd_leshouches1;                       // SLHA1 block is identical in SDECAY and HDECAY.
    slha_leshouches2_hdec->imod = sd_leshouches2->imod;             // model; 1, 1 = SUGRA.
    slha_leshouches2_hdec->smval = sd_leshouches2->smval;           // SMINPUTS
    slha_leshouches2_hdec->extval = sd_leshouches2->extval;         // EXTPAR
    slha_leshouches2_hdec->massval = sd_leshouches2->massval;       // MASS
    slha_leshouches2_hdec->nmixval = sd_leshouches2->nmixval;       // NMIX
    slha_leshouches2_hdec->vmixval = sd_leshouches2->vmixval;       // VMIX
    slha_leshouches2_hdec->umixval = sd_leshouches2->umixval;       // UMIX
    slha_leshouches2_hdec->stopmixval = sd_leshouches2->stopmixval; // STOPMIX
    slha_leshouches2_hdec->sbotmixval = sd_leshouches2->sbotmixval; // SBOTMIX
    slha_leshouches2_hdec->staumixval = sd_leshouches2->staumixval; // STAUMIX
    slha_leshouches2_hdec->alphaval = sd_leshouches2->alphaval;     // ALPHA
    slha_leshouches2_hdec->hmixval = sd_leshouches2->hmixval;       // HMIX
    slha_leshouches2_hdec->gaugeval = sd_leshouches2->gaugeval;     // GAUGE
    slha_leshouches2_hdec->msoftval = sd_leshouches2->msoftval;     // MSOFT
    slha_leshouches2_hdec->yuval = sd_leshouches2->yuval;           // YU
    slha_leshouches2_hdec->ydval = sd_leshouches2->ydval;           // YD
    slha_leshouches2_hdec->yeval = sd_leshouches2->yeval;           // YE
    for (int i=1; i<=3; ++i) for (int j=1; j<=3; ++j)
    {
      if (sd_leshouches2->auval(i,j) == unlikely())                 // AU
      {
        slha_leshouches2_hdec->auval(i,j) = 0.0;
      }
      else
      {
       slha_leshouches2_hdec->auval(i,j) = sd_leshouches2->auval(i,j);
      }
      if (sd_leshouches2->adval(i,j) == unlikely())                 // AD
      {
        slha_leshouches2_hdec->adval(i,j) = 0.0;
      }
      else
      {
       slha_leshouches2_hdec->adval(i,j) = sd_leshouches2->adval(i,j);
      }
      if (sd_leshouches2->aeval(i,j) == unlikely())                 // AE
      {
        slha_leshouches2_hdec->aeval(i,j) = 0.0;
      }
      else
      {
       slha_leshouches2_hdec->aeval(i,j) = sd_leshouches2->aeval(i,j);
      }
    }
    for (int i=1; i<=20; ++i) slha_leshouches2_hdec->qvalue(i) = sd_leshouches2->qvalue(i); // Q(GeV)

    // Run SUSY-HIT
    sdecay();

  }

}
END_BE_NAMESPACE


// Initialisation function (definition)
BE_INI_FUNCTION
{
  SLHAstruct slha;

  // If the user provides a file list, just read in SLHA files for debugging and ignore the MSSM_spectrum dependency.
  if (runOptions->hasKey("debug_SLHA_filenames"))
  {
    static unsigned int counter = 0;
    std::vector<str> filenames = runOptions->getValue<std::vector<str> >("debug_SLHA_filenames");
    logger() << "Reading SLHA file: " << filenames[counter] << EOM;
    std::ifstream ifs(filenames[counter]);
    if(!ifs.good()) backend_error().raise(LOCAL_INFO, "SLHA file not found.");
    ifs >> slha;
    ifs.close();
    counter++;
    if (counter >= filenames.size()) counter = 0;

    // If the CKM entries are missing from the SLHA file, fill them in with defaults.
    SLHAea_add(slha, "VCKMIN", 1, 0.22537, "lambda", false);
    SLHAea_add(slha, "VCKMIN", 2, 0.814, "A", false);
    SLHAea_add(slha, "VCKMIN", 3, 0.117, "rhobar", false);
    SLHAea_add(slha, "VCKMIN", 4, 0.353, "etabar", false);
    // If the b pole mass is missing from the SLHA file, fill it in with a default.
    SLHAea_add(slha, "MASS",   5, 4.87878, "mb (pole)", false);
  }
  else // Use the actual spectrum object.
  {
    // Make sure the spectrum object is at the SUSY scale
    double scale = Dep::MSSM_spectrum->get_HE().GetScale();
    double susy_scale = Dep::MSSM_spectrum->get(Par::mass1, "susy_scale");
    // For some reason the precision on these is different, so they won't be exactly the same
    if (fabs(scale - susy_scale) > 1e-10) backend_error().raise(LOCAL_INFO, "MSSM_spectrum dependency is not at the SUSY scale.");

    // Get an SLHA1 object. SUSY-HIT is not SLHA2-compliant, despite its ability to deal with FV stop decays.
    slha = Dep::MSSM_spectrum->getSLHAea(1);

    // Check the tolerances for off-diagonal sfermion mixing.  It's a bit inefficient to redo this here,
    // given that these conversions have been done already in getSLHAea, but it's neater than
    // providing a generic means to pass the different tolerances through the spectrum object.
    const static double gtol = runOptions->getValueOrDef<double>(1e-2, "gauge_mixing_tolerance");
    const static bool gpterror = runOptions->getValueOrDef<bool>(true, "gauge_mixing_tolerance_invalidates_point_only");
    const static double ftol = runOptions->getValueOrDef<double>(1e-2, "family_mixing_tolerance");
    const static bool fpterror = runOptions->getValueOrDef<bool>(true, "family_mixing_tolerance_invalidates_point_only");
    const static str SLHA1_states[18] = {"~e_L", "~e_R", "~mu_L", "~mu_R", "~u_L", "~u_R", "~d_L", "~d_R", "~c_L",
                                         "~c_R", "~s_L", "~s_R", "~t_1", "~t_2", "~b_1", "~b_2", "~tau_1", "~tau_2"};
    for (int i = 0; i < 12; i++) slhahelp::mass_es_from_gauge_es(SLHA1_states[i], Dep::MSSM_spectrum->get_HE(), gtol, LOCAL_INFO, gpterror);
    for (int i = 12; i < 18; i++) slhahelp::mass_es_closest_to_family(SLHA1_states[i], Dep::MSSM_spectrum->get_HE(), ftol, LOCAL_INFO, fpterror);
  }

  // Get the W and Z widths.
  double W_width = 0.5*(Dep::W_plus_decay_rates->width_in_GeV + Dep::W_minus_decay_rates->width_in_GeV);
  double Z_width = Dep::Z_decay_rates->width_in_GeV;

  // Calculate decay rates
  run_susy_hit(slha, W_width, Z_width);

}
END_BE_INI_FUNCTION

Updated on 2024-07-18 at 13:53:36 +0000