Code calculates three-body decay modes and prints out an SLHA format file with them in. For R-parity conserving NMSSM/MSSM. See arXiv://1703.09717. More...

#include "nmssmsoftsusy.h"
#include "softsusy.h"
#include "physpars.h"
#include "lowe.h"
#include "def.h"
#include "softpars.h"
#include "flavoursoft.h"
#include "susy.h"
#include "particle.h"
#include "twoBodyDecays.h"
#include <iostream>
#include <cstring>
#include <stdlib.h>
#include <cmath>
#include <vector>
#include <iomanip>
#include <complex>

Functions
void	printRowPDG (ostream &cout, double x)

void	OutputNoPWs (ostream &cout, Particle &P)
	Outputs the decay table into the leshouchesOutput file with no PWs (partial widths) given, just branching ratios.

void	OutputYesPWs (ostream &cout, Particle &P)
	Outputs the decay table into the leshouchesOutput file with PWs (partial widths) given after the comments column so as not to affect SLHA form.

double	Zsfintegralsum (double m1, double m2, double msf, double mf, double mz, double min, double max, double Nsteps, int adaptive, double approx)
	Integral calculating functions for 1->3 decays.

double	Aintegralsum (double m1, double m2, double mz, double mA, double mf, double min, double max, double Nsteps, int adaptive, double approx)

double	G3integralsum (double m1, double m2, double m3, double mt, double mb, double min, double max, double Nsteps, int adaptive, double approx)

double	G2integralsum (double m1, double m2, double m3, double mt, double mb, double min, double max, double Nsteps, int adaptive, double approx)

double	chiprimeintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	Yintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	Xintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	zetaintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	chiintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	rhointegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	xsiintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	phitildaintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	psitildaintegralsum (double m1, double m2, double m3, double m4, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	Zintegralsum (double m1, double m2, double mz, double mf, double min, double max, double Nsteps, int adaptive, double approx)

double	G1integralsum (double m1, double m2, double m3, double mq, double min, double max, double Nsteps, int adaptive, double approx)

double	G4integralsum (double m1, double m2, double m3, double m4, double mt, double mb, double min, double max, double Nsteps, int adaptive, double approx)

double	G5integralsum (double m1, double m2, double m3, double m4, double mt, double mb, double min, double max, double Nsteps, int adaptive, double approx)

double	G6integralsum (double m1, double m2, double m3, double m4, double mt, double mb, double min, double max, double Nsteps, int adaptive, double approx)

double	G7integralsum (double m1, double m2, double m3, double m4, double mt, double mb, double min, double max, double Nsteps, int adaptive, double approx)

double	G8integralsum (double m1, double m2, double m3, double m4, double mt, double mb, double min, double max, double Nsteps, int adaptive, double approx)

double	Jintegralsum (double m1, double m2, double msf, double mphi, double mf, double min, double max, double Nsteps, double approx, int adaptive, int AorhorH)

double	ghHintegral (double m1, double m2, double mf, double mh, double mH, double g, double gp, double alpha, int neutralinoi, int neutralinoj, DoubleMatrix &mixNeut, double E)

double	hHintegral (double m1, double m2, double mf, double mh, double mH, double fromE, double toE, double stepE, double g, double gp, double alpha, DoubleMatrix &mixNeut, int neutralinoi, int neutralinoj)

double	gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogen (double mgluino, double mneutralino, double msqL, double msqR, double mquark, double g, double gp, DoubleMatrix &mixNeut, double alphas, char uord, int neut, bool onetothree)
	Hintegral. More...

double	gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogenlimit (double mgluino, double mneutralino, double msqL, double msqR, double mquark, double g, double gp, DoubleMatrix &mixNeut, double alphas, char uord, int neut, bool onetothree)
	m1 is mgluino, m2 is neutralinoi mass, m3 is sqL mass, m4 is sqR mass, m5 is quark mass but assumed zero in calculation here, just used to check allowed for now; char uord tells us if the quark is u type 'u' or d type 'd', int neut tells us which neutralino it is More...

double	gluinoamplitudedecaydgaussneutralinottbar (double mgluino, double mst1, double mst2, double mneutralino, double mt, double mWboson, double g, double gp, double thetat, double beta, double alphas, DoubleMatrix &mixNeut, double runmt, int neutralino, bool onetothree, char torb)
	calculates PW for gluino -> neutralino + q qbar pair where q are t More...

double	gluinoamplitudedecaydgausschartbbar (double mgluino, double mst1, double mst2, double msb1, double msb2, double mtop, double mbottom, double mchar, double alphas, double thetat, double thetab, double MWBoson, double g, double gp, double gammaL, double gammaR, double beta, double runmt, double runmb, int chargino, bool onetothree)

double	neutralinoamplitudedecaydgaussneutralinoffbar (double mneutralinoi, double msf1, double msf2, double mZboson, double mhiggsl, double mhiggsH, double mhiggsA, double mneutralinoj, double mf, double alphas, double thetaq, double mWboson, double g, double gp, double alpha, double beta, double runmq, DoubleMatrix &mixNeut, int ineutralino, int jneutralino, bool onetothree, char uordornuorl)

double	neutralinoamplitudedecaycharginoffprimebar (double mneutralinoi, double msfp1, double msfp2, double msf1, double msf2, double mWboson, double mHP, double mcharginoj, double mfp, double mf, double thetaq, double thetaqp, double g, double gp, double alpha, double beta, double thetaL2, double thetaR2, double runmqp, double runmq, DoubleMatrix &mixNeut, int ineutralino, int jchargino, bool onetothree, char qorl, char norc)

double	fdgauss (double x)

double	gpsitildadgauss (double Et)

double	gphitildadgauss (double Et)

double	gxsidgauss (double Et)

double	grhodgauss (double Et)

double	gchidgauss (double Et)

double	gzetadgauss (double Et)

double	gXdgauss (double Et)

double	gYdgauss (double Et)

double	gchiprimedgauss (double Et)

double	gG1dgauss (double Et)
	m1 = mgluino, m2 = mstopi, m6 = mtop, m8 = mcharginoj

double	gG4dgauss (double Et)
	m1 = mgluino, m2 = mstopi, m4 = msbottomi, m6 = mtop, m7 = mbottom , m8 = mcharginoj

double	gG5dgauss (double Et)
	m1 = mgluino, m2 = mstopi, m4 = msbottomi, m6 = mtop, m7 = mbottom , m8 = mcharginoj

double	gG6dgauss (double Et)
	m1 = mgluino, m2 = mstopi, m4 = msbottomi, m6 = mtop, m7 = mbottom , m8 = mcharginoj

double	gG7dgauss (double Et)

double	gG8dgauss (double Et)

double	gG2dgauss (double Eb)

double	gG3dgauss (double Eb)

double	gZdgauss (double E)
	m1 = mZi, m4 = mZj, m2 = mstop1, m3 = mstop2, mq = mt, MZboson = mz

double	ghHdgauss (double E)
	m1 = mZi, m4 = mZj, m2 = mstop1, m3 = mstop2, mq = mt, MZboson = mz

double	gAdgauss (double E)
	m1 = mZi, m4 = mZj, m2 = mstop1, m3 = mstop2, mq = mt

double	gZsfdgauss (double s)
	m2 = msfi where i is the sfermion contribution interference with Z you are considering

double	gJdgauss (double s)

double	gneutineutjffZ1dgauss (double s)
	m1 = mneuti, m4 = mneutj, mq = mf, MZboson = mZboson More...

double	gneutineutjffZ2dgauss (double s)

double	gneutineutjffZ3dgauss (double s)

double	gneutineutjffZ4dgauss (double s)

double	gintegralhdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl

double	gintegralHdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mH = mhiggsH

double	gintegralh1dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl

double	gintegralh2dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl

double	gintegralh3dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl

double	gintegralh4dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl

double	gintegralH1dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mH = mhiggsH

double	gintegralH2dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mH = mhiggsH

double	gintegralH3dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mH = mhiggsH

double	gintegralH4dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl

double	gintegralhH1dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl, mH = mhiggsH

double	gintegralhH2dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl, mH = mhiggsH

double	gintegralhH3dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl, mH = mhiggsH

double	gintegralhH4dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mh = mhiggsl, mH = mhiggsH

double	gintegralA1dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mA = mhiggsA

double	gintegralA2dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mA = mhiggsA

double	gintegralA3dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mA = mhiggsA

double	gintegralA4dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mt, MZboson = mz, mA = mhiggsA

double	gintegral1Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral2Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral3Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral4Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral5Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral6Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral7Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral8Zsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral1hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral2hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral3hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral4hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral5hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral6hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral7hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral8hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral1Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral2Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral3Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral4Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral5Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral6Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral7Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral8Hsfdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2) More...

double	gintegral1ZAdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, m2 = msfi, MZboson = mz

double	gintegral2ZAdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, m2 = msfi, MZboson = mz

double	gintegral3ZAdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, m2 = msfi, MZboson = mz

double	gintegral4ZAdgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, m2 = msfi, MZboson = mz

double	gneutineutjffgA1dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, MZboson = mz

double	gneutineutjffgA2dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, MZboson = mz

double	gneutineutjffgA3dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, MZboson = mz

double	gneutineutjffgA4dgauss (double E)
	m1 = mZi, m4 = mZj, mq = mf, mA = mhiggsA, MZboson = mz

double	gneuticharjffpW1dgauss (double E)
	m1 = mZi, m2 = mWj, m3 = mf, m4 = mfp, MWboson = mw

double	gneuticharjffpW2dgauss (double E)
	m1 = mZi, m2 = mWj, m3 = mf, m4 = mfp, MWboson = mw

double	gneuticharjffpHpm1dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mHpm; //sometimes mass order differs as called in different cases

double	gneuticharjffpHpm2dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mHpm; //sometimes mass order differs as called in different cases

double	gneuticharjffpHpm3dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mHpm; //sometimes mass order differs as called in different cases

double	gneuticharjffpHpm4dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mHpm; //sometimes mass order differs as called in different cases

double	gneuticharjffp1sf1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2 More...

double	gneuticharjffp2sf1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2 More...

double	gneuticharjffp3sf1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2 More...

double	gneuticharjffp4sf1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2 More...

double	gneuticharjffp1sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp2sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp3sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp4sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp5sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp6sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp7sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp8sfp1sf2dgauss (double E)
	m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msfp1, m6 = msf2

double	gneuticharjffp1WHpmdgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHpm

double	gneuticharjffp2WHpmdgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHpm

double	gneuticharjffp3WHpmdgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHpm

double	gneuticharjffp4WHpmdgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHpm

double	gneuticharjffpW1Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpW2Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpW3Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpW4Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpW5Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpW6Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpW7Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpW8Sfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp

double	gneuticharjffpHg1dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHP

double	gneuticharjffpHg2dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHP

double	gneuticharjffpHg3dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHP

double	gneuticharjffpHg4dgauss (double E)
	m1 = mneutralinoi, m2 = mcharginoj, m3 = mf, m4 = mfp, m5 = mWboson, m6 = mHP

double	gneuticharjffp1gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp2gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp3gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp4gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp5gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp6gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp7gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp8gsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = mf, m4 = mcharginoj, m5 = mWboson, m6 = msfp1

double	gneuticharjffp1sfpsfpdgauss (double E)
	m1 = mneutralinoi, m2 = mfp, m3 = -mf, m4 = charginoj, m5 = msfp1, m6 = msfp2

double	gneuticharjffp21sfpsfpdgauss (double E)

double	gneuticharjffp31sfpsfpdgauss (double E)

double	gneuticharjffp41sfpsfpdgauss (double E)

DoubleVector	squarkmixcharginocouplings (double g, double theta, double beta, double gammaL, double gammaR, double runmt, double runmb, double mWboson, int torb)

DoubleVector	higgsphisamecharginocouplings (double alpha, double beta, double thetaL, double thetaR)
	calculates the couplings of a neutral higgs (h,H,A) to Wtildai Wtildai - i.e. to two charginos of the same type (mass)

DoubleVector	higgsphidifcharginocouplings (double alpha, double beta, double thetaL, double thetaR)
	calculates the couplings of a netural higgs (h,H,A) to Wtildai Wtildaj - i.e. to two charginos of different type (mass)

DoubleVector	higgshsquarksamehandcouplings (double mWboson, double g, double gp, double alpha, double beta, double mupq, double mdownq)
	calculates the couplings of light scalar higgs h to two squarks of same handedness

DoubleVector	higgshsquarkdiffhandcouplings (double mWboson, double g, double alpha, double beta, double mupq, double mdownq, double greekmu, double Aup, double Adown)
	calculates the couplings of light scalar higgs h to two squarks of different handedness

DoubleVector	higgsHsquarksamehandcouplings (double mWboson, double g, double gp, double alpha, double beta, double mupq, double mdownq)
	calculates the couplings of heavy scalar higgs H to two squarks of same handedness

DoubleVector	higgsHsquarkdiffhandcouplings (double mWboson, double g, double alpha, double beta, double mupq, double mdownq, double greekmu, double Aup, double Adown)
	calculates the couplings of heavy scalar higgs H to two squarks of different handedness

DoubleVector	higgshsleptonsamehandcouplings (double mWboson, double g, double gp, double alpha, double beta, double mel)
	calculates the couplings of light scalar higgs h to two sleptons of same handedness

DoubleVector	higgshsleptondiffhandcouplings (double mWboson, double g, double alpha, double beta, double mel, double greekmu, double Ae)
	calculates the couplings of light scalar higgs h to two sleptons of different handedness

DoubleVector	higgsHsleptonsamehandcouplings (double mWboson, double g, double gp, double alpha, double beta, double mel)
	calculates the couplings of heavy scalar higgs H to two sleptons of same handedness

DoubleVector	higgsHsleptondiffhandcouplings (double mWboson, double g, double alpha, double beta, double mel, double greekmu, double Ae)
	calculates the couplings of heavy scalar higgs H to two sleptons of different handedness

DoubleVector	higgsHplussquarkcouplings (double mWboson, double g, double beta, double mup, double mdown, double greekmu, double Au, double Ad)
	calculates the couplings of charged scalar higgs H+ to two squarks

DoubleVector	squarkmixcharginocouplings (double g, double theta, double beta, double gammaL, double gammaR, double runmt, double runmb, double mWboson, double mch1, double mch2, int torb)
	Functions that calculate the couplings:

DoubleVector	foftau (double mpart, double mcomp)
	Function Declarations used to calculate Partial Widths in decays.cpp. More...

DoubleVector	goftau (double mpart, double mcomp)
	g(tau) function for use in h->Z gamma

double	chToN2piInt (double qSq, const DoubleVector &v)

double	charginoToNeutralino2pion (const MssmSoftsusy *m)

Complex	fofqsq (double qSq)

Complex	bw (double mSq, double gamma, double qSq)

Variables
const double	accuracy = 0.01
	Approximate accuracy with which 3 body decays are calculated.

const double	mcpole = 1.40

const double	mspole = 0.19

const double	hadronicScale = 1.5

Detailed Description

Code calculates three-body decay modes and prints out an SLHA format file with them in. For R-parity conserving NMSSM/MSSM. See arXiv://1703.09717.

Project: SOFTSUSY
Author: Tom Cridge, Ben Allanach
Manual: hep-ph/0104145, Comp. Phys. Comm. 143 (2002) 305
Webpage: http://hepforge.cedar.ac.uk/softsusy/

Function Documentation

◆ fdgauss()

double fdgauss ( double x )

dgauss integrands for numerical integration of phase space in 1->3 decays Make sure you have set useInDGauss before calling this function

◆ foftau()

DoubleVector foftau	(	double	mpart,
		double	mcomp
	)

Function Declarations used to calculate Partial Widths in decays.cpp.

f(tau) function for use in h->gamma gamma or Z gamma

◆ gintegral1hsfdgauss()

double gintegral1hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral1Hsfdgauss()

double gintegral1Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral1Zsfdgauss()

double gintegral1Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral2hsfdgauss()

double gintegral2hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral2Hsfdgauss()

double gintegral2Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral2Zsfdgauss()

double gintegral2Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral3hsfdgauss()

double gintegral3hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral3Hsfdgauss()

double gintegral3Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral3Zsfdgauss()

double gintegral3Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral4hsfdgauss()

double gintegral4hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral4Hsfdgauss()

double gintegral4Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral4Zsfdgauss()

double gintegral4Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral5hsfdgauss()

double gintegral5hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral5Hsfdgauss()

double gintegral5Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral5Zsfdgauss()

double gintegral5Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral6hsfdgauss()

double gintegral6hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral6Hsfdgauss()

double gintegral6Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral6Zsfdgauss()

double gintegral6Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral7hsfdgauss()

double gintegral7hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral7Hsfdgauss()

double gintegral7Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral7Zsfdgauss()

double gintegral7Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral8hsfdgauss()

double gintegral8hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral8Hsfdgauss()

double gintegral8Hsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, mhiggsl = mh, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gintegral8Zsfdgauss()

double gintegral8Zsfdgauss ( double E )

m1 = mZi, m4 = mZj, mq = mf, MZboson = mz, m2 = msfi where i is the sfermion index you're considering for that interference (i = 1,2)

argument of the log

◆ gluinoamplitudedecaydgausschartbbar()

double gluinoamplitudedecaydgausschartbbar	(	double	mgluino,
		double	mst1,
		double	mst2,
		double	msb1,
		double	msb2,
		double	mtop,
		double	mbottom,
		double	mchar,
		double	alphas,
		double	thetat,
		double	thetab,
		double	MWBoson,
		double	g,
		double	gp,
		double	gammaL,
		double	gammaR,
		double	beta,
		double	runmt,
		double	runmb,
		int	chargino,
		bool	onetothree
	)

Different to T&B, they have + sin(2*thetat)*G8, I follow SPheno

Different to T&B, they have - sin(2*thetat)*G8, I follow SPheno

Global minus sign difference cf T&B, I follow SPheno

Different to T&B, they have - alphasb1ch*betasb1ch*G3sb1, I follow SPheno

Different to T&B, they have - alphasb2ch*betasb2ch*G3sb2, I follow SPheno

Taken mb -> 0 limit in squared matrix element (but not in phase space) so no sb1sb2 interference term

The sign in front of the G7st1sb1 and G5st1sb1 terms has been changed from that given in Baer and Tata in order to agree with SPheno.

◆ gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogen()

double gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogen	(	double	mgluino,
		double	mneutralino,
		double	msqL,
		double	msqR,
		double	mquark,
		double	g,
		double	gp,
		DoubleMatrix &	mixNeut,
		double	alphas,
		char	uord,
		int	neut,
		bool	onetothree
	)

Hintegral.

1 to 3 decay functions calling dgauss

Hintegral.

1->3 decay not relevant here if one to two decay open

◆ gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogenlimit()

double gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogenlimit	(	double	mgluino,
		double	mneutralino,
		double	msqL,
		double	msqR,
		double	mquark,
		double	g,
		double	gp,
		DoubleMatrix &	mixNeut,
		double	alphas,
		char	uord,
		int	neut,
		bool	onetothree
	)

m1 is mgluino, m2 is neutralinoi mass, m3 is sqL mass, m4 is sqR mass, m5 is quark mass but assumed zero in calculation here, just used to check allowed for now; char uord tells us if the quark is u type 'u' or d type 'd', int neut tells us which neutralino it is

1->3 decay not relevant here if one to two decay open

◆ gluinoamplitudedecaydgaussneutralinottbar()

double gluinoamplitudedecaydgaussneutralinottbar	(	double	mgluino,
		double	mst1,
		double	mst2,
		double	mneutralino,
		double	mt,
		double	mWboson,
		double	g,
		double	gp,
		double	thetat,
		double	beta,
		double	alphas,
		DoubleMatrix &	mixNeut,
		double	runmt,
		int	neutralino,
		bool	onetothree,
		char	torb
	)

calculates PW for gluino -> neutralino + q qbar pair where q are t

Note the effect of the complex couplings and the pm factors for negative masses cancel out as the couplings always appear in pairs so multiplying them gives an extra minus sign when they are purely imaginary, but this extra minus sign is cancelled out by the extra minus sign in the pm factor. Therefore really the additional minus signs come from the fact I've used the neutralino mass itself throughout my calculation (rather than its absolute value) which therefore naturally introduces additional minus signs.

Gammast1st2 terms as from T&B but with sign change for X term and also the factor of -(cos^2(thetat)-sin^2(thetat))(st1alpha1*st2alpha1+st1beta1*st2beta1) T&B has in the chiprime term has been changed to (cos^2(thetat)+sin^2(thetat))(st1alpha1*st2alpha1 - st1beta1*st2beta1) = 1*(st1alpha1*st2alpha1 - st1beta1*st2beta1), in order in both cases to agree with SPHENO

Extra terms SPHENO has in Gammast1st2 not present in T&B

◆ gneuticharjffp1sf1sf2dgauss()

double gneuticharjffp1sf1sf2dgauss ( double E )

m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2

Avoid small negative values of squareplus at upper boundary where theoretically s = 0 exactly.

◆ gneuticharjffp2sf1sf2dgauss()

double gneuticharjffp2sf1sf2dgauss ( double E )

m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2

Avoid small negative values of squareplus at upper boundary where theoretically s = 0 exactly.

◆ gneuticharjffp3sf1sf2dgauss()

double gneuticharjffp3sf1sf2dgauss ( double E )

m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2

Avoid small negative values of squareplus at upper boundary where theoretically s = 0 exactly.

◆ gneuticharjffp4sf1sf2dgauss()

double gneuticharjffp4sf1sf2dgauss ( double E )

m1 = mneutralinoi, m2 = mf, m3 = mcharginoj, m4 = mfp, m5 = msf1, m6 = msf2

Avoid small negative values of squareplus at upper boundary where theoretically s = 0 exactly.

◆ gneutineutjffZ1dgauss()

double gneutineutjffZ1dgauss ( double s )

m1 = mneuti, m4 = mneutj, mq = mf, MZboson = mZboson

Set to 0 to avoid numerical precision causing a very small but positive squareminus1 at smax when it should be exactly 0, this can cause problems as it would make lambda1 the sqrt of a negative number (as squareplus1 is negative), avoid this by setting to 0

◆ gneutineutjffZ2dgauss()

double gneutineutjffZ2dgauss ( double s )

Set to 0 to avoid numerical precision causing a very small but positive squareminus1 at smax when it should be exactly 0, this can cause problems as it would make lambda1 the sqrt of a negative number (as squareplus1 is negative), avoid this by setting to 0

◆ gneutineutjffZ3dgauss()

double gneutineutjffZ3dgauss ( double s )

Set to 0 to avoid numerical precision causing a very small but positive squareminus1 at smax when it should be exactly 0, this can cause problems as it would make lambda1 the sqrt of a negative number (as squareplus1 is negative), avoid this by setting to 0

◆ gneutineutjffZ4dgauss()

double gneutineutjffZ4dgauss ( double s )

Set to 0 to avoid numerical precision causing a very small but positive squareminus1 at smax when it should be exactly 0, this can cause problems as it would make lambda1 the sqrt of a negative number (as squareplus1 is negative), avoid this by setting to 0

◆ gpsitildadgauss()

double gpsitildadgauss ( double Et )

< this can happen erronesouly at very end of range due to finite precision used, squareplus should actually then be very very small and +ve

◆ neutralinoamplitudedecaycharginoffprimebar()

double neutralinoamplitudedecaycharginoffprimebar	(	double	mneutralinoi,
		double	msfp1,
		double	msfp2,
		double	msf1,
		double	msf2,
		double	mWboson,
		double	mHP,
		double	mcharginoj,
		double	mfp,
		double	mf,
		double	thetaq,
		double	thetaqp,
		double	g,
		double	gp,
		double	alpha,
		double	beta,
		double	thetaL2,
		double	thetaR2,
		double	runmqp,
		double	runmq,
		DoubleMatrix &	mixNeut,
		int	ineutralino,
		int	jchargino,
		bool	onetothree,
		char	qorl,
		char	norc
	)

Just usual yukawa for up type

Just usual yukawa for down type

correction factor for negative masses

correction factor for if it's chargino -> neutralino fp fbar rather than neutralino -> chargino fpbar f

No snuR exists so only one sfp contribution

GammaW

Hpm1 contribution (W+ goldstone):

Hpm2 contribution (Actual Hpm):

Sfp Sfp diagonal (Remember fp is u-type fermion)

Sf Sf diagonal

Sfp t sf u interference Sfp 1 sf 2

Sfp 1 Sf 1 interference

Sfp 2 Sf 2 interference

Sfp 2 Sf 1 interference

W-Hpm interference

W-goldstone interference

W Sfp 1 interference

W Sfp 2 interference

W Sf1 interference

W Sf2 interference

H+ goldstone interference

goldstone - sfp1 interference

swap which integral goes with which couplingcombo

goldstone - sfp2 interference

H+ - sfp1 interference

H+ - sfp2 interference

goldstone - sf1 interference

goldstone - sf2 interference

Hpm - sf2 interference

Sfp Sfp interference

Sf - Sf interference

◆ neutralinoamplitudedecaydgaussneutralinoffbar()

double neutralinoamplitudedecaydgaussneutralinoffbar	(	double	mneutralinoi,
		double	msf1,
		double	msf2,
		double	mZboson,
		double	mhiggsl,
		double	mhiggsH,
		double	mhiggsA,
		double	mneutralinoj,
		double	mf,
		double	alphas,
		double	thetaq,
		double	mWboson,
		double	g,
		double	gp,
		double	alpha,
		double	beta,
		double	runmq,
		DoubleMatrix &	mixNeut,
		int	ineutralino,
		int	jneutralino,
		bool	onetothree,
		char	uordornuorl
	)

correction factor for negative masses

First sf t sf u msf1-msf1, msf1-msf2, msf2-msf2 components:

the problem is in here - DEBUG

sf sf diag contribution 2x as t and u same

sf sf non diag 2x as t and u same

Pseudoscalar A contribution:

Have an additional pseudoscalar goldstone with mass of Z due to being in Feynman gauge, to add in the longitudinal polarisation of the Z

set so can use same integrals as for pseudoscalar A boson but now with Z mass

here set mh = mA as using the integrands written for h with mh -> mA

Goldstone - sf interference:

here set mh = mZ as using the integrands written for h with mh -> mA -> mZ

Goldstone-A interference

Coupling combos:

◆ printRowPDG()

void printRowPDG	(	ostream &	cout,
		double	x
	)

make it return a character when you've worked out the equivalent of printf

Traps -0.0

Variable Documentation

◆ hadronicScale

const double hadronicScale = 1.5

When mass splitting is below this scale, chargino decays involving hadrons are calculated rather than to quarks

◆ mcpole

const double mcpole = 1.40

Change mcpole to alter quark "pole" masses used in h -> qq with QCD corrections

◆ mspole

const double mspole = 0.19

Change mspole to alter quark "pole" masses used in h -> qq with QCD corrections

Functions

Variables

Detailed Description

Function Documentation

◆ fdgauss()

◆ foftau()

◆ gintegral1hsfdgauss()

◆ gintegral1Hsfdgauss()

◆ gintegral1Zsfdgauss()

◆ gintegral2hsfdgauss()

◆ gintegral2Hsfdgauss()

◆ gintegral2Zsfdgauss()

◆ gintegral3hsfdgauss()

◆ gintegral3Hsfdgauss()

◆ gintegral3Zsfdgauss()

◆ gintegral4hsfdgauss()

◆ gintegral4Hsfdgauss()

◆ gintegral4Zsfdgauss()

◆ gintegral5hsfdgauss()

◆ gintegral5Hsfdgauss()

◆ gintegral5Zsfdgauss()

◆ gintegral6hsfdgauss()

◆ gintegral6Hsfdgauss()

◆ gintegral6Zsfdgauss()

◆ gintegral7hsfdgauss()

◆ gintegral7Hsfdgauss()

◆ gintegral7Zsfdgauss()

◆ gintegral8hsfdgauss()

◆ gintegral8Hsfdgauss()

◆ gintegral8Zsfdgauss()

◆ gluinoamplitudedecaydgausschartbbar()

◆ gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogen()

◆ gluinoamplitudedecaydgaussneutralinoqqpbarfirsttwogenlimit()

◆ gluinoamplitudedecaydgaussneutralinottbar()

◆ gneuticharjffp1sf1sf2dgauss()

◆ gneuticharjffp2sf1sf2dgauss()

◆ gneuticharjffp3sf1sf2dgauss()

◆ gneuticharjffp4sf1sf2dgauss()

◆ gneutineutjffZ1dgauss()

◆ gneutineutjffZ2dgauss()

◆ gneutineutjffZ3dgauss()

◆ gneutineutjffZ4dgauss()

◆ gpsitildadgauss()

◆ neutralinoamplitudedecaycharginoffprimebar()

◆ neutralinoamplitudedecaydgaussneutralinoffbar()

◆ printRowPDG()

Variable Documentation

◆ hadronicScale

◆ mcpole

◆ mspole