file frontends/DarkSUSY_5_1_3.cpp
[No description available] More…
Functions
| Name | |
|---|---|
| if(scan_level ) | |
| if(runOptions-> hasKey"debug_SLHA_filenames") | |
| else | if(ModelInUse(“CMSSM”) and runOptions->getValueOrDef< bool >(false, “use_DS_isasugra”) ) Option use_DS_isasugra |
| else | if(ModelInUse(“MSSM63atQ”) |
| if((ModelInUse(“MSSM63atQ”) |
Attributes
| Name | |
|---|---|
| BE_NAMESPACE | |
| std::vector< double > | DSparticle_mass |
| std::vector< double > | GAMBITparticle_mass |
| END_BE_NAMESPACE | BE_INI_FUNCTION |
| bool | mssm_result |
Detailed Description
Author:
- Anders Kvellestad (anders.kvellestad@fys.uio.no)
- Pat Scott (p.scott@imperial.ac.uk)
- Christoph Weniger (c.weniger@uva.nl)
- Torsten Bringmann (torsten.bringmann@fys.uio.no)
- Lars A. Dal (l.a.dal@fys.uio.no)
- Joakim Edsjo (edsjo@fysik.su.se)
Date:
- 2013 Mar
- 2013 Apr 2015 Mar, Aug 2016 Feb
- 2013 Jul
- 2013 Jul, 2014 Mar, 2015 May, 2019 May
- 2014 Mar
- 2015 Aug, 2016 Mar
Frontend for DarkSUSY 5.1.3 backend
Authors (add name and date if you modify):
Functions Documentation
function if
if(
scan_level
)
function if
if(
runOptions-> hasKey"debug_SLHA_filenames"
)
Option debug_SLHA_filenames<std::vectorstd::string>: Optional override list of SLHA filenames used for backend initialization default
function if
else if(
ModelInUse("CMSSM") and runOptions->getValueOrDef< bool >(false, "use_DS_isasugra")
)
Option use_DS_isasugra
function if
else if(
ModelInUse("MSSM63atQ")||ModelInUse("CMSSM")
)
Option use_dsSLHAread
function if
if(
(ModelInUse("MSSM63atQ")||ModelInUse("CMSSM")) &&! mssm_result
)
Attributes Documentation
variable BE_NAMESPACE
BE_NAMESPACE {
const double min_DS_rwidth = 5.e-3;
variable DSparticle_mass
std::vector< double > DSparticle_mass;
variable GAMBITparticle_mass
std::vector< double > GAMBITparticle_mass;
variable BE_INI_FUNCTION
END_BE_NAMESPACE BE_INI_FUNCTION {
bool static scan_level = true;
variable mssm_result
bool mssm_result = false;
Source code
// GAMBIT: Global and Modular BSM Inference Tool
// *********************************************
/// \file
///
/// Frontend for DarkSUSY 5.1.3 backend
///
/// *********************************************
///
/// Authors (add name and date if you modify):
///
/// \author Anders Kvellestad
/// (anders.kvellestad@fys.uio.no)
/// \date 2013 Mar
///
/// \author Pat Scott
/// (p.scott@imperial.ac.uk)
/// \date 2013 Apr
/// 2015 Mar, Aug
/// 2016 Feb
///
/// \author Christoph Weniger
/// (c.weniger@uva.nl)
/// \date 2013 Jul
///
/// \author Torsten Bringmann
/// (torsten.bringmann@fys.uio.no)
/// \date 2013 Jul, 2014 Mar, 2015 May, 2019 May
///
/// \author Lars A. Dal
/// (l.a.dal@fys.uio.no)
/// \date 2014 Mar
///
/// \author Joakim Edsjo
/// (edsjo@fysik.su.se)
/// \date 2015 Aug, 2016 Mar
///
/// *********************************************
#include "gambit/Backends/frontend_macros.hpp"
#include "gambit/Backends/frontends/DarkSUSY_5_1_3.hpp"
#include "gambit/Utils/file_lock.hpp"
#include "gambit/Utils/mpiwrapper.hpp"
//#define DARKSUSY_DEBUG
// Some ad-hoc DarkSUSY global state.
BE_NAMESPACE
{
const double min_DS_rwidth = 5.e-3; // 0.5% to avoid numerical problems
std::vector<double> DSparticle_mass;
std::vector<double> GAMBITparticle_mass;
}
END_BE_NAMESPACE
// Initialisation function (definition)
BE_INI_FUNCTION
{
// Initialize DarkSUSY (only) if run for the first time
bool static scan_level = true;
if (scan_level)
{
// Do the call to dsinit one-by-one for each MPI process, as DarkSUSY loads up
// HiggsBounds, which writes files at init then reads them back in later.
Utils::ProcessLock mylock("DarkSUSY_" STRINGIFY(SAFE_VERSION) "_init");
mylock.get_lock();
dsinit();
mylock.release_lock();
dsrdinit();
// Initialize yield tables for use in cascade decays (initialize more if needed)
dshainit(151); // Initialize positron tables
dshainit(152); // Initialize gamma ray tables
dshainit(154); // Initialize antiproton tables
// Note, for DarkSUSY 5, anti-deuteron does not need special initialization
// as they are calculated from antiproton fluxes
// Call dshayield for first call initialization of variables
double tmp1 = 100.0;
double tmp2 = 10.0;
int tmp3 = 15;
int tmp4 = 152;
int tmp5 = 0;
dshayield(tmp1,tmp2,tmp3,tmp4,tmp5);
scan_level = false;
}
// Initialization function for a given MSSM point
// (previous capability DarkSUSY_PointInit)
bool mssm_result = false;
// 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;
logger() << LogTags::debug <<
"Initializing DarkSUSY via debug_SLHA_filenames option." << EOM;
/// Option debug_SLHA_filenames<std::vector<std::string>>: Optional override list of SLHA filenames used for backend initialization default
std::vector<str> filenames = runOptions->getValue<std::vector<str> >("debug_SLHA_filenames");
const char * filename = filenames[counter].c_str();
int len = filenames[counter].length();
int flag = 15;
dsSLHAread(byVal(filename),flag,byVal(len));
dsprep();
counter++;
if (counter >= filenames.size()) counter = 0;
mssm_result = true;
}
// CMSSM with DS-internal ISASUGRA (should be avoided, only for debugging)
/// Option use_DS_isasugra<bool>: Use DS internal isasugra for parameter running (false)
else if (ModelInUse("CMSSM") and runOptions->getValueOrDef<bool>(false, "use_DS_isasugra"))
{
// Setup mSUGRA model from CMSSM parameters
double am0 = *Param["M0"]; // m0
double amhf = *Param["M12"]; // m_1/2
double aa0 = *Param["A0"]; // A0
double asgnmu = *Param["SignMu"]; // sign(mu)
double atanbe = *Param["TanBeta"]; // tan(beta)
logger() << "Initializing DarkSUSY via dsgive_model_isasugra:" << EOM;
logger() << " m0 =" << am0 << std::endl;
logger() << " m_1/2 =" << amhf << std::endl;
logger() << " A0 =" << aa0 << std::endl;
logger() << " sign(mu) =" << asgnmu << std::endl;
logger() << " tan(beta) =" << atanbe << EOM;
dsgive_model_isasugra(am0, amhf, aa0, asgnmu, atanbe);
int unphys, hwarning;
dssusy_isasugra(unphys, hwarning);
if (unphys < 0)
{
backend_warning().raise(LOCAL_INFO, "Model point is theoretically inconsistent (DarkSUSY).");
invalid_point().raise("Model point is theoretically inconsistent (DarkSUSY).");
mssm_result = false;
}
else if (unphys > 0)
{
backend_warning().raise(LOCAL_INFO, "Neutralino is not the LSP (DarkSUSY).");
invalid_point().raise("Neutralino is not the LSP (DarkSUSY).");
mssm_result = false;
}
else if (hwarning != 0)
{
backend_warning().raise(LOCAL_INFO,
"Radiative corrections in Higgs sector "
"outside range of validity (DarkSUSY).");
mssm_result = true;
}
else
{
mssm_result = true;
}
}
else if (ModelInUse("MSSM63atQ") || ModelInUse("CMSSM"))
{
SLHAstruct mySLHA;
/// Option use_dsSLHAread<bool>: Use DS internal SLHA reader to initialize backend (false)
bool use_dsSLHAread = runOptions->getValueOrDef<bool>(false, "use_dsSLHAread");
int slha_version = 2;
const Spectrum& mySpec = *Dep::MSSM_spectrum;
try
{
mySLHA = mySpec.getSLHAea(2);
}
catch(Gambit::exception& e)
{
slha_version = 1;
mySLHA = mySpec.getSLHAea(1);
use_dsSLHAread = true;
}
// Use an actual SLHA file. DarkSUSY is on its own wrt (s)particle widths this way.
if (use_dsSLHAread || slha_version == 1)
{
if (!use_dsSLHAread) {backend_error().raise(LOCAL_INFO,
"An SLHA1 spectrum requires use of the DarkSUSY SLHA reader rather than the diskless\n"
"GAMBIT DarkSUSY initialization. To enable the DarkSUSY SLHA reader, set the option\n"
"use_dsSLHAread for the function DarkSUSY_PointInit_MSSM to true.");}
int rank = 0;
#ifdef WITH_MPI
if(GMPI::Is_initialized())
{
GMPI::Comm comm;
rank = comm.Get_rank();
}
#endif
// Add model select block to inform DS about 6x6 mixing
if (slha_version == 2)
{
SLHAea::Block modsel_block("MODSEL");
modsel_block.push_back("BLOCK MODSEL");
modsel_block.push_back("6 3 # FV");
mySLHA.push_back(modsel_block);
}
// Set filename
std::string fstr = "DarkBit_temp_";
fstr += std::to_string(rank) + ".slha";
// Dump SLHA onto disk
std::ofstream ofs(fstr);
ofs << mySLHA;
ofs.close();
// Initialize SUSY spectrum from SLHA
int len = fstr.size();
int flag = 15;
const char * filename = fstr.c_str();
logger() << LogTags::debug << "Initializing DarkSUSY via SLHA." << EOM;
dsSLHAread(byVal(filename),flag,byVal(len));
//The following used to be a separate capability dsprep up to DS5
//(as also specified in the manual)
dsprep();
mssm_result = true;
}
// Do pure diskless SLHA initialisation, including (s)particle widths from GAMBIT.
else
{
if (init_diskless(mySLHA, *Dep::decay_rates) == 0 )
{
logger() << LogTags::debug << "Using diskless SLHA interface to DarkSUSY." << EOM;
dsprep();
mssm_result = true;
}
}
}
if ( (ModelInUse("MSSM63atQ") || ModelInUse("CMSSM")) && !mssm_result )
{
backend_warning().raise(LOCAL_INFO,
"DarkSUSY point initialization failed.");
invalid_point().raise("DarkSUSY point initialization failed.");
}
}
END_BE_INI_FUNCTION
// Convenience functions (definitions)
BE_NAMESPACE
{
/// Shortcut for dealing with SLHA blocks
void required_block(const std::string& name, const SLHAea::Coll& slha)
{
if (slha.find(name) != slha.end()) return;
else backend_error().raise(LOCAL_INFO, "Sorry, DarkSUSY needs SLHA block: " + name + ".\n"
"If you tried to read in a debug SLHA file with missing entries, \n"
"then sort out your SLHA file so that it is readable by DarkSUSY! ");
}
/// Sets DarkSUSY's internal common blocks with some of the properties required to compute neutrino
/// yields for a generic WIMP. Remaining internal variables are internal to this frontend.
void dsgenericwimp_nusetup(const double (&annihilation_bf)[29], const double (&Higgs_decay_BFs_neutral)[29][3],
const double (&Higgs_decay_BFs_charged)[15], const double (&Higgs_masses_neutral)[3], const double &Higgs_mass_charged,
const double &mwimp)
{
// Transfer WIMP mass common block.
wabranch->wamwimp = mwimp;
// Transfer branching fractions to WIMP annihilation common blocks.
// For channel indices, see dswayieldone.f
for (int i=1; i<=29; i++)
{
wabranch->wabr(i) = annihilation_bf[i-1];
}
// Transfer Higgs decay branching fractions (not widths) to Higgs decay common blocks.
// The total width is not relevant, as all Higgs decay in flight eventually, so
// only the BFs are needed to calculate the yields into neutrinos from decays in flight.
for (int i=1; i<=3; i++) // Loop over the neutral Higgses
{
for (int j=1; j<=29; j++) // Loop over the known decay channels
{
wabranch->was0br(j,i) = Higgs_decay_BFs_neutral[j-1][i-1];
}
}
for (int i=1; i<=15; i++) // Loop over the known charged Higgs decay channels
{
wabranch->wascbr(i) = Higgs_decay_BFs_charged[i-1];
}
// Transfer Higgs masses to common blocks.
for (int i=1; i<=3; i++) // Loop over the neutral Higgses
{
wabranch->was0m(i) = Higgs_masses_neutral[i-1]; // Neutral Higgses
}
wabranch->wascm = Higgs_mass_charged; // Charged Higgses
// Tell DarkSUSY we've taken care of business.
wabranch->dswasetupcalled = true;
}
/// Returns neutrino yields at the top of the atmosphere,
/// in m^-2 GeV^-1 annihilation^-1. Provided here for
/// interfacing with nulike.
/// --> log10Enu log_10(neutrino energy/GeV)
/// --> p p=1 for neutrino yield, p=2 for nubar yield
/// --> context void pointer (ignored)
double neutrino_yield(const double& log10E, const int& ptype, void*&)
{
int istat = 0;
const char object[3] = "su";
double result = 1e-30 * dsntmuonyield(pow(10.0,log10E),10.0,object[0],3,1,ptype,istat);
if ((istat bitand 1) == 1)
{
if (not piped_warnings.inquire()) // Don't bother re-raising a warning if it's already been done since the last .check().
piped_warnings.request(LOCAL_INFO, "Neutrino yield from Sun is lower bound; likelihood will be conservative.");
}
if ((istat bitand 4) == 4)
{
if (not piped_warnings.inquire()) // Don't bother re-raising a warning if it's already been done since the last .check().
piped_warnings.request(LOCAL_INFO, "DarkSUSY's dswayield_int didn't converge. This occasionally happens "
"due to finite statistics in the nu yield tables from Pythia. "
"This is benign (the missing integrals are always negligible).");
}
if (istat > 4)
{
std::ostringstream err;
err << "Error from DarkSUSY::dswayield functions in neutrino flux calculation. istat = " << istat;
piped_errors.request(LOCAL_INFO, err.str());
}
return result;
}
/// Translates GAMBIT string identifiers to the SUSY
/// particle codes used internally in DS (as stored in common block /pacodes/)
// FIXME: add channel codes!
int DSparticle_code(const str& particleID)
{
int kpart;
if (particleID=="nu_e" or particleID=="nubar_e"){
kpart=1;
}else if (particleID=="e-_1" or particleID=="e+_1"){
kpart=2;
}else if (particleID=="nu_mu" or particleID=="nubar_mu"){
kpart=3;
}else if (particleID=="e-_2" or particleID=="e+_2"){
kpart=4;
}else if (particleID=="nu_tau" or particleID=="nubar_tau"){
kpart=5;
}else if (particleID=="e-_3" or particleID=="e+_3"){
kpart=6;
}else if (particleID=="u_1" or particleID=="ubar_1"){
kpart=7;
}else if (particleID=="d_1" or particleID=="dbar_1"){
kpart=8;
}else if (particleID=="u_2" or particleID=="ubar_2"){
kpart=9;
}else if (particleID=="d_2" or particleID=="dbar_2"){
kpart=10;
}else if (particleID=="u_3" or particleID=="ubar_3"){
kpart=11;
}else if (particleID=="d_3" or particleID=="dbar_3"){
kpart=12;
}else if (particleID=="gamma"){
kpart=13;
}else if (particleID=="W+" or particleID=="W-"){
kpart=14;
}else if (particleID=="Z0"){
kpart=15;
}else if (particleID=="g"){
kpart=16;
}else if (particleID=="h0_1"){
kpart=18;
}else if (particleID=="h0_2"){
kpart=17;
}else if (particleID=="A0"){
kpart=19;
}else if (particleID=="H+" or particleID=="H-"){
kpart=20;
}else if (particleID=="~nu_1" or particleID=="~nubar_1"){
kpart=21;
}else if (particleID=="~e-_1" or particleID=="~e+_1"){
kpart=22;
}else if (particleID=="~e-_4" or particleID=="~e+_4"){
kpart=23;
}else if (particleID=="~nu_2" or particleID=="~nubar_2"){
kpart=24;
}else if (particleID=="~e-_2" or particleID=="~e+_2"){
kpart=25;
}else if (particleID=="~e-_5" or particleID=="~e+_5"){
kpart=26;
}else if (particleID=="~nu_3" or particleID=="~nubar_3"){
kpart=27;
}else if (particleID=="~e-_3" or particleID=="~e+_3"){
kpart=28;
}else if (particleID=="~e-_6" or particleID=="~e+_6"){
kpart=29;
}else if (particleID=="~u_1" or particleID=="~ubar_1"){
kpart=30;
}else if (particleID=="~u_4" or particleID=="~ubar_4"){
kpart=31;
}else if (particleID=="~d_1" or particleID=="~dbar_1"){
kpart=32;
}else if (particleID=="~d_4" or particleID=="~dbar_4"){
kpart=33;
}else if (particleID=="~u_2" or particleID=="~ubar_2"){
kpart=34;
}else if (particleID=="~u_5" or particleID=="~ubar_5"){
kpart=35;
}else if (particleID=="~d_2" or particleID=="~dbar_2"){
kpart=36;
}else if (particleID=="~d_5" or particleID=="~dbar_5"){
kpart=37;
}else if (particleID=="~u_3" or particleID=="~ubar_3"){
kpart=38;
}else if (particleID=="~u_6" or particleID=="~ubar_6"){
kpart=39;
}else if (particleID=="~d_3" or particleID=="~dbar_3"){
kpart=40;
}else if (particleID=="~d_6" or particleID=="~dbar_6"){
kpart=41;
}else if (particleID=="~chi0_1"){
kpart=42;
}else if (particleID=="~chi0_2"){
kpart=43;
}else if (particleID=="~chi0_3"){
kpart=44;
}else if (particleID=="~chi0_4"){
kpart=45;
}else if (particleID=="~chi+_1" or particleID=="~chi-_1"){
kpart=46;
}else if (particleID=="~chi+_2" or particleID=="~chi-_2"){
kpart=47;
}else if (particleID=="~g"){
kpart=48;
} else{
std::ostringstream err;
err << "ERROR: translation into DS particle code not implemented "
<< "for string identifier " << particleID;
backend_error().raise(LOCAL_INFO, err.str());
kpart=-100;
}
return kpart;
}
/// Initialise an MSSM model in DarkSUSY from an SLHAea object and a DecayTable.
/// Closely mimics the DarkSUSY routine in dsfromslha.F, except that it
/// hands over a better b pole mass, explicit decay info, and a better
/// approximation of the CKM matrix from Wolfenstein parameters. Throughout
/// this routine, there are pieces of code commented out that would need to be
/// re-added to emulate dsfromslha.F exactly, if and only if the decays are not
/// set by GAMBIT.
int init_diskless(const SLHAstruct &mySLHA, const DecayTable &myDecays)
{
using SLHAea::to;
DS5_PACODES *DSpart = &(*pacodes);
// Define required blocks and raise an error if a block is missing
required_block("SMINPUTS", mySLHA);
required_block("VCKMIN", mySLHA);
required_block("MSOFT", mySLHA);
required_block("MASS", mySLHA);
required_block("NMIX", mySLHA);
required_block("VMIX", mySLHA);
required_block("UMIX", mySLHA);
required_block("ALPHA", mySLHA);
required_block("HMIX", mySLHA);
required_block("YU", mySLHA);
required_block("YD", mySLHA);
required_block("YE", mySLHA);
required_block("MSL2", mySLHA);
required_block("MSE2", mySLHA);
required_block("MSQ2", mySLHA);
required_block("MSD2", mySLHA);
required_block("MSU2", mySLHA);
required_block("TD", mySLHA);
required_block("TU", mySLHA);
required_block("TE", mySLHA);
required_block("USQMIX", mySLHA);
required_block("DSQMIX", mySLHA);
required_block("SELMIX", mySLHA);
required_block("SNUMIX", mySLHA);
// Make sure the b pole mass is present in the MASS block
if (mySLHA.at("MASS").find(initVector<str>("5")) == mySLHA.at("MASS").end())
backend_error().raise(LOCAL_INFO, "DarkSUSY init_diskless needs b pole mass entry (5) in SLHA(ea) MASS block.");
// Do some initial DarkSUSY housekeeping.
dsmssmzero(); // zero all the MSSM parameters and variables
mssmtype->modeltype = 0; // tell DarkSUSY that we are working in the general MSSM
// To match the DarkSUSY SLHAreader, you would need to also set
// mssmswitch->higwid = 1; // tell DarkSUSY not to use FeynHiggs for Higgs widths.
// Block SMINPUTS
couplingconstants->alphem = 1./to<double>(mySLHA.at("SMINPUTS").at(1).at(1)); // 1/alpha_{QED}
smruseful->alph3mz = to<double>(mySLHA.at("SMINPUTS").at(3).at(1)); // alpha_s @ MZ
smruseful->gfermi = to<double>(mySLHA.at("SMINPUTS").at(2).at(1)); // Fermi constant
mspctm->mass(DSparticle_code("Z0")) = to<double>(mySLHA.at("SMINPUTS").at(4).at(1)); // Z boson mass
// Here we set the masses to be used in DarkSUSY. Note that all masses in the mspctm->mass block
// must match those in the ProcessCatalog in DarkBit, as these are used to define the kinematic
// edges used in relic density integrations and similar within DarkSUSY. Mostly these should be
// pole masses, except in cases where that is not possible (i.e. light quarks).
// Lepton masses
mspctm->mass(DSpart->kl(1)) = to<double>(mySLHA.at("SMINPUTS").at(11).at(1)); // electron pole mass
mspctm->mass(DSpart->kl(2)) = to<double>(mySLHA.at("SMINPUTS").at(13).at(1)); // muon pole mass
mspctm->mass(DSpart->kl(3)) = to<double>(mySLHA.at("SMINPUTS").at(7).at(1)); // tau pole mass
mspctm->mass(DSpart->knu(1)) = to<double>(mySLHA.at("SMINPUTS").at(12).at(1)); // nu_1 pole mass
mspctm->mass(DSpart->knu(2)) = to<double>(mySLHA.at("SMINPUTS").at(14).at(1)); // nu_2 pole mass
mspctm->mass(DSpart->knu(3)) = to<double>(mySLHA.at("SMINPUTS").at(8).at(1)); // nu_3 pole mass
// Quark masses as defined in SLHA2
mspctm->mu2gev = to<double>(mySLHA.at("SMINPUTS").at(22).at(1)); // up quark mass @ 2 GeV
mspctm->md2gev = to<double>(mySLHA.at("SMINPUTS").at(21).at(1)); // down quark mass @ 2 GeV
mspctm->ms2gev = to<double>(mySLHA.at("SMINPUTS").at(23).at(1)); // strange mass @ 2 GeV
mspctm->mcmc = to<double>(mySLHA.at("SMINPUTS").at(24).at(1)); // charm mass at m_c
mspctm->mbmb = to<double>(mySLHA.at("SMINPUTS").at(5).at(1)); // bottom mass at m_b
mspctm->mass(DSpart->kqu(3)) = to<double>(mySLHA.at("SMINPUTS").at(6).at(1)); // top pole mass
// Do the DarkSUSY-style sin^2 theta_W calculation (will be overwritten later).
smruseful->s2thw=dsgf2s2thw(smruseful->gfermi, couplingconstants->alphem, mspctm->mass(DSparticle_code("Z0")), mspctm->mass(DSpart->kqu(3)),1);
// Set other internal quark masses for DarkSUSY
dsfindmtmt(); // top mass at mt
mspctm->mass(DSpart->kqu(1)) = mspctm->mu2gev; // use 2GeV u mass as proxy for pole
mspctm->mass(DSpart->kqd(1)) = mspctm->md2gev; // use 2GeV d mass as proxy for pole
mspctm->mass(DSpart->kqd(2)) = mspctm->ms2gev; // use 2GeV s mass as proxy for pole
mspctm->mass(DSpart->kqu(2)) = dsmqpole4loop(DSpart->kqu(2),mspctm->mcmc); // use DarkSUSY internal routine to get mc pole
mspctm->mass(DSpart->kqd(3)) = to<double>(mySLHA.at("MASS").at(5).at(1)); // the GAMBIT way to get the bottom pole mass
//mspctm->mass(DSpart->kqd(3)) = dsmqpole4loop(DSpart->kqd(3),mspctm->mbmb); // the DarkSUSY SLHAreader way to get mb pole
// Block MINPAR we skip, it is not needed
// Block MSOFT
mssmpar->m1 = to<double>(mySLHA.at("MSOFT").at(1).at(1));
mssmpar->m2 = to<double>(mySLHA.at("MSOFT").at(2).at(1));
mssmpar->m3 = to<double>(mySLHA.at("MSOFT").at(3).at(1));
mssmpar->mass2l(1) = pow(to<double>(mySLHA.at("MSOFT").at(31).at(1)),2);
mssmpar->mass2l(2) = pow(to<double>(mySLHA.at("MSOFT").at(32).at(1)),2);
mssmpar->mass2l(3) = pow(to<double>(mySLHA.at("MSOFT").at(33).at(1)),2);
mssmpar->mass2e(1) = pow(to<double>(mySLHA.at("MSOFT").at(34).at(1)),2);
mssmpar->mass2e(2) = pow(to<double>(mySLHA.at("MSOFT").at(35).at(1)),2);
mssmpar->mass2e(3) = pow(to<double>(mySLHA.at("MSOFT").at(36).at(1)),2);
mssmpar->mass2q(1) = pow(to<double>(mySLHA.at("MSOFT").at(41).at(1)),2);
mssmpar->mass2q(2) = pow(to<double>(mySLHA.at("MSOFT").at(42).at(1)),2);
mssmpar->mass2q(3) = pow(to<double>(mySLHA.at("MSOFT").at(43).at(1)),2);
mssmpar->mass2u(1) = pow(to<double>(mySLHA.at("MSOFT").at(44).at(1)),2);
mssmpar->mass2u(2) = pow(to<double>(mySLHA.at("MSOFT").at(45).at(1)),2);
mssmpar->mass2u(3) = pow(to<double>(mySLHA.at("MSOFT").at(46).at(1)),2);
mssmpar->mass2d(1) = pow(to<double>(mySLHA.at("MSOFT").at(47).at(1)),2);
mssmpar->mass2d(2) = pow(to<double>(mySLHA.at("MSOFT").at(48).at(1)),2);
mssmpar->mass2d(3) = pow(to<double>(mySLHA.at("MSOFT").at(49).at(1)),2);
// Block HMIX
mssmpar->mu = to<double>(mySLHA.at("HMIX").at(1).at(1));
mssmpar->tanbe = to<double>(mySLHA.at("HMIX").at(2).at(1));
// A boson mass
mspctm->mass(DSparticle_code("A0")) = to<double>(mySLHA.at("MASS").at(36).at(1));
mssmpar->ma = mspctm->mass(DSparticle_code("A0"));
// Now set up some defaults (part of it will be overwritten later)
dssuconst();
mssmiuseful->lsp = DSpart->kn(1);
mssmiuseful->kln = 1;
// To match the DS SLHAreader, you would then need to call
// int unphys, hwarning;
// dsspectrum(unphys, hwarning);
// CKM matrix read from VCKMIN - Wolfenstein parameters.
double lambda = to<double>(mySLHA.at("VCKMIN").at(1).at(1)); // Wolfenstein lambda
double A = to<double>(mySLHA.at("VCKMIN").at(2).at(1)); // Wolfenstein A
double rhobar = to<double>(mySLHA.at("VCKMIN").at(3).at(1)); // Wolfenstein rhobar
double etabar = to<double>(mySLHA.at("VCKMIN").at(4).at(1)); // Wolfenstein etabar
// Use Wolfenstein converter to get the VCKM matrix.
mixing->ckm(1,1) = Spectrum::Wolf2V_ud(lambda,A,rhobar,etabar);
mixing->ckm(1,2) = Spectrum::Wolf2V_us(lambda,A,rhobar,etabar);
mixing->ckm(1,3) = Spectrum::Wolf2V_ub(lambda,A,rhobar,etabar);
mixing->ckm(2,1) = Spectrum::Wolf2V_cd(lambda,A,rhobar,etabar);
mixing->ckm(2,2) = Spectrum::Wolf2V_cs(lambda,A,rhobar,etabar);
mixing->ckm(2,3) = Spectrum::Wolf2V_cb(lambda,A,rhobar,etabar);
mixing->ckm(3,1) = Spectrum::Wolf2V_td(lambda,A,rhobar,etabar);
mixing->ckm(3,2) = Spectrum::Wolf2V_ts(lambda,A,rhobar,etabar);
mixing->ckm(3,3) = Spectrum::Wolf2V_tb(lambda,A,rhobar,etabar);
// The lower-order DarkSUSY SLHAreader way of doing it; not as good, but this is what you'd need to use to match DS exactly.
//sckm->ckms12 = lambda;
//sckm->ckms23 = A*pow(sckm->ckms12,2);
//std::complex<double> aux(rhobar, etabar);
//aux = aux * (sckm->ckms23/pow(sckm->ckms12,2)) * pow(sckm->ckms12,3);
//sckm->ckms13 = std::norm(aux);
//sckm->ckmdelta = std::arg(aux);
//dssuconst_ckm();
// In principle, we might want to change to VCKM instead of VCKMIN, if VCKM is present. Like this:
// sckm->ckms12 = to<double>(mySLHA.at("VCKM").at(1).at(1));
// sckm->ckms23 = to<double>(mySLHA.at("VCKM").at(2).at(1))*sckm.ckms12**2;
// sckm->ckmdelta = 0;
// for (int i=1; i<=3; i++) for (int j=1; j<=3; j++)
// {
// mixing->ckm(i,j) = to<double_complex>(mySLHA.at("VCKM").at(i,j).at(2));
// }
// OK, we now have to enforce the tree-level condition for unitarity.
// We then have a choice of calculating both sin^2 theta_W and MW
// from alpha, MZ and GF as we normally do in DarkSUSY. This line would
// enforce that:
// mspctm->mass(DSparticle_code("W+"))=mspctm->mass(DSparticle_code("Z0"))*sqrt(1.0-smruseful->s2thw);
// However, it is more prudent to take the value of MW from the SLHA file
// as given (read in earlier), and instead enforce the tree-level condition
// by redefining sin^2 theta_W. That we do here:
mspctm->mass(DSparticle_code("W+")) = to<double>(mySLHA.at("MASS").at(24).at(1)); // W boson mass
smruseful->s2thw=1.0-pow(mspctm->mass(DSparticle_code("W+"))/mspctm->mass(DSparticle_code("Z0")),2);
// Higgs masses. Note h0_1 is the lightest CP-even neutral higgs, and h2_0 the heavier.
mspctm->mass(DSparticle_code("h0_1")) = to<double>(mySLHA.at("MASS").at(25).at(1));
mspctm->mass(DSparticle_code("h0_2")) = to<double>(mySLHA.at("MASS").at(35).at(1));
mspctm->mass(DSparticle_code("H+")) = to<double>(mySLHA.at("MASS").at(37).at(1));
mspctm->mass(DSparticle_code("A0")) = to<double>(mySLHA.at("MASS").at(36).at(1));
// SUSY particles
mspctm->mass(DSpart->ksnu(1)) = to<double>(mySLHA.at("MASS").at(1000012).at(1));
mspctm->mass(DSpart->ksnu(2)) = to<double>(mySLHA.at("MASS").at(1000014).at(1));
mspctm->mass(DSpart->ksnu(3)) = to<double>(mySLHA.at("MASS").at(1000016).at(1));
mspctm->mass(DSpart->ksl(1)) = to<double>(mySLHA.at("MASS").at(1000011).at(1));
mspctm->mass(DSpart->ksl(2)) = to<double>(mySLHA.at("MASS").at(1000013).at(1));
mspctm->mass(DSpart->ksl(3)) = to<double>(mySLHA.at("MASS").at(1000015).at(1));
mspctm->mass(DSpart->ksl(4)) = to<double>(mySLHA.at("MASS").at(2000011).at(1));
mspctm->mass(DSpart->ksl(5)) = to<double>(mySLHA.at("MASS").at(2000013).at(1));
mspctm->mass(DSpart->ksl(6)) = to<double>(mySLHA.at("MASS").at(2000015).at(1));
mspctm->mass(DSpart->ksqu(1)) = to<double>(mySLHA.at("MASS").at(1000002).at(1));
mspctm->mass(DSpart->ksqu(2)) = to<double>(mySLHA.at("MASS").at(1000004).at(1));
mspctm->mass(DSpart->ksqu(3)) = to<double>(mySLHA.at("MASS").at(1000006).at(1));
mspctm->mass(DSpart->ksqu(4)) = to<double>(mySLHA.at("MASS").at(2000002).at(1));
mspctm->mass(DSpart->ksqu(5)) = to<double>(mySLHA.at("MASS").at(2000004).at(1));
mspctm->mass(DSpart->ksqu(6)) = to<double>(mySLHA.at("MASS").at(2000006).at(1));
mspctm->mass(DSpart->ksqd(1)) = to<double>(mySLHA.at("MASS").at(1000001).at(1));
mspctm->mass(DSpart->ksqd(2)) = to<double>(mySLHA.at("MASS").at(1000003).at(1));
mspctm->mass(DSpart->ksqd(3)) = to<double>(mySLHA.at("MASS").at(1000005).at(1));
mspctm->mass(DSpart->ksqd(4)) = to<double>(mySLHA.at("MASS").at(2000001).at(1));
mspctm->mass(DSpart->ksqd(5)) = to<double>(mySLHA.at("MASS").at(2000003).at(1));
mspctm->mass(DSpart->ksqd(6)) = to<double>(mySLHA.at("MASS").at(2000005).at(1));
// Neutralinos carry the sign of the eigenvalue, as we need it for NMIX later
mspctm->mass(DSpart->kn(1)) = to<double>(mySLHA.at("MASS").at(1000022).at(1));
mspctm->mass(DSpart->kn(2)) = to<double>(mySLHA.at("MASS").at(1000023).at(1));
mspctm->mass(DSpart->kn(3)) = to<double>(mySLHA.at("MASS").at(1000025).at(1));
mspctm->mass(DSpart->kn(4)) = to<double>(mySLHA.at("MASS").at(1000035).at(1));
// Charginos
mspctm->mass(DSpart->kcha(1)) = to<double>(mySLHA.at("MASS").at(1000024).at(1));
mspctm->mass(DSpart->kcha(2)) = to<double>(mySLHA.at("MASS").at(1000037).at(1));
// Gluino
mspctm->mass(DSparticle_code("~g")) = to<double>(mySLHA.at("MASS").at(1000021).at(1));
// Gravitino (not implemented in DarkSUSY)
// mspctm->mass(DSpart->k...) = to<double>(mySLHA.at("MASS").at(1000039).at(1));
// Block NMIX
for (int i=1; i<=4; i++)
{
for (int j=1; j<=4; j++)
{
mssmmixing->neunmx(i,j)=to<double>(mySLHA.at("NMIX").at(i,j).at(2));
}
}
// Make NMIX imaginary if mass eigenvalue is negative
for (int i=1; i<=4; i++)
{
for (int j=1; j<=4; j++)
{
if (mspctm->mass(DSpart->kn(i)) < 0)
{
mssmmixing->neunmx(i,j).im = mssmmixing->neunmx(i,j).re;
mssmmixing->neunmx(i,j).re = 0.0;
}
}
mspctm->mass(DSpart->kn(i)) = std::abs(mspctm->mass(DSpart->kn(i)));
}
// Block UMIX
for (int i=1; i<=2; i++)
{
for (int j=1; j<=2; j++)
{
mssmmixing->chaumx(i,j)=to<double>(mySLHA.at("UMIX").at(i,j).at(2));
}
}
// Block VMIX
for (int i=1; i<=2; i++)
{
for (int j=1; j<=2; j++)
{
mssmmixing->chavmx(i,j)=to<double>(mySLHA.at("VMIX").at(i,j).at(2));
}
}
// Block ALPHA
mssmmixing->alpha = to<double>(mySLHA.at("ALPHA").back().at(0)); // Higgs mixing angle
//If you want to exactly match the DarkSUSY SLHAreader, you should also now run
//dssuconst_yukawa_running();
// Block YE, YU, YD - Yukawas
for (int i=1; i<=3; i++)
{
couplingconstants->yukawa(DSpart->kl(i))=to<double>(mySLHA.at("YE").at(i,i).at(2));
couplingconstants->yukawa(DSpart->kqu(i))=to<double>(mySLHA.at("YU").at(i,i).at(2));
couplingconstants->yukawa(DSpart->kqd(i))=to<double>(mySLHA.at("YD").at(i,i).at(2));
}
// Block MSL2, MSE2, MSQ2, MSU2, MSD2
for (int i=1; i<=3; i++)
{
mssmpar->mass2l(i) = to<double>(mySLHA.at("MSL2").at(i,i).at(2));
mssmpar->mass2e(i) = to<double>(mySLHA.at("MSE2").at(i,i).at(2));
mssmpar->mass2q(i) = to<double>(mySLHA.at("MSQ2").at(i,i).at(2));
mssmpar->mass2u(i) = to<double>(mySLHA.at("MSU2").at(i,i).at(2));
mssmpar->mass2d(i) = to<double>(mySLHA.at("MSD2").at(i,i).at(2));
}
// BLOCK TE, TU and TD. I read these instead of AE, AU, AD.
for (int i=1; i<=3; i++)
{
mssmpar->asofte(i)=to<double>(mySLHA.at("TE").at(i,i).at(2))/couplingconstants->yukawa(DSpart->kl(i));
mssmpar->asoftu(i)=to<double>(mySLHA.at("TU").at(i,i).at(2))/couplingconstants->yukawa(DSpart->kqu(i));
mssmpar->asoftd(i)=to<double>(mySLHA.at("TD").at(i,i).at(2))/couplingconstants->yukawa(DSpart->kqd(i));
}
// Block SNUMIX
for (int i=1; i<=3; i++)
{
for (int j=1; j<=3; j++)
{
mssmmixing->slulmx(i,j) = to<double>(mySLHA.at("SNUMIX").at(i,j).at(2));
}
}
// Block SELMIX
for (int i=1; i<=6; i++)
{
for (int j=1; j<=3; j++)
{
mssmmixing->sldlmx(i,j) = to<double>(mySLHA.at("SELMIX").at(i,j).at(2));
mssmmixing->sldrmx(i,j) = to<double>(mySLHA.at("SELMIX").at(i,j+3).at(2));
}
}
// Block USQMIX
for (int i=1; i<=6; i++)
{
for (int j=1; j<=3; j++)
{
mssmmixing->squlmx(i,j) = to<double>(mySLHA.at("USQMIX").at(i,j).at(2));
mssmmixing->squrmx(i,j) = to<double>(mySLHA.at("USQMIX").at(i,j+3).at(2));
}
}
// Block DSQMIX
for (int i=1; i<=6; i++)
{
for (int j=1; j<=3; j++)
{
mssmmixing->sqdlmx(i,j) = to<double>(mySLHA.at("DSQMIX").at(i,j).at(2));
mssmmixing->sqdrmx(i,j) = to<double>(mySLHA.at("DSQMIX").at(i,j+3).at(2));
}
}
// Do flavour reordering for SLHA2 compatibility
dsorder_flavour();
// Set up SUSY vertices
dsvertx();
// At this point, if you wanted to match the DarkSUSY SLHAreader, you would also call
// dshigwid();
// dsspwid();
// which just fudge a few widths...but we won't do that, because we can get real decay widths from DecayBit.
// Set up Higgs widths. h1_0 is the lightest CP even Higgs in GAMBIT (opposite to DS).
widths->width(DSparticle_code("h0_1")) = myDecays.at(std::pair<int,int>(25,0)).width_in_GeV;
widths->width(DSparticle_code("h0_2")) = myDecays.at(std::pair<int,int>(35,0)).width_in_GeV;
widths->width(DSparticle_code("A0")) = myDecays.at(std::pair<int,int>(36,0)).width_in_GeV;
widths->width(DSparticle_code("H+")) = myDecays.at(std::pair<int,int>(37,0)).width_in_GeV;
// Set up Higgs partial widths.
const static std::vector< std::vector<str> > charged_channels = DS_charged_h_decay_channels();
const static std::vector< std::vector<str> > neutral_channels = DS_neutral_h_decay_channels();
const static std::vector<str> sister_chan = initVector<str>("W+", "H-");
const static std::vector<str> missing_chan = initVector<str>("W-", "H+");
const DecayTable::Entry& h01 = myDecays.at(std::pair<int,int>(25,0));
const DecayTable::Entry& h02 = myDecays.at(std::pair<int,int>(35,0));
const DecayTable::Entry& A0 = myDecays.at(std::pair<int,int>(36,0));
const DecayTable::Entry& Hpm = myDecays.at(std::pair<int,int>(37,0));
for (unsigned int i = 0; i < neutral_channels.size(); i++)
{
const std::vector<str>& chan = neutral_channels[i];
mssmwidths->hdwidth(i+1,2) = (h01.has_channel(chan) ? widths->width(DSparticle_code("h0_1")) * h01.BF(chan) : 0.0);
mssmwidths->hdwidth(i+1,1) = (h02.has_channel(chan) ? widths->width(DSparticle_code("h0_2")) * h02.BF(chan) : 0.0);
mssmwidths->hdwidth(i+1,3) = (A0.has_channel(chan) ? widths->width(DSparticle_code("A0")) * A0.BF(chan) : 0.0);
if (neutral_channels[i] == sister_chan)
{
// Add the missing W-H+ contributions.
mssmwidths->hdwidth(i+1,2) = (h01.has_channel(missing_chan) ? widths->width(DSparticle_code("h0_1")) * h01.BF(missing_chan) : 0.0);
mssmwidths->hdwidth(i+1,1) = (h02.has_channel(missing_chan) ? widths->width(DSparticle_code("h0_2")) * h02.BF(missing_chan) : 0.0);
mssmwidths->hdwidth(i+1,3) = (A0.has_channel(missing_chan) ? widths->width(DSparticle_code("A0")) * A0.BF(missing_chan) : 0.0);
}
}
for (unsigned int i = 0; i < charged_channels.size(); i++)
{
mssmwidths->hdwidth(i+1,4) = (Hpm.has_channel(charged_channels[i]) ? widths->width(DSparticle_code("H+")) * Hpm.BF(charged_channels[i]) : 0.0);
}
// Set up SM fermion widths
widths->width(DSparticle_code("u_3")) = myDecays.at(std::pair<int,int>(6,1)).width_in_GeV;
widths->width(DSparticle_code("e-_2")) = myDecays.at(std::pair<int,int>(13,1)).width_in_GeV;
widths->width(DSparticle_code("e-_3")) = myDecays.at(std::pair<int,int>(15,1)).width_in_GeV;
// Set up SM gauge boson widths
widths->width(DSparticle_code("W+")) = myDecays.at(std::pair<int,int>(24,0)).width_in_GeV;
widths->width(DSparticle_code("Z0")) = myDecays.at(std::pair<int,int>(23,0)).width_in_GeV;
// Set up sfermion widths
widths->width(DSpart->ksnu(1)) = myDecays.at(std::pair<int,int>(1000012,0)).width_in_GeV;
widths->width(DSpart->ksnu(2)) = myDecays.at(std::pair<int,int>(1000014,0)).width_in_GeV;
widths->width(DSpart->ksnu(3)) = myDecays.at(std::pair<int,int>(1000016,0)).width_in_GeV;
widths->width(DSpart->ksl(1)) = myDecays.at(std::pair<int,int>(1000011,0)).width_in_GeV;
widths->width(DSpart->ksl(2)) = myDecays.at(std::pair<int,int>(1000013,0)).width_in_GeV;
widths->width(DSpart->ksl(3)) = myDecays.at(std::pair<int,int>(1000015,0)).width_in_GeV;
widths->width(DSpart->ksl(4)) = myDecays.at(std::pair<int,int>(2000011,0)).width_in_GeV;
widths->width(DSpart->ksl(5)) = myDecays.at(std::pair<int,int>(2000013,0)).width_in_GeV;
widths->width(DSpart->ksl(6)) = myDecays.at(std::pair<int,int>(2000015,0)).width_in_GeV;
widths->width(DSpart->ksqu(1)) = myDecays.at(std::pair<int,int>(1000002,0)).width_in_GeV;
widths->width(DSpart->ksqu(2)) = myDecays.at(std::pair<int,int>(1000004,0)).width_in_GeV;
widths->width(DSpart->ksqu(3)) = myDecays.at(std::pair<int,int>(1000006,0)).width_in_GeV;
widths->width(DSpart->ksqu(4)) = myDecays.at(std::pair<int,int>(2000002,0)).width_in_GeV;
widths->width(DSpart->ksqu(5)) = myDecays.at(std::pair<int,int>(2000004,0)).width_in_GeV;
widths->width(DSpart->ksqu(6)) = myDecays.at(std::pair<int,int>(2000006,0)).width_in_GeV;
widths->width(DSpart->ksqd(1)) = myDecays.at(std::pair<int,int>(1000001,0)).width_in_GeV;
widths->width(DSpart->ksqd(2)) = myDecays.at(std::pair<int,int>(1000003,0)).width_in_GeV;
widths->width(DSpart->ksqd(3)) = myDecays.at(std::pair<int,int>(1000005,0)).width_in_GeV;
widths->width(DSpart->ksqd(4)) = myDecays.at(std::pair<int,int>(2000001,0)).width_in_GeV;
widths->width(DSpart->ksqd(5)) = myDecays.at(std::pair<int,int>(2000003,0)).width_in_GeV;
widths->width(DSpart->ksqd(6)) = myDecays.at(std::pair<int,int>(2000005,0)).width_in_GeV;
// Set up neutralino widths. Note that the zero neutralino width is taken care of below.
widths->width(DSpart->kn(1)) = myDecays.at(std::pair<int,int>(1000022,0)).width_in_GeV;
widths->width(DSpart->kn(2)) = myDecays.at(std::pair<int,int>(1000023,0)).width_in_GeV;
widths->width(DSpart->kn(3)) = myDecays.at(std::pair<int,int>(1000025,0)).width_in_GeV;
widths->width(DSpart->kn(4)) = myDecays.at(std::pair<int,int>(1000035,0)).width_in_GeV;
// Set up chargino widths.
widths->width(DSpart->kcha(1)) = myDecays.at(std::pair<int,int>(1000024,0)).width_in_GeV;
widths->width(DSpart->kcha(2)) = myDecays.at(std::pair<int,int>(1000037,0)).width_in_GeV;
// Gluino width.
widths->width(DSparticle_code("~g")) = myDecays.at(std::pair<int,int>(1000021,0)).width_in_GeV;
// Gravitino width (not implemented in DS).
//widths->width(DSparticle_code("~G")) = ;
// Integration routines in DS cannot handle very small sparticle widths.
// Make sure not to fall below minimal value in order to avoid numerical issues.
for (std::size_t i=21; i<49; i++)
{
if (widths->width(i)<min_DS_rwidth *mspctm->mass(i))
{
widths->width(i)=min_DS_rwidth *mspctm->mass(i);
}
}
#ifdef DARKSUSY_DEBUG
// Spit out spectrum and width files for debug purposes.
int u1 = 50;
int u2 = 100050;
#ifdef WITH_MPI
int rank = GMPI::Comm().Get_rank();
u1 += rank;
u2 += rank;
#endif
dswspectrum(u1);
dswwidth(u2);
#endif
return 0; // everything OK (hah. maybe.)
}
/// Returns direct detection couplings gps,gns,gpa,gna
/// (proton/neutron scalar/axial four-couplings)
std::vector<double> DD_couplings()
{
double gps,gns,gpa,gna;
dsddgpgn(gps,gns,gpa,gna);
std::vector<double> couplings;
couplings.clear();
couplings.push_back(gps);
couplings.push_back(gns);
couplings.push_back(gpa);
couplings.push_back(gna);
return couplings;
}
/// Returns the vector of neutral Higgs decay channels in DarkSUSY
std::vector< std::vector<str> > DS_neutral_h_decay_channels()
{
return initVector< std::vector<str> >
(initVector<str>("h0_2", "h0_2"),
initVector<str>("h0_1", "h0_2"),
initVector<str>("h0_1", "h0_1"),
initVector<str>("A0", "A0"),
initVector<str>("h0_2", "A0"),
initVector<str>("h0_1", "A0"),
initVector<str>("H+", "H-"),
initVector<str>("Z0", "h0_2"),
initVector<str>("Z0", "h0_1"),
initVector<str>("Z0", "A0"),
// actually supposed to be W+H- and W-H+
initVector<str>("W+", "H-"),
initVector<str>("Z0", "Z0"),
initVector<str>("W+", "W-"),
initVector<str>("nu_e", "nubar_e"),
initVector<str>("e+_1", "e-_1"),
initVector<str>("nu_mu", "nubar_mu"),
initVector<str>("e+_2", "e-_2"),
initVector<str>("nu_tau", "nubar_tau"),
initVector<str>("e+_3", "e-_3"),
initVector<str>("u_1", "ubar_1"),
initVector<str>("d_1", "dbar_1"),
initVector<str>("u_2", "ubar_2"),
initVector<str>("d_2", "dbar_2"),
initVector<str>("u_3", "ubar_3"),
initVector<str>("d_3", "dbar_3"),
initVector<str>("g", "g"),
// actually qqg (not implemented in DS though)
initVector<str>("d_3", "dbar_3", "g"),
initVector<str>("gamma", "gamma"),
initVector<str>("Z0", "gamma")
);
}
/// Returns the vector of charged Higgs decay channels in DarkSUSY
std::vector< std::vector<str> > DS_charged_h_decay_channels()
{
return initVector< std::vector<str> >
(initVector<str>("u_1", "dbar_1"),
initVector<str>("u_1", "dbar_2"),
initVector<str>("u_1", "dbar_3"),
initVector<str>("u_2", "dbar_1"),
initVector<str>("u_2", "dbar_2"),
initVector<str>("u_2", "dbar_3"),
initVector<str>("u_3", "dbar_1"),
initVector<str>("u_3", "dbar_2"),
initVector<str>("u_3", "dbar_3"),
initVector<str>("e+_1", "nu_e"),
initVector<str>("e+_2", "nu_mu"),
initVector<str>("e+_3", "nu_tau"),
initVector<str>("W+", "h0_2"),
initVector<str>("W+", "h0_1"),
initVector<str>("W+", "A0")
);
}
}
END_BE_NAMESPACE
Updated on 2024-07-18 at 13:53:36 +0000