gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
Gambit::DarkBit Namespace Reference

Namespaces

 DarkBit_utils
 
 DecayChain
 

Classes

struct  nudata
 Neutrino telescope data container. More...
 
struct  nuyield_info
 Neutrino telescope yield info container. More...
 
struct  RD_coannihilating_particle
 
struct  RD_spectrum_type
 
struct  SimpleHist
 Histogram class for cascade decays. More...
 
struct  SimYieldChannel
 Annihilation/decay channel. More...
 
class  SimYieldTable
 Channel container Object containing tabularized yields for particle decay and two-body final states. More...
 
struct  TH_Channel
 All data on a single annihilation or decay channel, e.g. More...
 
struct  TH_ParticleProperty
 A container for the mass and spin of a particle. More...
 
struct  TH_Process
 A container for a single process. More...
 
struct  TH_ProcessCatalog
 A container holding all annihilation and decay initial states relevant for DarkBit. More...
 
struct  TH_Resonance
 A single resonance of a given width at a given energy (both in GeV) More...
 
struct  TH_resonances_thresholds
 Location of resonances and thresholds in energy (GeV) More...
 
struct  Wstruct
 

Typedefs

typedef std::map< str, std::map< str, Gambit::DarkBit::SimpleHist > > simpleHistContainter
 
typedef std::map< str, intstringIntMap
 
typedef std::map< str, daFunk::FunkstringFunkMap
 

Enumerations

enum  InterpolationOptions1D { InterpolationOptions1D::linear, InterpolationOptions1D::cspline }
 Generic one-dimensional integration container for linear interpolation and cubic splines. More...
 
enum  InterpolationOptions2D { InterpolationOptions2D::bilinear, InterpolationOptions2D::bicubic }
 Two-dimensional integration container for bilinear interpolation and bicubic splines. More...
 
enum  cascadeMC_SpecialEvents { MC_INIT =-1, MC_NEXT_STATE =-2, MC_FINALIZE =-3 }
 Special events for event loop. More...
 

Functions

void CMC_dummy (DarkBit::stringFunkMap &result)
 
void createSpectrum (Spectrum &outSpec)
 
void createDecays (DecayTable &outDecays)
 
void createSLHA1Names (mass_es_pseudonyms &names)
 
void TH_ProcessCatalog_WIMP (TH_ProcessCatalog &result)
 
void DarkMatter_ID_WIMP (std::string &result)
 
void DarkMatterConj_ID_WIMP (std::string &result)
 
void DD_couplings_WIMP (DM_nucleon_couplings &result)
 
errorDarkBit_error ()
 
warningDarkBit_warning ()
 
double parabola (double x, const double params[])
 H.E.S.S.-likelihood-related interpolation routines. More...
 
double intersect_parabola_line (double a, double b, double sign, const double pparams[])
 
double rho_integrand (double rho, void *params)
 
double rad_integrand (double rad, void *params)
 
double erg_integrand (double erg, void *params)
 
double alt_erg_integrand (double erg, void *params)
 
double gaussian_nuisance_lnL (double theo, double obs, double sigma)
 
double gStar (double T)
 Various capabilities and functions to provide SM physics as well as QCD input for axions. More...
 
double gStar_S (double T)
 
void QCDAxion_ZeroTemperatureMass_Nuisance_lnL (double &result)
 
void QCDAxion_AxionPhotonConstant_Nuisance_lnL (double &result)
 
double log_chi (double T, double beta, double Tchi)
 
void QCDAxion_TemperatureDependence_Nuisance_lnL (double &result)
 
double ALPS1_signal_general (double power, double nm1, double m_ax, double gagg)
 Likelihoods for ALPS 1 (LSW), CAST (helioscopes), and ADMX, UF, RBF (haloscopes). More...
 
void calc_ALPS1_signal_vac (double &result)
 
void calc_ALPS1_signal_gas (double &result)
 
double ALPS1_lnL_general (double s, double mu, double sigma)
 
void calc_lnL_ALPS1 (double &result)
 
void calc_CAST2007_signal_vac (std::vector< double > &result)
 
void calc_CAST2017_signal_vac (std::vector< std::vector< double >> &result)
 
double CAST_lnL_general (std::vector< double > s, const std::vector< double > bkg_counts, const std::vector< int > sig_counts)
 
void calc_lnL_CAST2007 (double &result)
 
void calc_lnL_CAST2017 (double &result)
 
void calc_Haloscope_signal (double &result)
 
void calc_lnL_Haloscope_ADMX1 (double &result)
 
void calc_lnL_Haloscope_ADMX2 (double &result)
 
void calc_lnL_Haloscope_UF (double &result)
 
void calc_lnL_Haloscope_RBF (double &result)
 
double SpecialFun1 (double T)
 Capabilities relating to axion cosmology. Currently only provides the energy density in axions today due to the realignment mechanism. More...
 
double SpecialFun3 (double T)
 
double hubble_rad_dom (double T)
 
double axion_mass_temp (double T, double beta, double Tchi)
 
double equation_Tosc (double T, void *params)
 
void calc_AxionOscillationTemperature (double &result)
 
int scal_field_eq (double tau, const double y[], double f[], void *params)
 
int scal_field_eq_jac (double tau, const double y[], double *dfdy, double dfdt[], void *params)
 
void RD_oh2_Axions (double &result)
 
void calc_RParameter (double &result)
 Capabilities relating to astrophysical observations (R-parameter, H.E.S.S. telescope search, cooling hints). More...
 
void calc_lnL_RParameter (double &result)
 
void calc_lnL_WDVar_G117B15A (double &result)
 
void calc_lnL_WDVar_R548 (double &result)
 
void calc_lnL_WDVar_PG1351489 (double &result)
 
void calc_lnL_WDVar_L192 (double &result)
 
void calc_lnL_SN1987A (double &result)
 
void calc_PhotonFluence_SN1987A_Conversion (double &result)
 
void calc_lnL_HESS_GCMF (double &result)
 
void calc_lnL_XENON1T_Anomaly (double &result)
 Capability for the XENON1T likelihood from 2006.10035. More...
 
double dRdE (double E, void *params)
 
void calc_lnL_XENON1T_DM_Anomaly (double &result)
 
void calc_lnL_XENON1T_Anomaly_NuisanceParameters (double &result)
 
void cascadeMC_FinalStates (std::vector< std::string > &list)
 Function for retrieving list of final states for cascade decays. More...
 
void cascadeMC_DecayTable (DarkBit::DecayChain::DecayTable &table)
 Function setting up the decay table used in decay chains. More...
 
void cascadeMC_LoopManager ()
 Loop manager for cascade decays. More...
 
void cascadeMC_InitialState (std::string &pID)
 Function selecting initial state for decay chain. More...
 
void cascadeMC_EventCount (std::map< std::string, int > &counts)
 Event counter for cascade decays. More...
 
void cascadeMC_GenerateChain (DarkBit::DecayChain::ChainContainer &chain)
 Function for generating decay chains. More...
 
void cascadeMC_sampleSimYield (const SimYieldTable &table, const DarkBit::DecayChain::ChainParticle *endpoint, std::string finalState, const TH_ProcessCatalog &catalog, std::map< std::string, std::map< std::string, SimpleHist > > &histList, std::string initialState, double weight, int cMC_numSpecSamples)
 Function for sampling SimYieldTables (tabulated spectra). More...
 
void cascadeMC_Histograms (std::map< std::string, std::map< std::string, SimpleHist > > &result)
 Function responsible for histogramming, and evaluating end conditions for event loop. More...
 
void cascadeMC_fetchSpectra (std::map< std::string, daFunk::Funk > &spectra, std::string finalState, const std::vector< std::string > &ini, const std::vector< std::string > &fin, const std::map< std::string, std::map< std::string, SimpleHist > > &h, const std::map< std::string, int > &eventCounts)
 Convenience function for getting a daFunk::Funk object of a given spectrum. More...
 
void cascadeMC_gammaSpectra (std::map< std::string, daFunk::Funk > &spectra)
 Function requesting and returning gamma ray spectra from cascade decays. More...
 
void mwimp_generic (double &result)
 Retrieve the DM mass in GeV for generic models (GeV) More...
 
void sigmav_late_universe (double &result)
 Retrieve the total thermally-averaged annihilation cross-section for indirect detection (cm^3 / s). More...
 
void sigmav_late_universe_MicrOmegas (double &result)
 
void DM_process_from_ProcessCatalog (std::string &result)
 Information about the nature of the DM process in question (i.e. More...
 
double profile_gNFW (double rhos, double rs, double alpha, double beta, double gamma, double r)
 Generalized NFW dark matter halo profile function. More...
 
double profile_Einasto (double rhos, double rs, double alpha, double r)
 Einasto dark matter halo profile function. More...
 
void GalacticHalo_gNFW (GalacticHaloProperties &result)
 Module function to generate GalacticHaloProperties for gNFW profile. More...
 
void GalacticHalo_Einasto (GalacticHaloProperties &result)
 Module function to generate GalacticHaloProperties for Einasto profile. More...
 
void ExtractLocalMaxwellianHalo (LocalMaxwellianHalo &result)
 Module function providing local density and velocity dispersion parameters. More...
 
void UnitTest_DarkBit (int &result)
 Central unit test routine. More...
 
void DarkMatter_ID_DiracSingletDM (std::string &result)
 
void DarkMatterConj_ID_DiracSingletDM (std::string &result)
 
void DD_couplings_DiracSingletDM_Z2 (DM_nucleon_couplings_fermionic_HP &result)
 Direct detection couplings for the DiracSingletDM_Z2 model. More...
 
void TH_ProcessCatalog_DiracSingletDM_Z2 (DarkBit::TH_ProcessCatalog &result)
 Set up process catalog for the DiracSingletDM_Z2 model. More...
 
void DD_couplings_DarkSUSY_DS5 (DM_nucleon_couplings &result)
 Get direct detection couplings from initialized DarkSUSY 5. More...
 
void DD_couplings_DarkSUSY_MSSM (DM_nucleon_couplings &result)
 Get direct detection couplings from DarkSUSY 6 initialized with MSSM module. More...
 
void DD_couplings_MicrOmegas (DM_nucleon_couplings &result)
 Get direct detection couplings from initialized MicrOmegas. More...
 
void sigma_SI_p_simple (double &result)
 Simple calculator of the spin-independent WIMP-proton cross-section. More...
 
void sigma_SI_n_simple (double &result)
 Simple calculator of the spin-independent WIMP-neutron cross-section. More...
 
void sigma_SD_p_simple (double &result)
 Simple calculator of the spin-dependent WIMP-proton cross-section. More...
 
void sigma_SD_n_simple (double &result)
 Simple calculator of the spin-dependent WIMP-neutron cross-section. More...
 
void sigma_SI_vnqn (map_intpair_dbl &result)
 Calculation of SI and SD cross sections at a reference momentum q0 for the fermionic Higgs portal models. More...
 
void sigma_SD_vnqn (map_intpair_dbl &result)
 
void GA_missingFinalStates (std::vector< std::string > &result)
 Identification of final states that are not yet tabulated. More...
 
daFunk::Funk boost_dNdE (daFunk::Funk dNdE, double gamma, double mass)
 Boosts an energy spectrum of isotropic particles into another frame (and isotropizes again). Parameters: gamma: Lorentz boost factor dNdE: Spectrum mass: mass of particle. More...
 
daFunk::Funk getYield (TH_Process process, double Ecm, double mass, TH_ProcessCatalog catalog, SimYieldTable table, double line_width, stringFunkMap cascadeMC_gammaSpectra)
 Helper function returning yield from a given DM process. More...
 
void GA_AnnYield_General (daFunk::Funk &result)
 General routine to derive annihilation yield. More...
 
void GA_DecayYield_General (daFunk::Funk &result)
 General routine to derive yield from decaying DM. More...
 
void SimYieldTable_DS5 (SimYieldTable &result)
 SimYieldTable based on DarkSUSY5 tabulated results. (DS6 below) More...
 
void SimYieldTable_DarkSUSY (SimYieldTable &result)
 SimYieldTable based on DarkSUSY6 tabulated results. More...
 
void SimYieldTable_MicrOmegas (SimYieldTable &result)
 SimYieldTable based on MicrOmegas tabulated results. More...
 
void SimYieldTable_PPPC (SimYieldTable &)
 SimYieldTable based on PPPC4DMID Cirelli et al. 2010. More...
 
void DarkMatter_ID_MajoranaSingletDM (std::string &result)
 
void DarkMatterConj_ID_MajoranaSingletDM (std::string &result)
 
void DD_couplings_MajoranaSingletDM_Z2 (DM_nucleon_couplings_fermionic_HP &result)
 Direct detection couplings for the MajoranaSingletDM_Z2 model. More...
 
void TH_ProcessCatalog_MajoranaSingletDM_Z2 (DarkBit::TH_ProcessCatalog &result)
 Set up process catalog for the MajoranaSingletDM_Z2 model. More...
 
double DSgamma3bdy (double(*IBfunc)(int &, double &, double &), int IBch, double Eg, double E1, double M_DM, double m_1, double m_2)
 Fully initialize DarkSUSY to the current model point. More...
 
void TH_ProcessCatalog_DS5_MSSM (DarkBit::TH_ProcessCatalog &result)
 Initialization of Process Catalog based on DarkSUSY 5 calculations. More...
 
void TH_ProcessCatalog_DS_MSSM (DarkBit::TH_ProcessCatalog &result)
 Initialization of Process Catalog based on DarkSUSY 6 calculations. More...
 
void DarkMatter_ID_MSSM (std::string &result)
 
void DarkMatterConj_ID_MSSM (std::string &result)
 
void RD_spectrum_MSSM (RD_spectrum_type &result)
 Collects spectrum information about coannihilating particles, resonances and threshold energies. More...
 
void RD_spectrum_SUSY_DS5 (RD_spectrum_type &result)
 Collects spectrum information about coannihilating particles, resonances and threshold energies – directly from DarkSUSY 5. More...
 
void RD_spectrum_from_ProcessCatalog (RD_spectrum_type &result)
 Collects information about resonances and threshold energies directly from the ProcessCatalog [NB: this assumes no coannihilating particles!]. More...
 
void RD_spectrum_ordered_func (RD_spectrum_type &result)
 Order RD_spectrum object and derive coannihilation thresholds. More...
 
void RD_annrate_DS5prep_func (int &result)
 Some helper function to prepare evaluation of Weff from DarkSUSY 5. More...
 
void RD_annrate_DSprep_MSSM_func (int &result)
 Some helper function to prepare evaluation of Weff from DarkSUSY 6. More...
 
void RD_eff_annrate_DS_MSSM (double(*&result)(double &))
 Get Weff directly from initialized DarkSUSY. Note that these functions do not (and should not) correct Weff for non-self-conjugate dark matter. More...
 
void RD_eff_annrate_DS5_MSSM (double(*&result)(double &))
 
 DEF_FUNKTRAIT (RD_EFF_ANNRATE_FROM_PROCESSCATALOG_TRAIT) void RD_eff_annrate_from_ProcessCatalog(double(*&result)(double &))
 Infer Weff from process catalog. More...
 
void RD_oh2_DS_general (double &result)
 General routine for calculation of relic density, using DarkSUSY 6+ Boltzmann solver. More...
 
void RD_oh2_DS5_general (double &result)
 General routine for calculation of relic density, using DarkSUSY 5 Boltzmann solver. More...
 
void RD_oh2_Xf_MicrOmegas (ddpair &result)
 Relic density directly from a call of initialized MicrOmegas. More...
 
void RD_oh2_DarkSUSY_DS5 (double &result)
 Relic density directly from a call of initialized DarkSUSY 5. More...
 
void RD_oh2_MicrOmegas (double &result)
 
void Xf_MicrOmegas (double &result)
 
void print_channel_contributions_MicrOmegas (double &result)
 
void get_semi_ann_MicrOmegas (double &result)
 
void vSigma_freezeout_MicrOmegas (double &result)
 Return the thermally averaged cross-section at T_freezeout. More...
 
void RD_fraction_one (double &result)
 
void RD_fraction_leq_one (double &result)
 
void RD_fraction_rescaled (double &result)
 
void DarkMatter_ID_ScalarSingletDM (std::string &result)
 
void DarkMatterConj_ID_ScalarSingletDM (std::string &result)
 
void get_ScalarSingletDM_DD_couplings (const Spectrum &spec, DM_nucleon_couplings &result, Models::safe_param_map< safe_ptr< const double > > &Param)
 Common code for different scalar singlet direct detection coupling routines. More...
 
void DD_couplings_ScalarSingletDM_Z2 (DM_nucleon_couplings &result)
 Direct detection couplings for Z2 scalar singlet DM. More...
 
void DD_couplings_ScalarSingletDM_Z3 (DM_nucleon_couplings &result)
 Direct detection couplings for Z3 scalar singlet DM. More...
 
void TH_ProcessCatalog_ScalarSingletDM_Z2 (DarkBit::TH_ProcessCatalog &result)
 Set up process catalog for Z2 scalar singlet DM. More...
 
void TH_ProcessCatalog_ScalarSingletDM_Z3 (DarkBit::TH_ProcessCatalog &result)
 Set up process catalog for Z3 scalar singlet DM. More...
 
void set_gamLike_GC_halo (bool &result)
 Fermi LAT dwarf likelihoods, based on arXiv:1108.2914. More...
 
void lnL_FermiLATdwarfs_gamLike (double &result)
 Fermi LAT dwarf likelihoods, using gamLike backend. More...
 
void lnL_HESSGC_gamLike (double &result)
 
void lnL_CTAGC_gamLike (double &result)
 
void lnL_FermiGC_gamLike (double &result)
 Fermi LAT galactic center likelihoods, using gamLike backend. More...
 
void lnL_oh2_Simple (double &result)
 Likelihood for cosmological relic density constraints. More...
 
void lnL_oh2_upperlimit (double &result)
 Likelihood for cosmological relic density constraints, implemented as an upper limit only Default data: Omega_c h^2 = 0.1188 +/- 0.0010 (1 sigma), Gaussian. More...
 
void lnL_sigmas_sigmal (double &result)
 Likelihoods for spin independent nuclear parameters. More...
 
void lnL_deltaq (double &result)
 Likelihoods for spin dependent nuclear parameters. More...
 
void lnL_rho0_lognormal (double &result)
 Likelihoods for halo parameters. More...
 
void lnL_vrot_gaussian (double &result)
 
void lnL_v0_gaussian (double &result)
 
void lnL_vesc_gaussian (double &result)
 
void dump_GammaSpectrum (double &result)
 Helper function to dump gamma-ray spectra. More...
 
void capture_rate_Sun_const_xsec_DS5 (double &result)
 Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1). DarkSUSY 5 version. More...
 
void capture_rate_Sun_const_xsec (double &result)
 Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1). DarkSUSY 6 version. More...
 
void capture_rate_Sun_const_xsec_capgen (double &result)
 Alternative to the darkSusy fct, using captn_specific from capgen instead. More...
 
void capture_rate_Sun_vnqn (double &result)
 Capture rate for v^n and q^n-dependent cross sections. More...
 
void equilibration_time_Sun (double &result)
 Equilibration time for capture and annihilation of dark matter in the Sun (s) More...
 
void annihilation_rate_Sun (double &result)
 Annihilation rate of dark matter in the Sun (s^-1) More...
 
void nuyield_from_DS (nuyield_info &result)
 Neutrino yield function pointer and setup. More...
 
void IC79_loglike (double &result)
 Composite IceCube 79-string likelihood function. More...
 
void IC_loglike (double &result)
 Complete composite IceCube likelihood function. More...
 
void DarkSUSY5_PointInit_LocalHalo_func (bool &result)
 Function to set Local Halo Parameters in DarkSUSY (DS5 only) More...
 
void DarkSUSY_PointInit_LocalHalo_func (bool &result)
 Function to set Local Halo Parameters in DarkSUSY (DS 6) More...
 
void DarkMatter_ID_VectorSingletDM (std::string &result)
 
void DarkMatterConj_ID_VectorSingletDM (std::string &result)
 
void DD_couplings_VectorSingletDM_Z2 (DM_nucleon_couplings &result)
 Direct detection couplings for the VectorSingletDM_Z2 model. More...
 
void TH_ProcessCatalog_VectorSingletDM_Z2 (DarkBit::TH_ProcessCatalog &result)
 Set up process catalog for the VectorSingletDM_Z2 model. More...
 
void IC22_full (nudata &result)
 Likelihood calculators for different IceCube event samples These functions all include the likelihood of the background-only model for the respective sample. More...
 
void IC79WH_full (nudata &result)
 79-string IceCube WH sample: predicted signal and background counts, observed counts and likelihoods. More...
 
void IC79WL_full (nudata &result)
 79-string IceCube WL sample: predicted signal and background counts, observed counts and likelihoods. More...
 
void IC79SL_full (nudata &result)
 79-string IceCube SL sample: predicted signal and background counts, observed counts and likelihoods. More...
 
void IC22_signal (double &result)
 22-string extractor module functions More...
 
void IC22_bg (double &result)
 
void IC22_nobs (int &result)
 
void IC22_loglike (double &result)
 
void IC22_bgloglike (double &result)
 
void IC22_pvalue (double &result)
 
void IC79WH_signal (double &result)
 79-string WH extractor module functions More...
 
void IC79WH_bg (double &result)
 
void IC79WH_nobs (int &result)
 
void IC79WH_loglike (double &result)
 
void IC79WH_bgloglike (double &result)
 
void IC79WH_pvalue (double &result)
 
void IC79WL_signal (double &result)
 79-string WL extractor module functions More...
 
void IC79WL_bg (double &result)
 
void IC79WL_nobs (int &result)
 
void IC79WL_loglike (double &result)
 
void IC79WL_bgloglike (double &result)
 
void IC79WL_pvalue (double &result)
 
void IC79SL_signal (double &result)
 79-string SL extractor module functions More...
 
void IC79SL_bg (double &result)
 
void IC79SL_nobs (int &result)
 
void IC79SL_loglike (double &result)
 
void IC79SL_bgloglike (double &result)
 
void IC79SL_pvalue (double &result)
 

Variables

const double gagg_conversion = 1.0E-9
 Supporting classes and functions for the axion module. More...
 
const double gaee_conversion = 1.0E+13
 
const std::map< InterpolationOptions1D, std::string > int_type_name = { { InterpolationOptions1D::linear, "linear" }, { InterpolationOptions1D::cspline, "cspline"} }
 
const std::map< InterpolationOptions2D, std::string > int_2d_type_name = { { InterpolationOptions2D::bilinear, "bilinear" }, { InterpolationOptions2D::bicubic, "bicubic"} }
 
const double abs_prec = 1.0E-1
 
const double rel_prec = 1.0E-6
 
const int method = 5
 

Typedef Documentation

◆ simpleHistContainter

Definition at line 99 of file SimpleHist.hpp.

◆ stringFunkMap

Definition at line 101 of file SimpleHist.hpp.

◆ stringIntMap

Definition at line 100 of file SimpleHist.hpp.

Enumeration Type Documentation

◆ cascadeMC_SpecialEvents

Special events for event loop.

Enumerator
MC_INIT 
MC_NEXT_STATE 
MC_FINALIZE 

Definition at line 38 of file Cascades.cpp.

◆ InterpolationOptions1D

Generic one-dimensional integration container for linear interpolation and cubic splines.

Enumerator
linear 
cspline 

Definition at line 80 of file Axions.cpp.

◆ InterpolationOptions2D

Two-dimensional integration container for bilinear interpolation and bicubic splines.

Enumerator
bilinear 
bicubic 

Definition at line 197 of file Axions.cpp.

Function Documentation

◆ ALPS1_lnL_general()

double Gambit::DarkBit::ALPS1_lnL_general ( double  s,
double  mu,
double  sigma 
)

Definition at line 1238 of file Axions.cpp.

Referenced by calc_lnL_ALPS1().

1239  {
1240  // Propagate uncertainty from efficiency in chi^2-likelihood.
1241  const double rel_error = 0.05/0.82;
1242  return -0.5*gsl_pow_2(s-mu)/(gsl_pow_2(rel_error*s)+gsl_pow_2(sigma));
1243  }
const double sigma
Definition: SM_Z.hpp:43
const double mu
Definition: SM_Z.hpp:42
Here is the caller graph for this function:

◆ ALPS1_signal_general()

double Gambit::DarkBit::ALPS1_signal_general ( double  power,
double  nm1,
double  m_ax,
double  gagg 
)

Likelihoods for ALPS 1 (LSW), CAST (helioscopes), and ADMX, UF, RBF (haloscopes).

Supported models: GeneralALP

Definition at line 1191 of file Axions.cpp.

References Gambit::gev2cm, and Gambit::pi.

Referenced by calc_ALPS1_signal_gas(), and calc_ALPS1_signal_vac().

1192  {
1193  const double eVm = gev2cm*1E7;
1194  // Photon energy in eV.
1195  const double erg = 2.0*pi*eVm/532.0E-9;
1196  const double len = 4.2;
1197  // We include the uncertainty of the detection efficiency eff = 0.82(5) in the likelihood.
1198  const double eff = 0.82;
1199 
1200  double result = 0.0;
1201 
1202  // CAVE: Approximations/conversion only valid/possible for m_ax << 2.33 eV (532 nm).
1203  if (m_ax < 1.0)
1204  {
1205  // Effective photon mass and momentum transfer.
1206  double m_ph_sq = 2.0*erg*erg*nm1;
1207  double q = 0.5*(m_ax*m_ax+m_ph_sq)/(eVm*erg);
1208  double factor = gsl_pow_4((gagg*1E17)*gsl_sf_sinc(0.5*q*len/pi));
1209 
1210  // Prefactor: 1096 W * 1 h * (10^-17/eV * 4.98 T * 4.2 m)^4 / 16.
1211  result = 0.00282962979*eff*factor*(power/1096.0)/erg;
1212  }
1213 
1214  return result;
1215  }
const double pi
const double gev2cm
Here is the caller graph for this function:

◆ alt_erg_integrand()

double Gambit::DarkBit::alt_erg_integrand ( double  erg,
void params 
)

Definition at line 741 of file Axions.cpp.

References abs_prec, DarkBit_error(), Gambit::EOM, erg_integrand(), Gambit::Utils::file_exists(), Gambit::FlavBit::G(), Gambit::gev2cm, Gambit::LogTags::info, daFunk::interp(), LOCAL_INFO, Gambit::logger(), Gambit::pi, Gambit::Scanner::pow(), rel_prec, and rs.

742  {
743  const double eVm = gev2cm*1E7;
744  const double L = 9.26/eVm;
745  struct SolarModel_params4 * p4 = (struct SolarModel_params4 *)params;
746  double m_ax = p4->ma0;
747 
748  double argument = 0.25*1.0E-3*L*m_ax*m_ax/erg;
749  double temp = gsl_pow_2(gsl_sf_sinc(argument/pi));
750  double exposure = p4->eff_exp->interpolate(erg);
751  double gaee_flux = p4->gaee_flux->interpolate(erg);
752 
753  return temp*exposure*gaee_flux;
754  }
const double pi
const double gev2cm
Here is the call graph for this function:

◆ annihilation_rate_Sun()

void Gambit::DarkBit::annihilation_rate_Sun ( double result)

Annihilation rate of dark matter in the Sun (s^-1)

Definition at line 215 of file SunNeutrinos.cpp.

References equilibration_time_Sun(), and Gambit::Scanner::pow().

Referenced by IC22_full(), IC79SL_full(), IC79WH_full(), IC79WL_full(), and main().

216  {
217  using namespace Pipes::annihilation_rate_Sun;
218  double tt_sun = 1.5e17 / *Dep::equilibration_time_Sun;
219  result = *Dep::capture_rate_Sun * 0.5 * pow(tanh(tt_sun),2.0);
220  }
void equilibration_time_Sun(double &result)
Equilibration time for capture and annihilation of dark matter in the Sun (s)
void annihilation_rate_Sun(double &result)
Annihilation rate of dark matter in the Sun (s^-1)
double pow(const double &a)
Outputs a^i.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ axion_mass_temp()

double Gambit::DarkBit::axion_mass_temp ( double  T,
double  beta,
double  Tchi 
)

Definition at line 1648 of file Axions.cpp.

References beta, Gambit::Scanner::pow(), Tchi, and thetai.

Referenced by equation_Tosc(), RD_oh2_Axions(), scal_field_eq(), and scal_field_eq_jac().

1649  {
1650  double res = 1.0;
1651  if (T > Tchi) { res = pow(T/Tchi,-0.5*beta); }
1652  return res;
1653  }
START_MODEL beta
Definition: Axions.hpp:36
START_MODEL Tchi
Definition: Axions.hpp:36
double pow(const double &a)
Outputs a^i.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ boost_dNdE()

daFunk::Funk Gambit::DarkBit::boost_dNdE ( daFunk::Funk  dNdE,
double  gamma,
double  mass 
)

Boosts an energy spectrum of isotropic particles into another frame (and isotropizes again). Parameters: gamma: Lorentz boost factor dNdE: Spectrum mass: mass of particle.

Definition at line 154 of file GamYields.cpp.

References DarkBit_error(), Gambit::FlavBit::LoopFunctions::E(), LOCAL_INFO, and daFunk::var().

Referenced by getYield().

155  {
156  if ( gamma < 1.0 + .02 ) // Ignore less than 2% boosts
157  {
158  if (gamma < 1.0)
159  DarkBit_error().raise(LOCAL_INFO,
160  "boost_dNdE: Requested Lorentz boost with gamma < 1");
161  return dNdE;
162  }
163  double betaGamma = sqrt(gamma*gamma-1);
164  daFunk::Funk E = daFunk::var("E");
165  daFunk::Funk lnE = daFunk::var("lnE");
166  daFunk::Funk Ep = daFunk::var("Ep");
167  daFunk::Funk halfBox_int = betaGamma*sqrt(E*E-mass*mass);
168  daFunk::Funk halfBox_bound = betaGamma*sqrt(Ep*Ep-mass*mass);
169  daFunk::Funk integrand = dNdE/(2*halfBox_int);
170  return integrand->gsl_integration("E", Ep*gamma-halfBox_bound, Ep*gamma+halfBox_bound)
171  ->set_epsabs(0)->set_limit(100)->set_epsrel(1e-3)->set_use_log_fallback(true)->set("Ep", daFunk::var("E"));
172  //
173  // Note: integration over lnE causes problems in the WIMP example (3) as the singularity is dropped.
174  // return (integrand*E)->set("E", exp(lnE))->gsl_integration("lnE", log(Ep*gamma-halfBox_bound), log(Ep*gamma+halfBox_bound))
175  // ->set_epsabs(0)->set_epsrel(1e-3)->set("Ep", daFunk::var("E"));
176  }
error & DarkBit_error()
#define LOCAL_INFO
Definition: local_info.hpp:34
Funk var(std::string arg)
Definition: daFunk.hpp:921
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
double E(const double x, const double y)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ calc_ALPS1_signal_gas()

void Gambit::DarkBit::calc_ALPS1_signal_gas ( double result)

Definition at line 1228 of file Axions.cpp.

References ALPS1_signal_general(), and Gambit::FlavBit::LoopFunctions::E().

1229  {
1230  using namespace Pipes::calc_ALPS1_signal_gas;
1231  double m_ax = *Param["ma0"];
1232  double gagg = gagg_conversion*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1233 
1234  result = ALPS1_signal_general(1044.0, 5.0E-8, m_ax, gagg);
1235  }
const double gagg_conversion
Supporting classes and functions for the axion module.
Definition: Axions.cpp:70
void calc_ALPS1_signal_gas(double &result)
Definition: Axions.cpp:1228
double ALPS1_signal_general(double power, double nm1, double m_ax, double gagg)
Likelihoods for ALPS 1 (LSW), CAST (helioscopes), and ADMX, UF, RBF (haloscopes). ...
Definition: Axions.cpp:1191
double E(const double x, const double y)
Here is the call graph for this function:

◆ calc_ALPS1_signal_vac()

void Gambit::DarkBit::calc_ALPS1_signal_vac ( double result)

Definition at line 1218 of file Axions.cpp.

References ALPS1_signal_general().

1219  {
1220  using namespace Pipes::calc_ALPS1_signal_vac;
1221  double m_ax = *Param["ma0"];
1222  double gagg = gagg_conversion*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1223 
1224  result = ALPS1_signal_general(1096.0, 0.0, m_ax, gagg);
1225  }
const double gagg_conversion
Supporting classes and functions for the axion module.
Definition: Axions.cpp:70
void calc_ALPS1_signal_vac(double &result)
Definition: Axions.cpp:1218
double ALPS1_signal_general(double power, double nm1, double m_ax, double gagg)
Likelihoods for ALPS 1 (LSW), CAST (helioscopes), and ADMX, UF, RBF (haloscopes). ...
Definition: Axions.cpp:1191
Here is the call graph for this function:

◆ calc_AxionOscillationTemperature()

void Gambit::DarkBit::calc_AxionOscillationTemperature ( double result)

Definition at line 1676 of file Axions.cpp.

References beta, Gambit::FlavBit::LoopFunctions::E(), equation_Tosc(), Gambit::m_planck_red, Gambit::pi, Gambit::Scanner::pow(), r, and Tchi.

1677  {
1679 
1680  double ma0 = *Param["ma0"];
1681  double beta = *Param["beta"];
1682  double Tchi = *Param["Tchi"];
1683  // m_pl/10^12 eV = m_pl/10^3 GeV
1684  const double m_pl = m_planck_red*1.0E-3;
1685 
1686  // Initialising the parameter structure with known and dummy values.
1687  AxionEDT_params params = {ma0, beta, Tchi, 0.0, 0.0};
1688 
1689  // Use gsl implementation Brent's method to determine the oscillation temperature.
1690  gsl_function F;
1691  F.function = &equation_Tosc;
1692  F.params = &params;
1693 
1694  // Set counters and initialise equation solver.
1695  int status;
1696  int iter = 0, max_iter = 1E6;
1697  gsl_root_fsolver *s;
1698  s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
1699  double r, r_up, r_lo;
1700 
1701  // Calculate first estimate for the root bracketing [r_lo, r_up].
1702  // Calculate best estimate for comparison. g(Tchi)^-0.25 = 0.49, g(Tchi)^-...=0.76
1703  r = 0.49*pow((10.0/(pi*pi)) * gsl_pow_2(m_pl*ma0), 0.25);
1704  // Compare to decide which branch of the equation is valid; T1 > Tchi or T1 < Tchi
1705  if ( (r > Tchi) && (beta > 1.0E-10) )
1706  {
1707  r = 0.76*pow((10.0/(pi*pi)) * gsl_pow_2(m_pl*ma0) * pow(Tchi, beta), 1.0/(4.0+beta));
1708  }
1709  // Find appropriate values for r_lo and r_up
1710  r_up = r;
1711  r_lo = r;
1712  while (GSL_FN_EVAL(&F,r_up) > 0.0) { r_up = 2.0*r_up; }
1713  while (GSL_FN_EVAL(&F,r_lo) < 0.0) { r_lo = 0.5*r_lo; }
1714 
1715  // Execute equation solver until we reach 10^-6 absolute precision.
1716  gsl_root_fsolver_set(s, &F, r_lo, r_up);
1717  do
1718  {
1719  iter++;
1720  status = gsl_root_fsolver_iterate (s);
1721  r = gsl_root_fsolver_root (s);
1722  r_lo = gsl_root_fsolver_x_lower (s);
1723  r_up = gsl_root_fsolver_x_upper (s);
1724  status = gsl_root_test_interval (r_lo, r_up, 1.0E-6, 0.0);
1725  } while (status == GSL_CONTINUE && iter < max_iter);
1726 
1727  gsl_root_fsolver_free (s);
1728 
1729  result = r;
1730  }
START_MODEL beta
Definition: Axions.hpp:36
const double m_planck_red
const double pi
START_MODEL dNur_CMB r
START_MODEL Tchi
Definition: Axions.hpp:36
double pow(const double &a)
Outputs a^i.
double equation_Tosc(double T, void *params)
Definition: Axions.cpp:1661
double E(const double x, const double y)
void calc_AxionOscillationTemperature(double &result)
Definition: Axions.cpp:1676
Here is the call graph for this function:

◆ calc_CAST2007_signal_vac()

void Gambit::DarkBit::calc_CAST2007_signal_vac ( std::vector< double > &  result)

Definition at line 1274 of file Axions.cpp.

References Gambit::Scanner::pow().

1275  {
1276  using namespace Pipes::calc_CAST2007_signal_vac;
1277  double m_ax = *Param["ma0"];
1278  double gagg = gagg_conversion*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1279  double gaee = std::fabs(*Param["gaee"]);
1280 
1281  // Initialise the Solar model calculator and get the reference counts for a given mass.
1282  // Get Solar model we are working with; set default value here
1283  static std::string solar_model_gagg = runOptions->getValueOrDef<std::string> ("AGSS09met", "solar_model_gagg");
1284  static std::string solar_model_gaee = runOptions->getValueOrDef<std::string> ("AGSS09met_old", "solar_model_gaee");
1285  static CAST_SolarModel_Interpolator lg_ref_counts (solar_model_gagg, solar_model_gaee, "CAST2007");
1286  std::vector<double> lg_ref_counts_gagg = lg_ref_counts.evaluate_gagg_contrib(m_ax);
1287  std::vector<double> lg_ref_counts_gaee = lg_ref_counts.evaluate_gaee_contrib(m_ax);
1288  static int n_bins = lg_ref_counts_gagg.size();
1289 
1290  std::vector<double> counts;
1291  double dummy;
1292  for (int i = 0; i < n_bins; i++)
1293  {
1294  dummy = gsl_pow_2(gagg*1E19)*pow(10,lg_ref_counts_gagg[i]) + gsl_pow_2(gaee*gaee_conversion)*pow(10,lg_ref_counts_gaee[i]);
1295  counts.push_back(gsl_pow_2(gagg*1E19)*dummy);
1296  }
1297 
1298  result = counts;
1299  }
const double gagg_conversion
Supporting classes and functions for the axion module.
Definition: Axions.cpp:70
void calc_CAST2007_signal_vac(std::vector< double > &result)
Definition: Axions.cpp:1274
const double gaee_conversion
Definition: Axions.cpp:71
double pow(const double &a)
Outputs a^i.
Here is the call graph for this function:

◆ calc_CAST2017_signal_vac()

void Gambit::DarkBit::calc_CAST2017_signal_vac ( std::vector< std::vector< double >> &  result)

Definition at line 1302 of file Axions.cpp.

References Gambit::Scanner::pow().

1303  {
1304  using namespace Pipes::calc_CAST2017_signal_vac;
1305  double m_ax = *Param["ma0"];
1306  double gagg = gagg_conversion*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1307  double gaee = std::fabs(*Param["gaee"]);
1308  std::vector<std::vector<double>> res;
1309 
1310  // Initialise the Solar model calculator and get the reference counts for a given mass.
1311  // Get Solar model we are working with; set default value here
1312  static std::string solar_model_gagg = runOptions->getValueOrDef<std::string> ("AGSS09met", "solar_model_gagg");
1313  static std::string solar_model_gaee = runOptions->getValueOrDef<std::string> ("AGSS09met_old", "solar_model_gaee");
1314 
1315  const int n_exps = 12;
1316  const int n_bins = 10;
1317  const std::string exp_names [n_exps] = {"A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"};
1318  static std::vector<CAST_SolarModel_Interpolator> lg_ref_counts;
1319  static bool lg_ref_counts_not_calculated = true;
1320  if (lg_ref_counts_not_calculated)
1321  {
1322  for (int e = 0; e < n_exps; e++)
1323  {
1324  CAST_SolarModel_Interpolator dummy (solar_model_gagg, solar_model_gaee, "CAST2017_"+exp_names[e]);
1325  lg_ref_counts.push_back(std::move(dummy));
1326  }
1327  }
1328  lg_ref_counts_not_calculated = false;
1329 
1330  for (int e = 0; e < n_exps; e++)
1331  {
1332  std::vector<double> lg_ref_counts_gagg = lg_ref_counts[e].evaluate_gagg_contrib(m_ax);
1333  std::vector<double> lg_ref_counts_gaee = lg_ref_counts[e].evaluate_gaee_contrib(m_ax);
1334 
1335  std::vector<double> counts;
1336  double dummy;
1337  for (int bin = 0; bin < n_bins; bin++)
1338  {
1339  dummy = gsl_pow_2(gagg*1E19)*pow(10,lg_ref_counts_gagg[bin]) + gsl_pow_2(gaee*gaee_conversion)*pow(10,lg_ref_counts_gaee[bin]);
1340  counts.push_back(gsl_pow_2(gagg*1E19)*dummy);
1341  }
1342 
1343  res.push_back(counts);
1344  }
1345 
1346  result = res;
1347  }
const double gagg_conversion
Supporting classes and functions for the axion module.
Definition: Axions.cpp:70
void calc_CAST2017_signal_vac(std::vector< std::vector< double >> &result)
Definition: Axions.cpp:1302
const double gaee_conversion
Definition: Axions.cpp:71
double pow(const double &a)
Outputs a^i.
Here is the call graph for this function:

◆ calc_Haloscope_signal()

void Gambit::DarkBit::calc_Haloscope_signal ( double result)

Definition at line 1447 of file Axions.cpp.

References Gambit::FlavBit::LoopFunctions::E(), and Gambit::LocalMaxwellianHalo::rho0.

1448  {
1449  using namespace Pipes::calc_Haloscope_signal;
1450  double gagg = gagg_conversion*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1451  // Get the DM fraction in axions and the local DM density.
1452  double fraction = *Dep::RD_fraction;
1453  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
1454  double rho0 = LocalHaloParameters.rho0;
1455 
1456  // Signal relative to a reference coupling and local DM density.
1457  double s = gsl_pow_2(gagg/1.0E-24) * fraction * (rho0/0.45);
1458 
1459  result = s;
1460  }
const double gagg_conversion
Supporting classes and functions for the axion module.
Definition: Axions.cpp:70
void calc_Haloscope_signal(double &result)
Definition: Axions.cpp:1447
double E(const double x, const double y)
Here is the call graph for this function:

◆ calc_lnL_ALPS1()

void Gambit::DarkBit::calc_lnL_ALPS1 ( double result)

Definition at line 1246 of file Axions.cpp.

References ALPS1_lnL_general().

1247  {
1248  using namespace Pipes::calc_lnL_ALPS1;
1249  double s1 = *Dep::ALPS1_signal_vac;
1250  double s2 = *Dep::ALPS1_signal_gas;
1251 
1252  // ALPS Collaboration results (limits from this data published in 1004.1313).
1253  // ALPS Collaboration results, vacuum, 5 frames.
1254  const double exp_sig_mu_v1 = -4.01, exp_sig_std_v1 = 3.01;
1255  // ALPS Collaboration results, vacuum, 6 frames.
1256  const double exp_sig_mu_v2 = -2.35, exp_sig_std_v2 = 3.44;
1257  // ALPS Collaboration results, vacuum combined(!), 11 frames (we keep them seperated).
1258  //const double exp_sig_mu_vc = -3.29, exp_sig_std_vc = 2.27;
1259  // ALPS Collaboration results, gas, 8 frames (P = 0.18 mbar).
1260  const double exp_sig_mu_g = 3.98, exp_sig_std_g = 2.45;
1261 
1262  double l1 = ALPS1_lnL_general(s1, exp_sig_mu_v1, exp_sig_std_v1);
1263  double l2 = ALPS1_lnL_general(s1, exp_sig_mu_v2, exp_sig_std_v2);
1264  double l3 = ALPS1_lnL_general(s2, exp_sig_mu_g, exp_sig_std_g);
1265 
1266  result = l1 + l2 + l3;
1267  }
double ALPS1_lnL_general(double s, double mu, double sigma)
Definition: Axions.cpp:1238
void calc_lnL_ALPS1(double &result)
Definition: Axions.cpp:1246
Here is the call graph for this function:

◆ calc_lnL_CAST2007()

void Gambit::DarkBit::calc_lnL_CAST2007 ( double result)

Definition at line 1365 of file Axions.cpp.

References CAST_lnL_general().

1366  {
1367  using namespace Pipes::calc_lnL_CAST2007;
1368  std::vector<double> sig_vac = *Dep::CAST2007_signal_vac;
1369 
1370  // CAST CCD 2004 vacuum data (based on hep-ex/0702006).
1371  const int n_bins = 20;
1372  const std::vector<int> dat_vac {1, 3, 1, 1, 1, 2, 1, 2, 0, 2, 0, 1, 0, 2, 2, 0, 2, 1, 2, 2};
1373  const std::vector<double> bkg_vac {2.286801272, 1.559182673, 2.390746817, 1.559182673, 2.598637835, 1.039455092, 0.727618599, 1.559182673, 1.247346181, 1.455237199, 1.871019235, 0.831564073, 1.663128217, 1.247346181, 1.143400636, 1.663128217,
1374  1.247346181, 1.247346181, 2.286801272, 1.247346181};
1375 
1376  // Only calculate norm once.
1377  static double norm = 0.0;
1378  static bool norm_not_calculated = true;
1379  if (norm_not_calculated)
1380  {
1381  for (int i = 0; i < n_bins; i++) { norm += gsl_sf_lnfact(dat_vac[i]); }
1382  }
1383  norm_not_calculated = false;
1384 
1385  result = CAST_lnL_general(sig_vac, bkg_vac, dat_vac) - norm;
1386  }
void calc_lnL_CAST2007(double &result)
Definition: Axions.cpp:1365
double CAST_lnL_general(std::vector< double > s, const std::vector< double > bkg_counts, const std::vector< int > sig_counts)
Definition: Axions.cpp:1350
Here is the call graph for this function:

◆ calc_lnL_CAST2017()

void Gambit::DarkBit::calc_lnL_CAST2017 ( double result)

Definition at line 1389 of file Axions.cpp.

References CAST_lnL_general().

1390  {
1391  using namespace Pipes::calc_lnL_CAST2017;
1392  std::vector<std::vector<double>> sig_vac = *Dep::CAST2017_signal_vac;
1393 
1394  // CAST 2017 vacuum data (naming scheme based on the 10 data sets published in 1705.02290).
1395  const int n_bins = 10;
1396  const int n_exps = 12;
1397 
1398  const std::vector<std::vector<int>> dat_vac_all { {0, 3, 3, 0, 0, 1, 3, 3, 3, 3},
1399  {5, 5, 5, 3, 3, 0, 5, 2, 2, 1},
1400  {3, 3, 1, 2, 2, 2, 4, 5, 4, 3},
1401  {1, 5, 5, 2, 1, 2, 2, 5, 4, 0},
1402  {2, 3, 2, 2, 2, 1, 0, 2, 1, 1},
1403  {3, 5, 1, 4, 1, 2, 0, 3, 2, 2},
1404  {3, 4, 4, 5, 1, 2, 3, 2, 3, 2},
1405  {2, 1, 0, 1, 3, 2, 2, 3, 0, 1},
1406  {1, 2, 2, 1, 3, 0, 0, 1, 4, 0},
1407  {2, 1, 3, 1, 1, 0, 1, 1, 5, 5},
1408  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
1409  {0, 2, 1, 0, 0, 0, 0, 0, 0, 0} };
1410 
1411  const std::vector<std::vector<double>> bkg_vac_all { {0.926256, 1.96148, 1.79803, 1.30766, 1.30766, 1.96148, 2.61531, 2.77877, 2.94223, 2.07045},
1412  {3.68151, 4.86486, 4.99634, 3.55003, 2.49817, 3.28707, 2.89262, 3.68151, 3.48429, 3.41855},
1413  {2.54573, 3.18216, 4.45502, 2.86394, 2.29116, 2.29116, 3.30945, 3.75495, 3.62766, 3.56402},
1414  {2.72482, 5.5794, 3.95748, 2.40044, 2.27069, 2.33556, 3.37359, 3.43847, 3.24384, 3.11408},
1415  {1.44613, 2.30066, 2.43213, 1.70906, 1.97199, 1.24893, 1.24893, 2.23493, 2.16919, 2.23493},
1416  {1.30963, 2.94666, 2.35733, 2.55377, 2.02992, 1.50607, 2.16088, 2.75022, 2.29185, 2.29185},
1417  {2.33334, 2.74167, 2.21667, 2.80001, 2.21667, 1.75001, 2.62501, 2.21667, 2.80001, 2.33334},
1418  {1.74724, 2.37125, 2.68326, 1.62243, 2.05924, 1.74724, 1.49763, 1.74724, 1.18563, 2.24645},
1419  {1.72998, 3.45995, 1.79405, 1.72998, 1.9222, 1.72998, 2.69107, 2.24256, 1.98627, 2.11442},
1420  {1.89627, 2.25182, 2.96292, 1.4222, 1.65924, 1.65924, 1.95553, 2.1333, 1.71849, 2.07404},
1421  {0.0150685, 0.0568493, 0.060274, 0.0150685, 0.0150685, 0.00753425, 0.0267123, 0.0150685, 0.0267123, 0.0116438},
1422  {0.0409574, 0.226904, 0.243287, 0.0532447, 0.0188404, 0.0344043, 0.0417766, 0.0409574, 0.0409574, 0.0286702} };
1423 
1424  // Only calculate norm once.
1425  static double norm = 0.0;
1426  static bool norm_not_calculated = true;
1427  if (norm_not_calculated)
1428  {
1429  for (int bin = 0; bin < n_bins; bin++)
1430  {
1431  for (int e = 0; e < n_exps; e++) { norm += gsl_sf_lnfact(dat_vac_all[e][bin]); }
1432  }
1433  }
1434  norm_not_calculated = false;
1435 
1436  result = 0.0;
1437  for (int e = 0; e < n_exps; e++) { result = result + CAST_lnL_general(sig_vac[e], bkg_vac_all[e], dat_vac_all[e]); }
1438  result = result - norm;
1439  }
void calc_lnL_CAST2017(double &result)
Definition: Axions.cpp:1389
double CAST_lnL_general(std::vector< double > s, const std::vector< double > bkg_counts, const std::vector< int > sig_counts)
Definition: Axions.cpp:1350
Here is the call graph for this function:

◆ calc_lnL_Haloscope_ADMX1()

void Gambit::DarkBit::calc_lnL_Haloscope_ADMX1 ( double result)

Approximated likelihood for the AxionDarkMatterEXperiment (ADMX).

Definition at line 1466 of file Axions.cpp.

References b, f, Gambit::hbar, combine_hdf5::index, and Gambit::pi.

1467  {
1468  using namespace Pipes::calc_lnL_Haloscope_ADMX1;
1469  double m_ax = *Param["ma0"];
1470  // Calculate equivalent frequency in MHz.
1471  const double eV_to_MHz = 1.0E-15/(2.0*pi*hbar);
1472  double freq = m_ax*eV_to_MHz;
1473  double l = 0.0;
1474  // Initialise GSL histogram and flag.
1475  static gsl_histogram *h = gsl_histogram_alloc(89);
1476  static bool init_flag = false;
1477 
1478  // Unless initialised already, read in digitised limits from 0910.5914.
1479  if (not(init_flag))
1480  {
1481  FILE * f = fopen(GAMBIT_DIR "/DarkBit/data/ADMXLimitsHistogram.dat", "r");
1482  gsl_histogram_fscanf (f, h);
1483  fclose(f);
1484  init_flag = true;
1485  }
1486 
1487  // Likelihood shape parameters based on limits from astro-ph/9801286.
1488  const double a = 0.013060890;
1489  const double b = 0.455482976;
1490 
1491  if ((freq > gsl_histogram_min(h)) && (freq < gsl_histogram_max(h)))
1492  {
1493  size_t index;
1494  gsl_histogram_find(h, freq, &index);
1495  double s = *Dep::Haloscope_signal;
1496  double s_ref = gsl_pow_2(gsl_histogram_get(h, index));
1497  double s_rel = s/s_ref;
1498  // Only apply contraints for a signal > threshold a.
1499  if (s_rel > a) { l = -0.5 * gsl_pow_2( (s_rel - a)/b ); }
1500  }
1501 
1502  result = l;
1503  }
START_MODEL b
Definition: demo.hpp:270
void calc_lnL_Haloscope_ADMX1(double &result)
Definition: Axions.cpp:1466
const double hbar
const double pi

◆ calc_lnL_Haloscope_ADMX2()

void Gambit::DarkBit::calc_lnL_Haloscope_ADMX2 ( double result)

Definition at line 1506 of file Axions.cpp.

1507  {
1508  using namespace Pipes::calc_lnL_Haloscope_ADMX2;
1509  // Rescale the axion mass to ueV.
1510  double m_ax = 1.0E+6*(*Param["ma0"]);
1511  double l = 0.0;
1512 
1513  // ADMX 2018 90% C.L. exclusion limits; digitised from Fig. 4, 1804.05750.
1514  static AxionInterpolator g_limits (GAMBIT_DIR "/DarkBit/data/ADMX2018Limits.dat");
1515 
1516  // If we are within the avialable data range, calculate the limit.
1517  if ( (m_ax > g_limits.lower()) && (m_ax < g_limits.upper()) )
1518  {
1519  double s = *Dep::Haloscope_signal;
1520  // Get limit and rescale it to 1 sigma from the appropriate number of sigmas for 90% C.L. (1 d.o.f.).
1521  double sigma_exp = gsl_pow_2(g_limits.interpolate(m_ax))/1.644817912489;
1522  // Add systematics of 13% according to 1804.05750.
1523  double var_exp = gsl_pow_2(sigma_exp);
1524  double var_theo = gsl_pow_2(0.13*s);
1525  double var_tot = var_exp + var_theo;
1526  l = -0.5*gsl_pow_2(s)/var_tot;
1527  }
1528 
1529  result = l;
1530  }
void calc_lnL_Haloscope_ADMX2(double &result)
Definition: Axions.cpp:1506

◆ calc_lnL_Haloscope_RBF()

void Gambit::DarkBit::calc_lnL_Haloscope_RBF ( double result)

Definition at line 1557 of file Axions.cpp.

References combine_hdf5::index, and Gambit::DecayBit::SM_Z::sigma.

1558  {
1559  using namespace Pipes::calc_lnL_Haloscope_RBF;
1560  // Rescale the axion mass to ueV.
1561  double m = 1.0E+6*(*Param["ma0"]);
1562  double l = 0.0;
1563  // Results from Phys. Rev. D 40, 3153 (1989)
1564  // The bins and results below (from Table I) appear to be inconsitent with Figure 14.
1565  //const std::vector<double> bins = {4.507875, 5.037241, 5.459079, 5.851967, 5.996715, 6.257262, 6.617065, 6.976868, 7.113345, 7.295314, 7.857765, 8.631134, 8.684898, 9.259755, 9.760171, 10.173737,
1566  // 11.298638, 11.583999, 12.845377, 13.234130, 15.301962, 16.2655809};
1567  // N_sigma values as defined in the paper.
1568  //const double N_sigma [21] = {5.0, 5.0, 5.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 0.0, 5.0, 4.0, 4.0, 4.0, 4.0};
1569  // Proportionality factors ("sigma") inferred from Table I in in units of GeV^-1/cm^3.
1570  //const double eff_sigma [21] = {8.3030524510E+01, 3.5789241789E+02, 5.3457189090E+02, 8.3673921774E+02, 7.3205267295E+02, 7.1850356207E+02, 7.0099211538E+02, 9.3243407987E+02, 1.3132694610E+03, 1.9447760075E+03, 2.4028734743E+03,
1571  // 3.5992849457E+03, 5.8433323192E+03, 1.2415907565E+03, 1.1487509033E+03, 1.0000000000E+99, 2.6768234439E+03, 9.1546564260E+04, 1.7208310692E+04, 4.2462784870E+04, 2.8794160844E+04};
1572  // The results below are derived from Fig. 14 (assuming N_sigma = 4 for all values).
1573  const std::vector<double> bins = {4.400841, 4.960600, 5.209095, 5.668611, 6.934590, 7.445686, 8.041207, 8.898392, 9.570607, 10.067396, 11.213613, 11.626834, 12.773085, 13.450179, 14.704884, 16.170394};
1574  const double N_sigma = 4.0;
1575  const double eff_sigma [15] = {7.794388E+01, 3.808827E+02, 5.328136E+02, 6.765588E+02, 1.772892E+03, 2.752458E+03, 5.945156E+03, 2.025315E+03, 1.546855E+03, 1.022957E+13, 5.464075E+03, 9.621171E+04, 2.023187E+04, 5.201449E+04, 3.597168E+04};
1576 
1577  // Likelihood shape parameters based on information from PRD 40 (1989) 3153.
1578  // Note that there are no limits between 10.1 and 11.2 ueV; the tabulated value in that bin is just a placeholder, which is never used.
1579  if ( ((m > bins.front()) && (m < 10.067396)) || ((m > 11.213613) && (m < bins.back())))
1580  {
1581  // Use the standard search algorthim to identify the bin index and use the appropriate values for the likelihood.
1582  auto index = upper_bound(bins.begin(), bins.end(), m);
1583  double sigma = eff_sigma[index-bins.begin()-1];
1584  // Uncomment and comment out lines below to swap between implementations using Table I and Figure 14, respectively.
1585  //double offset = sigma*N_sigma[index-bins.begin()-1];
1586  double offset = sigma*N_sigma;
1587  // The "signal" needs to be rescaled to 0.3 GeV/cm^3, which is their reference value.
1588  double s = (0.45/0.3)*(*Dep::Haloscope_signal);
1589  if (s > offset) {
1590  l = -0.5 * gsl_pow_2( (s - offset)/sigma );
1591  }
1592  }
1593 
1594  result = l;
1595  }
const double sigma
Definition: SM_Z.hpp:43
void calc_lnL_Haloscope_RBF(double &result)
Definition: Axions.cpp:1557

◆ calc_lnL_Haloscope_UF()

void Gambit::DarkBit::calc_lnL_Haloscope_UF ( double result)

Definition at line 1533 of file Axions.cpp.

References Gambit::DecayBit::SM_Z::sigma.

1534  {
1535  using namespace Pipes::calc_lnL_Haloscope_UF;
1536  // Rescale the axion mass to ueV.
1537  double m = 1.0E+6*(*Param["ma0"]);
1538  double l = 0.0;
1539 
1540  // There are only limits between 5.4 and 5.9 ueV.
1541  if ((m > 5.4) && (m < 5.9)) {
1542  // Likelihood parameters based on information from Phys. Rev. D 42, 1297(R) (1990); correspond to power in 10^-22 W.
1543  const double sigma = 2.859772;
1544  // The "signal" needs to be rescaled to 0.2804 GeV/cm^3 (their reference value)
1545  // and also to the reference DFSZ coupling strength gDFSZ^2 = 0.6188 x 10^-30 GeV^-2
1546  // const double PowerDFSZ = 6.92;
1547  //double s = (0.45/0.2804)*(PowerDFSZ/0.6188)*(*Dep::Haloscope_signal);
1548  //double s = 0.035083106*(0.45/0.2804)*(*Dep::Haloscope_signal);
1549  double s = 0.0273012*(*Dep::Haloscope_signal);
1550  l = -0.5 * gsl_pow_2(s/sigma);
1551  }
1552 
1553  result = l;
1554  }
const double sigma
Definition: SM_Z.hpp:43
void calc_lnL_Haloscope_UF(double &result)
Definition: Axions.cpp:1533

◆ calc_lnL_HESS_GCMF()

void Gambit::DarkBit::calc_lnL_HESS_GCMF ( double result)

Definition at line 2095 of file Axions.cpp.

References daFunk::interp().

2096  {
2097  using namespace Pipes::calc_lnL_HESS_GCMF;
2098  double m_ax = *Param["ma0"];
2099  double gagg = gagg_conversion*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
2100 
2101  // Compute the domensionless parameters Epsilon and Gamma from the axion mass and axion-photon coupling (see 1311.3148).
2102  const double c_epsilon = 0.071546787;
2103  const double c_gamma = 0.015274036*370.0/sqrt(37.0);
2104  double epsilon = log10(m_ax*c_epsilon) + 5.0;
2105  double gamma = log10(gagg*c_gamma) + 20.0;
2106 
2107  // Initialise the interpolation and extrapolation routies for the H.E.S.S. results.
2108  static HESS_Interpolator interp (GAMBIT_DIR "/DarkBit/data/HESS_GCMF_Table.dat");
2109 
2110  result = interp.lnL(epsilon, gamma);
2111  }
const double gagg_conversion
Supporting classes and functions for the axion module.
Definition: Axions.cpp:70
MATH_OPERATION(Dif,-) MATH_OPERATION(pow) MATH_OPERATION(fmin) MATH_OPERATION(fmax) class FunkInterp shared_ptr< FunkInterp > interp(T f, std::vector< double > x, std::vector< double > y)
Definition: daFunk.hpp:1349
void calc_lnL_HESS_GCMF(double &result)
Definition: Axions.cpp:2095
Here is the call graph for this function:

◆ calc_lnL_RParameter()

void Gambit::DarkBit::calc_lnL_RParameter ( double result)

Definition at line 1899 of file Axions.cpp.

References linear, Gambit::pi, Gambit::Scanner::pow(), and x2.

1900  {
1901  using namespace Pipes::calc_lnL_RParameter;
1902  double Rtheo = *Dep::RParameter;
1903 
1904  // Observed R values from astro-ph/0403600.
1905  const double Robs = 1.39;
1906  const double RobsErr = 0.03;
1907  // Value for He-abundance Y from 1503.08146: <Y> = 0.2515(17).
1908  const double YobsErrContrib = 7.3306*0.0017;
1909 
1910  result = -0.5*gsl_pow_2(Rtheo - Robs)/(RobsErr*RobsErr+YobsErrContrib*YobsErrContrib);
1911  }
void calc_lnL_RParameter(double &result)
Definition: Axions.cpp:1899
Here is the call graph for this function:

◆ calc_lnL_SN1987A()

void Gambit::DarkBit::calc_lnL_SN1987A ( double result)

Definition at line 2061 of file Axions.cpp.

2062  {
2063  using namespace Pipes::calc_lnL_SN1987A;
2064  double f_10s = *Dep::PhotonFluence_SN1987A_Conversion;
2065 
2066  // Standard devations of the null observation.
2067  const double sigma_10s = 0.2;
2068 
2069  double ratio = f_10s/sigma_10s;
2070 
2071  result = -0.5*ratio*ratio;
2072  }
void calc_lnL_SN1987A(double &result)
Definition: Axions.cpp:2061

◆ calc_lnL_WDVar_G117B15A()

void Gambit::DarkBit::calc_lnL_WDVar_G117B15A ( double result)

Definition at line 1952 of file Axions.cpp.

References x2.

1953  {
1954  using namespace Pipes::calc_lnL_WDVar_G117B15A;
1955  // Rescale coupling to be used in their model prediction.
1956  double x2 = (1.0e+14 * std::fabs(*Param["gaee"]))/2.8;
1957  x2 = x2*x2;
1958 
1959  // Values for the model prediction provided by the authors.
1960  const std::vector<double> x2vals = {0.0, 1.0, 6.25, 25.0, 56.25, 100.0, 156.25, 225.0, 306.25, 404.0, 506.25, 625.0, 756.25, 900.0};
1961  const std::vector<double> dPidts = {1.235687, 1.244741, 1.299579, 1.470017, 1.796766, 2.260604, 2.795575, 3.484570, 4.232738, 5.056075, 6.113390, 7.342085, 8.344424, 9.775156};
1962  const double err2 = 0.09*0.09;
1963  const double obs = 4.19;
1964  const double obs_err2 = 0.73*0.73;
1965 
1966  static WDInterpolator dPidt (x2vals, dPidts, GAMBIT_DIR "/DarkBit/data/Axions_WDCorrection_G117B15A.dat");
1967 
1968  // Use interpolation for the finite-mass correction.
1969  const double internal_temperature_keV = 1.19698;
1970  double mrel = 0.001 * (*Param["ma0"]) / internal_temperature_keV;
1971 
1972  double pred = dPidt.evaluate(mrel, x2);
1973 
1974  result = -0.5 * gsl_pow_2(obs - pred) / (obs_err2 + err2);
1975  }
START_MODEL x2
Definition: demo.hpp:42
void calc_lnL_WDVar_G117B15A(double &result)
Definition: Axions.cpp:1952

◆ calc_lnL_WDVar_L192()

void Gambit::DarkBit::calc_lnL_WDVar_L192 ( double result)

Definition at line 2030 of file Axions.cpp.

References x2.

2031  {
2032  using namespace Pipes::calc_lnL_WDVar_L192;
2033  // Rescale coupling to be used in their model prediction.
2034  double x2 = (1.0E+14 * std::fabs(*Param["gaee"]))/2.8;
2035  x2 = x2*x2;
2036 
2037  // Values for the model prediction provided by the authors.
2038  const std::vector<double> x2vals = {0.0, 1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0, 121.0, 144.0, 169.0, 196.0, 225.0, 256.0, 289.0, 324.0, 361.0, 400.0, 441.0, 484.0, 529.0, 576.0, 625.0, 676.0, 729.0, 784.0, 841.0, 900.0};
2039  const std::vector<double> dPidts = {2.41, 2.40, 2.44, 2.42, 2.50, 2.57, 2.63, 2.74, 2.83, 2.99, 3.15, 3.32, 3.52, 3.70, 3.90, 4.08, 4.42, 4.69, 4.98, 5.34, 5.62, 6.02, 6.27, 6.62, 7.04, 7.38, 7.89, 8.09, 8.65, 9.16, 9.62};
2040  const double err2 = 0.85*0.85;
2041  const double obs = 3.0;
2042  const double obs_err2 = 0.6*0.6;
2043 
2044  static WDInterpolator dPidt (x2vals, dPidts, GAMBIT_DIR "/DarkBit/data/Axions_WDCorrection_L192.dat");
2045 
2046  // Use interpolation for the finite-mass correction.
2047  const double internal_temperature_keV = 1.04931;
2048  double mrel = 0.001 * (*Param["ma0"]) / internal_temperature_keV;
2049 
2050  double pred = dPidt.evaluate(mrel, x2);
2051 
2052  result = -0.5 * gsl_pow_2(obs - pred) / (obs_err2 + err2);
2053  }
void calc_lnL_WDVar_L192(double &result)
Definition: Axions.cpp:2030
START_MODEL x2
Definition: demo.hpp:42

◆ calc_lnL_WDVar_PG1351489()

void Gambit::DarkBit::calc_lnL_WDVar_PG1351489 ( double result)

Definition at line 2004 of file Axions.cpp.

References x2.

2005  {
2006  using namespace Pipes::calc_lnL_WDVar_PG1351489;
2007  // Rescale coupling to be used in their model prediction.
2008  double x2 = (1.0E+14 * std::fabs(*Param["gaee"]))/2.8;
2009  x2 = x2*x2;
2010 
2011  // Values for the model prediction provided by the authors.
2012  const std::vector<double> x2vals = {0.0, 4.0, 16.0, 36.0, 64.0, 100.0, 144.0, 196.0, 256.0, 324.0, 400.0};
2013  const std::vector<double> dPidts = {0.90878126, 0.96382008, 1.2022906, 1.5712931, 2.1220619, 2.8002354, 3.6172605, 4.5000560, 5.5256592, 6.5055283, 7.5341296};
2014  const double err2 = 0.5*0.5;
2015  const double obs = 2.0;
2016  const double obs_err2 = 0.9*0.9;
2017 
2018  static WDInterpolator dPidt (x2vals, dPidts, GAMBIT_DIR "/DarkBit/data/Axions_WDCorrection_PG1351489.dat");
2019 
2020  // Use interpolation for the finite-mass correction.
2021  const double internal_temperature_keV = 2.64273;
2022  double mrel = 0.001 * (*Param["ma0"]) / internal_temperature_keV;
2023 
2024  double pred = dPidt.evaluate(mrel, x2);
2025 
2026  result = -0.5 * gsl_pow_2(obs - pred) / (obs_err2 + err2);
2027  }
START_MODEL x2
Definition: demo.hpp:42
void calc_lnL_WDVar_PG1351489(double &result)
Definition: Axions.cpp:2004

◆ calc_lnL_WDVar_R548()

void Gambit::DarkBit::calc_lnL_WDVar_R548 ( double result)

Definition at line 1978 of file Axions.cpp.

References x2.

1979  {
1980  using namespace Pipes::calc_lnL_WDVar_R548;
1981  // Rescale coupling to be used in their model prediction.
1982  double x2 = (1.0E+14 * std::fabs(*Param["gaee"]))/2.8;
1983  x2 = x2*x2;
1984 
1985  // Values for the model prediction provided by the authors.
1986  const std::vector<double> x2vals = {0.0, 1.0, 6.25, 25.0, 56.25, 100.0, 156.25, 225.0, 306.25, 400.0, 506.25, 625.0, 756.25, 900.0};
1987  const std::vector<double> dPidts = {1.075373, 1.095319, 1.123040, 1.289434, 1.497666, 1.869437, 2.300523, 2.844954, 3.379978, 4.086028, 4.847149, 5.754807, 6.714841, 7.649140};
1988  const double err2 = 0.09*0.09;
1989  const double obs = 3.3;
1990  const double obs_err2 = 1.1*1.1;
1991 
1992  static WDInterpolator dPidt (x2vals, dPidts, GAMBIT_DIR "/DarkBit/data/Axions_WDCorrection_R548.dat");
1993 
1994  // Use interpolation for the finite-mass correction.
1995  const double internal_temperature_keV = 1.11447;
1996  double mrel = 0.001 * (*Param["ma0"]) / internal_temperature_keV;
1997 
1998  double pred = dPidt.evaluate(mrel, x2);
1999 
2000  result = -0.5 * gsl_pow_2(obs - pred) / (obs_err2 + err2);
2001  }
START_MODEL x2
Definition: demo.hpp:42
void calc_lnL_WDVar_R548(double &result)
Definition: Axions.cpp:1978

◆ calc_lnL_XENON1T_Anomaly()

void Gambit::DarkBit::calc_lnL_XENON1T_Anomaly ( double result)

Capability for the XENON1T likelihood from 2006.10035.

The signal model consists of 3 components: Primakoff, ABC, and Fe57.

Definition at line 2118 of file Axions.cpp.

References Gambit::DecayBit::SM_Z::sigma.

2119  {
2120  using namespace Pipes::calc_lnL_XENON1T_Anomaly;
2121 
2122  // Rescale couplings and 3H abundance to the reference values used in 2006.10035 for convenience.
2123  double gae = std::fabs(*Param["gaee"]) / 5.0e-12;
2124  double gagamma = std::fabs(*Param["gagg"]) / 2.0e-10;
2125  double gaN = std::fabs(*Param["gaN"]) / 1.0e-6;
2126  double x_3H = *Param["x_3H"] / 6.2e-25;
2127  double bkg_scale = 1.0 + *Param["delta_bkg"];
2128  double eff = 1.0 + *Param["delta_eff"];
2129 
2130  static bool include_inverse_Primakoff = runOptions->getValueOrDef<bool> (true, "include_inverse_Primakoff");
2131 
2132  // XENON1T 2020 data (based on 2006.10035 and using an exposure of 0.65 tonne-years)
2133  static const Eigen::ArrayXd observed = (Eigen::ArrayXd(29) <<
2134  26., 61., 55., 47., 49.,
2135  47., 44., 41., 40., 37.,
2136  51., 41., 42., 51., 47.,
2137  48., 24., 43., 42., 34.,
2138  42., 40., 38., 53., 41.,
2139  57., 39., 46., 35.).finished();
2140 
2141  static const Eigen::ArrayXd bkg_tritium = (Eigen::ArrayXd(29) <<
2142  4.54543769e+00, 8.60509728e+00, 8.94256482e+00, 8.61684767e+00, 8.02464466e+00,
2143  7.29168481e+00, 6.48068011e+00, 5.65037508e+00, 4.81438779e+00, 3.97835836e+00,
2144  3.17827210e+00, 2.43987288e+00, 1.78022901e+00, 1.21426641e+00, 7.54437371e-01,
2145  4.13276721e-01, 1.95253807e-01, 7.58276424e-02, 2.37741495e-02, 1.17262241e-02,
2146  4.76851304e-03, 1.65434695e-04, 5.42752619e-05, 5.22947241e-05, 5.03141863e-05,
2147  4.83336485e-05, 4.63531107e-05, 4.43725729e-05, 4.23920351e-05).finished();
2148 
2149  static const Eigen::ArrayXd bkg_other = (Eigen::ArrayXd(29) <<
2150  22.07404375, 39.45186703, 41.83417331, 42.46968003, 42.78761694, 42.96532578,
2151  43.13847573, 43.56505579, 44.1162301 , 44.04330642, 43.60594679, 43.40248223,
2152  43.45031774, 43.49263918, 43.57084528, 43.66270961, 43.75990478, 43.85453193,
2153  43.95159076, 44.04782058, 44.14510208, 44.24450247, 44.3406822 , 44.43638726,
2154  44.5331988, 44.62865958, 44.72654689, 44.82382807, 44.91842725).finished();
2155 
2156  static const Eigen::ArrayXd signal_ref_ABC = (Eigen::ArrayXd(29) <<
2157  7.46283683e+01, 9.11300502e+01, 4.91874199e+01, 3.53433982e+01, 3.97196350e+01,
2158  3.57128137e+01, 2.27540737e+01, 1.19536450e+01, 6.29278747e+00, 3.30948412e+00,
2159  1.61495065e+00, 6.21479171e-01, 1.55434142e-01, 2.13046184e-02, 1.63362970e-03,
2160  5.94631877e-05, 6.22346656e-07, 0, 0, 0,
2161  0, 0, 0, 0, 0,
2162  0, 0, 0, 0).finished();
2163 
2164  static const Eigen::ArrayXd signal_ref_primakoff = (Eigen::ArrayXd(29) <<
2165  8.52269995e+00, 2.00539997e+01, 1.92433543e+01, 2.22113223e+01, 2.89487211e+01,
2166  2.48618807e+01, 1.56554086e+01, 8.76011034e+00, 4.77088790e+00, 2.59966792e+00,
2167  1.43054774e+00, 7.84191139e-01, 4.15162321e-01, 2.10029374e-01, 1.02847245e-01,
2168  5.27333566e-02, 2.93001979e-02, 1.70759274e-02, 1.08781046e-02, 7.76743790e-03,
2169  6.24916022e-03, 5.51837192e-03, 5.15693524e-03, 4.97968515e-03, 4.88717867e-03,
2170  4.84242242e-03, 4.82038071e-03, 2.38442778e-03, 0).finished();
2171 
2172  static const Eigen::ArrayXd signal_ref_fe57 = (Eigen::ArrayXd(29) <<
2173  4.64697658e-22, 1.14417236e-18, 1.48421451e-15, 1.01555196e-12, 3.67061998e-10,
2174  7.02067458e-08, 7.12116995e-06, 3.84028003e-04, 1.10431167e-02, 1.69885040e-01,
2175  1.40292258e+00, 6.23942228e+00, 1.49856121e+01, 1.94718769e+01, 1.36959639e+01,
2176  5.21070939e+00, 1.07020697e+00, 1.18324090e-01, 7.01902932e-03, 2.22639151e-04,
2177  3.76396817e-06, 3.38186408e-08, 1.61083936e-10, 4.05908328e-13, 5.40175385e-16,
2178  3.79105023e-19, 1.40152641e-22, 2.72678128e-26, 2.78978087e-30).finished();
2179 
2180  static const Eigen::ArrayXd signal_ref_ABC_inv_p = (Eigen::ArrayXd(29) <<
2181  6.90095435e+00, 1.39022260e+01, 1.15151093e+01, 7.97641672e+00, 5.52073353e+00,
2182  4.15542096e+00, 2.85016524e+00, 1.77391643e+00, 1.07580666e+00, 6.28453363e-01,
2183  3.26090878e-01, 1.31638303e-01, 3.77207469e-02, 7.52278937e-03, 1.06834615e-03,
2184  1.12228106e-04, 9.07094259e-06, 5.85004243e-07, 3.10665723e-08, 1.39581067e-09,
2185  5.42752207e-11, 1.86190569e-12, 5.72719620e-14, 1.60150437e-15, 4.11906863e-17,
2186  9.84235043e-19, 2.20372014e-20, 4.65786516e-22, 9.35349973e-24).finished();
2187 
2188  static const Eigen::ArrayXd signal_ref_primakoff_inv_p = (Eigen::ArrayXd(29) <<
2189  8.80070525e-01, 3.29923191e+00, 4.56900949e+00, 4.56444853e+00, 3.82216381e+00,
2190  2.86024934e+00, 1.98175269e+00, 1.29748018e+00, 8.13307028e-01, 4.92544613e-01,
2191  2.90070188e-01, 1.66927792e-01, 9.42159402e-02, 5.23046562e-02, 2.86265445e-02,
2192  1.54743414e-02, 8.27421198e-03, 4.38184245e-03, 2.30069622e-03, 1.19872791e-03,
2193  6.20254770e-04, 3.18925270e-04, 1.63051423e-04, 8.29261227e-05, 4.19736901e-05,
2194  2.11518187e-05, 1.06157318e-05, 5.30780091e-06, 2.64457315e-06).finished();
2195 
2196  static const Eigen::ArrayXd signal_ref_fe57_inv_p = (Eigen::ArrayXd(29) <<
2197  2.68534579e-23, 1.02403478e-19, 1.64401745e-16, 1.34346994e-13, 5.65333857e-11,
2198  1.23357889e-08, 1.40118942e-06, 8.30825834e-05, 2.57975431e-03, 4.20954642e-02,
2199  3.62322333e-01, 1.65087221e+00, 3.99391811e+00, 5.14060546e+00, 3.52226638e+00,
2200  1.28362583e+00, 2.48258237e-01, 2.54009765e-02, 1.36995132e-03, 3.88033329e-05,
2201  5.75235843e-07, 4.44951177e-09, 1.79118821e-11, 3.74451860e-14, 4.05797696e-17,
2202  2.27644923e-20, 6.60289329e-24, 9.89294778e-28, 7.65058090e-32).finished();
2203 
2204  static const double asimov = (observed * observed.log() - observed).sum();
2205 
2206  const Eigen::ArrayXd bkg = x_3H * bkg_tritium + bkg_other;
2207  Eigen::ArrayXd signal = gae * gae * (
2208  signal_ref_ABC * gae * gae +
2209  signal_ref_primakoff * gagamma * gagamma +
2210  signal_ref_fe57 * gaN * gaN);
2211 
2212  if (include_inverse_Primakoff)
2213  {
2214  signal = signal + gagamma * gagamma * (
2215  signal_ref_ABC_inv_p * gae * gae +
2216  signal_ref_primakoff_inv_p * gagamma * gagamma +
2217  signal_ref_fe57_inv_p * gaN * gaN);
2218  }
2219 
2220  const Eigen::ArrayXd expected = eff * (bkg_scale * bkg + signal);
2221 
2222  result = (observed * expected.log() - expected).sum() - asimov;
2223  }
void calc_lnL_XENON1T_Anomaly(double &result)
Capability for the XENON1T likelihood from 2006.10035.
Definition: Axions.cpp:2118

◆ calc_lnL_XENON1T_Anomaly_NuisanceParameters()

void Gambit::DarkBit::calc_lnL_XENON1T_Anomaly_NuisanceParameters ( double result)

Definition at line 2314 of file Axions.cpp.

2315  {
2317 
2318  result = -0.5 * ( gsl_pow_2(*Param["delta_bkg"]/0.026) + gsl_pow_2(*Param["delta_eff"]/0.030) );
2319  }
void calc_lnL_XENON1T_Anomaly_NuisanceParameters(double &result)
Definition: Axions.cpp:2314

◆ calc_lnL_XENON1T_DM_Anomaly()

void Gambit::DarkBit::calc_lnL_XENON1T_DM_Anomaly ( double result)

Definition at line 2236 of file Axions.cpp.

References dRdE(), f, Gambit::pi, Gambit::LocalMaxwellianHalo::rho0, and Gambit::DecayBit::SM_Z::sigma.

2237  {
2238  using namespace Pipes::calc_lnL_XENON1T_DM_Anomaly;
2239 
2240  result = 0;
2241 
2242  // Rescale couplings and 3H abundance to the reference values used in 2006.10035 for convenience.
2243  double gae = std::fabs(*Param["gaee"]);
2244  double ma = *Param["ma0"] / 1.0e3;
2245  double x_3H = *Param["x_3H"] / 6.2e-25;
2246  double bkg_scale = 1.0 + *Param["delta_bkg"];
2247  double rel_eff = 1.0 + *Param["delta_eff"];
2248  double dm_fraction = *Param["eta"];
2249  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
2250  double rho0 = LocalHaloParameters.rho0;
2251 
2252  // XENON1T 2020 data (based on 2006.10035 and using an exposure of 0.65 tonne-years)
2253  static const Eigen::ArrayXd observed = (Eigen::ArrayXd(29) <<
2254  26., 61., 55., 47., 49.,
2255  47., 44., 41., 40., 37.,
2256  51., 41., 42., 51., 47.,
2257  48., 24., 43., 42., 34.,
2258  42., 40., 38., 53., 41.,
2259  57., 39., 46., 35.).finished();
2260 
2261  static const Eigen::ArrayXd bkg_tritium = (Eigen::ArrayXd(29) <<
2262  4.54543769e+00, 8.60509728e+00, 8.94256482e+00, 8.61684767e+00, 8.02464466e+00,
2263  7.29168481e+00, 6.48068011e+00, 5.65037508e+00, 4.81438779e+00, 3.97835836e+00,
2264  3.17827210e+00, 2.43987288e+00, 1.78022901e+00, 1.21426641e+00, 7.54437371e-01,
2265  4.13276721e-01, 1.95253807e-01, 7.58276424e-02, 2.37741495e-02, 1.17262241e-02,
2266  4.76851304e-03, 1.65434695e-04, 5.42752619e-05, 5.22947241e-05, 5.03141863e-05,
2267  4.83336485e-05, 4.63531107e-05, 4.43725729e-05, 4.23920351e-05).finished();
2268 
2269  static const Eigen::ArrayXd bkg_other = (Eigen::ArrayXd(29) <<
2270  22.07404375, 39.45186703, 41.83417331, 42.46968003, 42.78761694,
2271  42.96532578, 43.13847573, 43.56505579, 44.1162301 , 44.04330642,
2272  43.60594679, 43.40248223, 43.45031774, 43.49263918, 43.57084528,
2273  43.66270961, 43.75990478, 43.85453193, 43.95159076, 44.04782058,
2274  44.14510208, 44.24450247, 44.3406822 , 44.43638726, 44.5331988 ,
2275  44.62865958, 44.72654689, 44.82382807, 44.91842725).finished();
2276 
2277  // Photoelectric cross section for Xe from https://dx.doi.org/10.18434/T48G6X
2278  // Columns: Photon energy [MeV] | Photoelectric absorption [10^-28 m^2/atom]
2279  static AxionInterpolator sigma_pe (GAMBIT_DIR "/DarkBit/data/XENON1T/photoelectric.txt");
2280 
2281  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
2282  if ( (ma >= 1.0) && (ma <= 30.0) )
2283  {
2284  // Energy resolution from 2003.03825
2285  double energy_resolution = 0.15 + 31.71/sqrt(ma);
2286  double sigma = ma * energy_resolution / 100.0;
2287  const double exposure = 0.647309514*1000.0*365.0;
2288  const double photoel_eff_conversion = 1.5e19/131.0;
2289  double amplitude = dm_fraction * (rho0/0.3) * exposure * gae*gae * ma * photoel_eff_conversion*sigma_pe.interpolate(ma/1000.0);
2290  gsl_function f;
2291  struct dRdE_params params = {ma, sigma};
2292  f.function = &dRdE;
2293  f.params = &params;
2294  double dRdE_result, error;
2295  std::vector<double> signal_vec;
2296  for (int i = 0; i < 29; i++)
2297  {
2298  gsl_integration_qags (&f, 1.+i, 2.+i, 0, 1e-7, 1000, w, &dRdE_result, &error);
2299  double s = amplitude * 1/sqrt(2*pi)/sigma * dRdE_result;
2300  signal_vec.push_back(s);
2301  }
2302  gsl_integration_workspace_free (w);
2303 
2304  static const double asimov = (observed * observed.log() - observed).sum();
2305 
2306  const Eigen::ArrayXd bkg = x_3H * bkg_tritium + bkg_other;
2307  const Eigen::ArrayXd signal = Eigen::ArrayXd::Map(signal_vec.data(), signal_vec.size());
2308  const Eigen::ArrayXd expected = rel_eff * (bkg_scale * bkg + signal);
2309 
2310  result = (observed * expected.log() - expected).sum() - asimov;
2311  }
2312  }
double dRdE(double E, void *params)
Definition: Axions.cpp:2227
const double sigma
Definition: SM_Z.hpp:43
const double pi
void calc_lnL_XENON1T_DM_Anomaly(double &result)
Definition: Axions.cpp:2236
Here is the call graph for this function:

◆ calc_PhotonFluence_SN1987A_Conversion()

void Gambit::DarkBit::calc_PhotonFluence_SN1987A_Conversion ( double result)

Definition at line 2080 of file Axions.cpp.

References g, and Gambit::Scanner::pow().

2081  {
2083  double m = (1.0E+10*(*Param["ma0"]))/5.433430;
2084  double g = (1.0E+12*std::fabs(*Param["gagg"]))/5.339450;
2085 
2086  result = 0.570589*gsl_pow_4(g);
2087  if (m > 1.0) { result = result*pow(m, -4.021046); }
2088  }
double pow(const double &a)
Outputs a^i.
void calc_PhotonFluence_SN1987A_Conversion(double &result)
Definition: Axions.cpp:2080
Here is the call graph for this function:

◆ calc_RParameter()

void Gambit::DarkBit::calc_RParameter ( double result)

Capabilities relating to astrophysical observations (R-parameter, H.E.S.S. telescope search, cooling hints).

Supported models: GeneralALP

Definition at line 1873 of file Axions.cpp.

References linear, and Gambit::Scanner::pow().

1874  {
1875  using namespace Pipes::calc_RParameter;
1876  double ma0 = *Param["ma0"];
1877  double gaee2 = gsl_pow_2(gaee_conversion*std::fabs(*Param["gaee"]));
1878  double gagg = 1.0E+10*std::fabs(*Param["gagg"]); // gagg needs to be in 10^-10 GeV^-1.
1879  // Value for He-abundance Y from 1503.08146: <Y> = 0.2515(17).
1880  const double Y = 0.2515;
1881  // Use interpolation for the finite-mass correction.
1882  static AxionInterpolator correction (GAMBIT_DIR "/DarkBit/data/Axions_RParameterCorrection.dat", InterpolationOptions1D::linear);
1883  // Initialise an effective axion-photon coupling, valid for low masses.
1884  double geff = gagg;
1885  // Apply correction for higher mass values...
1886  static double m_min = pow(10,correction.lower());
1887  static double m_max = pow(10,correction.upper());
1888  if ((ma0 > m_min) && (ma0 < m_max)) { geff *= pow(10, 0.5*correction.interpolate(log10(ma0))); }
1889  // ... or set to zero if mass is too high.
1890  if (ma0 >= m_max) { geff = 0.0; }
1891  // Expressions only valid for gaee2 < 35.18 but limits should become stronger for gaee2 > 35.18 (but perhaps not gaee2 >> 35.18).
1892  // Conservative approach: Constrain gaee2 > 35.18 at the level of gaee2 = 35.18.
1893  if (gaee2 > 35.18) { gaee2 = 35.18; }
1894 
1895  result = -0.421824 - 0.0948659*(-4.675 + sqrt(21.8556 + 21.0824*geff)) - 0.00533169*gaee2 - 0.0386834*(-1.23 - 0.137991*pow(gaee2,0.75) + sqrt(1.5129 + gaee2)) + 7.3306*Y;
1896  }
void calc_RParameter(double &result)
Capabilities relating to astrophysical observations (R-parameter, H.E.S.S. telescope search...
Definition: Axions.cpp:1873
const double gaee_conversion
Definition: Axions.cpp:71
double pow(const double &a)
Outputs a^i.
Here is the call graph for this function:

◆ capture_rate_Sun_const_xsec()

void Gambit::DarkBit::capture_rate_Sun_const_xsec ( double result)

Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1). DarkSUSY 6 version.

Definition at line 66 of file SunNeutrinos.cpp.

References DarkBit_error(), LOCAL_INFO, and Gambit::LocalMaxwellianHalo::rho0.

Referenced by main().

67  {
68  using namespace Pipes::capture_rate_Sun_const_xsec;
69 
70  if (BEreq::cap_Sun_v0q0_isoscalar.origin()=="DarkSUSY")
71  if(!(*Dep::DarkSUSY_PointInit_LocalHalo))
72  DarkBit_error().raise(LOCAL_INFO,"DarkSUSY halo model not initialized!");
73 
74  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
75  double rho0 = LocalHaloParameters.rho0;
76  double rho0_eff = (*Dep::RD_fraction)*rho0;
77 
78  // When calculating the solar capture rate, DarkSUSY assumes that the
79  // proton and neutron scattering cross-sections are the same; we
80  // assume that whichever backend has been hooked up here does so too.
81 
82  result = BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp, rho0_eff, *Dep::sigma_SI_p, *Dep::sigma_SD_p);
83 
84  }
error & DarkBit_error()
#define LOCAL_INFO
Definition: local_info.hpp:34
void capture_rate_Sun_const_xsec(double &result)
Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1)...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ capture_rate_Sun_const_xsec_capgen()

void Gambit::DarkBit::capture_rate_Sun_const_xsec_capgen ( double result)

Alternative to the darkSusy fct, using captn_specific from capgen instead.

Definition at line 87 of file SunNeutrinos.cpp.

88  {
90  double resultSD;
91  double resultSI;
92  double maxcap;
93 
94  BEreq::cap_sun_saturation(*Dep::mwimp,maxcap);
95  BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp,*Dep::sigma_SD_p,*Dep::sigma_SI_p,resultSD,resultSI);
96  result = resultSI + resultSD;
97 
98  if (maxcap < result)
99  {
100  result = maxcap;
101  }
102 
103  }
void capture_rate_Sun_const_xsec_capgen(double &result)
Alternative to the darkSusy fct, using captn_specific from capgen instead.

◆ capture_rate_Sun_const_xsec_DS5()

void Gambit::DarkBit::capture_rate_Sun_const_xsec_DS5 ( double result)

Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1). DarkSUSY 5 version.

Definition at line 47 of file SunNeutrinos.cpp.

References DarkBit_error(), and LOCAL_INFO.

Referenced by main().

48  {
50 
51  if (BEreq::cap_Sun_v0q0_isoscalar.origin()=="DarkSUSY")
52  if(!(*Dep::DarkSUSY5_PointInit_LocalHalo))
53  DarkBit_error().raise(LOCAL_INFO,"DarkSUSY halo model not initialized!");
54 
55  // When calculating the solar capture rate, DarkSUSY assumes that the
56  // proton and neutron scattering cross-sections are the same; we
57  // assume that whichever backend has been hooked up here does so too.
58 
59  result = BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp, *Dep::sigma_SI_p, *Dep::sigma_SD_p);
60 
61  }
error & DarkBit_error()
#define LOCAL_INFO
Definition: local_info.hpp:34
void capture_rate_Sun_const_xsec_DS5(double &result)
Capture rate of regular dark matter in the Sun (no v-dependent or q-dependent cross-sections) (s^-1)...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ capture_rate_Sun_vnqn()

void Gambit::DarkBit::capture_rate_Sun_vnqn ( double result)

Capture rate for v^n and q^n-dependent cross sections.

Isoscalar (same proton/neutron coupling) SD only couples to Hydrogen. See DirectDetection.cpp to see how to define the cross sections sigma_SD_p, sigma_SI_pi

Definition at line 109 of file SunNeutrinos.cpp.

References Gambit::EOM, and Gambit::logger().

110  {
111  using namespace Pipes::capture_rate_Sun_vnqn;
112 
113  double resultSD;
114  double resultSI;
115  double capped;
116  int qpow;
117  int vpow;
118  const int nelems = 29;
119  double maxcap;
120 
121  BEreq::cap_sun_saturation(*Dep::mwimp,maxcap);
122 
123  resultSI = 0e0;
124  resultSD = 0e0;
125  typedef map_intpair_dbl::const_iterator it_type;
126 
127  //Spin-dependent:
128  for(it_type iterator = Dep::sigma_SD_p->begin();
129  iterator != Dep::sigma_SD_p->end();
130  iterator++)
131  {
132  //don't capture anything if cross section is zero or all the DM is already capped
133  if((iterator->second > 1e-90) && (resultSD < maxcap))
134  {
135  qpow = (iterator->first).first/2 ;
136  vpow = (iterator->first).second/2;
137 
138  //Capture
139  BEreq::cap_Sun_vnqn_isoscalar(*Dep::mwimp,iterator->second,1,qpow,vpow,capped);
140  resultSD = resultSD+capped;
141  }
142  }
143 
144  //Spin independent:
145  for(it_type iterator = Dep::sigma_SI_p->begin();
146  iterator != Dep::sigma_SI_p->end();
147  iterator++)
148  {
149  if((iterator->second > 1e-90) && (resultSI+resultSD < maxcap))
150  {
151  qpow = (iterator->first).first/2 ;
152  vpow = (iterator->first).second/2;
153 
154  //Capture
155  BEreq::cap_Sun_vnqn_isoscalar(*Dep::mwimp,iterator->second,nelems,qpow,vpow,capped);
156  resultSI = resultSI+capped;
157  }
158  }
159  result = resultSI+resultSD;
160 
161  logger() << "Capgen captured: SI: " << resultSI << " SD: " << resultSD << " total: " << result << "max = " << maxcap << "\n" << EOM;
162 
163  // If capture is above saturation, return saturation value.
164  if (maxcap < result)
165  {
166  result = maxcap;
167  }
168 
169  //cout << "capture rate via capture_rate_Sun_vnqn = " << result << "\n";
170 
171  }
void capture_rate_Sun_vnqn(double &result)
Capture rate for v^n and q^n-dependent cross sections.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
Here is the call graph for this function:

◆ cascadeMC_DecayTable()

void Gambit::DarkBit::cascadeMC_DecayTable ( DarkBit::DecayChain::DecayTable table)

Function setting up the decay table used in decay chains.

Definition at line 58 of file Cascades.cpp.

References DarkBit_error(), Gambit::LogTags::err, and Gambit::DarkBit::DecayChain::DecayTable::printTable().

Referenced by cascadeMC_GenerateChain(), dumpSpectrum(), and main().

59  {
60  using namespace DecayChain;
61  using namespace Pipes::cascadeMC_DecayTable;
62  std::set<std::string> disabled;
63  // Note: One could add to "disabled" particles decays in that are in the
64  // process catalog but should for some reason not propagate to the FCMC
65  // DecayTable.
66  try
67  {
68  table = DecayTable(*Dep::TH_ProcessCatalog, *Dep::SimYieldTable, disabled);
69  }
70  catch(Piped_exceptions::description err)
71  {
72  DarkBit_error().raise(err.first,err.second);
73  }
74 #ifdef DARKBIT_DEBUG
75  table.printTable();
76 #endif
77  }
error & DarkBit_error()
void cascadeMC_DecayTable(DarkBit::DecayChain::DecayTable &table)
Function setting up the decay table used in decay chains.
Definition: Cascades.cpp:58
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cascadeMC_EventCount()

void Gambit::DarkBit::cascadeMC_EventCount ( std::map< std::string, int > &  counts)

Event counter for cascade decays.

Definition at line 175 of file Cascades.cpp.

References cascadeMC_InitialState(), MC_FINALIZE, MC_INIT, and MC_NEXT_STATE.

Referenced by cascadeMC_gammaSpectra(), and main().

176  {
177  using namespace Pipes::cascadeMC_EventCount;
178  static std::map<std::string, int> counters;
179  switch(*Loop::iteration)
180  {
181  case MC_INIT:
182  counters.clear();
183  break;
184  case MC_NEXT_STATE:
185  counters[*Dep::cascadeMC_InitialState] = 0;
186  break;
187  case MC_FINALIZE:
188  // For performance, only return the actual result once finished
189  counts=counters;
190  break;
191  default:
192 #pragma omp atomic
193  counters[*Dep::cascadeMC_InitialState]++;
194  }
195  }
void cascadeMC_InitialState(std::string &pID)
Function selecting initial state for decay chain.
Definition: Cascades.cpp:140
void cascadeMC_EventCount(std::map< std::string, int > &counts)
Event counter for cascade decays.
Definition: Cascades.cpp:175
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cascadeMC_fetchSpectra()

void Gambit::DarkBit::cascadeMC_fetchSpectra ( std::map< std::string, daFunk::Funk > &  spectra,
std::string  finalState,
const std::vector< std::string > &  ini,
const std::vector< std::string > &  fin,
const std::map< std::string, std::map< std::string, SimpleHist > > &  h,
const std::map< std::string, int > &  eventCounts 
)

Convenience function for getting a daFunk::Funk object of a given spectrum.

This function has no associated capability. Function retrieving specific spectra (like cascadeMC_gammaSpectra) should call this function.

Definition at line 625 of file Cascades.cpp.

References Gambit::DarkBit::SimpleHist::divideByBinSize(), Gambit::FlavBit::LoopFunctions::E(), Gambit::DarkBit::SimpleHist::getBinCenters(), Gambit::DarkBit::SimpleHist::getBinValues(), and daFunk::zero().

Referenced by cascadeMC_gammaSpectra().

631  {
632  spectra.clear();
633  // Check if final state has been calculated
634  bool calculated = (
635  std::find(fin.begin(), fin.end(), finalState) != fin.end());
636  for(std::vector<std::string>::const_iterator it = ini.begin();
637  it != ini.end(); ++it )
638  {
639 #ifdef DARKBIT_DEBUG
640  std::cout << "Trying to get cascade spectra for initial state: " << *it << std::endl;
641 #endif
642  if(calculated)
643  {
644 #ifdef DARKBIT_DEBUG
645  std::cout << finalState << "...was calculated!" << std::endl;
646  std::cout << eventCounts.at(*it) << " events generated" << std::endl;
647 #endif
648  SimpleHist hist = h.at(*it).at(finalState);
649  hist.divideByBinSize();
650  std::vector<double> E = hist.getBinCenters();
651  std::vector<double> dN_dE = hist.getBinValues();
652  // Normalize to per-event spectrum
653  int i = 0;
654  for (std::vector<double>::iterator it2=dN_dE.begin();
655  it2!=dN_dE.end();++it2)
656  {
657  *it2 /= eventCounts.at(*it);
658 #ifdef DARKBIT_DEBUG
659  std::cout << E[i] << " " << *it2 << std::endl;
660 #endif
661  i++;
662  }
663  // Default values provide 1-2% accuracy for singular integrals
664  // Make this optional.
665  spectra[*it] = daFunk::Funk(new daFunk::FunkInterp("E", E, dN_dE, "lin"));
666 
667  for (size_t i = 1; i<E.size()-1; i++)
668  {
669  if (dN_dE[i]/(dN_dE[i-1]+dN_dE[i+1]+dN_dE[i]*1e-4) > 1e2)
670  {
671 #ifdef DARKBIT_DEBUG
672  std::cout << "Set singularity at " << E[i] << " with width "
673  << E[i+1]-E[i] << endl;
674 #endif
675  spectra[*it]->set_singularity("E", E[i], (E[i+1]-E[i]));
676  }
677  }
678  }
679  else
680  {
681  spectra[*it] = daFunk::zero("E");
682  }
683  }
684  }
Funk zero(Args... argss)
Definition: daFunk.hpp:604
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
double E(const double x, const double y)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cascadeMC_FinalStates()

void Gambit::DarkBit::cascadeMC_FinalStates ( std::vector< std::string > &  list)

Function for retrieving list of final states for cascade decays.

Option cMC_finalStates<std::vector<std::string>>: List of final states to simulate (default is ["gamma"])

Definition at line 41 of file Cascades.cpp.

Referenced by cascadeMC_gammaSpectra(), cascadeMC_Histograms(), dumpSpectrum(), and main().

42  {
43  list.clear();
44  list.push_back("gamma");
45  using namespace Pipes::cascadeMC_FinalStates;
47  list = runOptions->getValueOrDef<std::vector<std::string>>(list, "cMC_finalStates");
48  #ifdef DARKBIT_DEBUG
49  std::cout << "Final states to generate: " << list.size() << std::endl;
50  for(size_t i=0; i < list.size(); i++)
51  {
52  std::cout << " " << list[i] << std::endl;
53  }
54  #endif
55  }
void cascadeMC_FinalStates(std::vector< std::string > &list)
Function for retrieving list of final states for cascade decays.
Definition: Cascades.cpp:41
Here is the caller graph for this function:

◆ cascadeMC_gammaSpectra()

void Gambit::DarkBit::cascadeMC_gammaSpectra ( std::map< std::string, daFunk::Funk > &  spectra)

Function requesting and returning gamma ray spectra from cascade decays.

Definition at line 687 of file Cascades.cpp.

References cascadeMC_EventCount(), cascadeMC_fetchSpectra(), cascadeMC_FinalStates(), cascadeMC_gammaSpectra, cascadeMC_Histograms(), Gambit::FlavBit::LoopFunctions::E(), f, and GA_missingFinalStates().

688  {
689  using namespace Pipes::cascadeMC_gammaSpectra;
693 #ifdef DARKBIT_DEBUG
694  std::cout << "Retrieving cascade spectra for gamma final states" << std::endl;
695  std::cout << "Number of simulated final states: " << spectra.size() << std::endl;
696  for ( auto it = spectra.begin(); it != spectra.end(); it ++ )
697  {
698  std::cout << "Particle: " << it->first << std::endl;
699  auto f= it->second;
700  for ( double E = 0.1; E < 1000; E*=1.5 )
701  {
702  std::cout << " " << E << " " << f->bind("E")->eval(E) << std::endl;
703  }
704  std::cout << " Integrated spectrum: " << f->gsl_integration("E", 0, 1000)->bind()->eval() << std::endl;
705  }
706 #endif
707  }
void cascadeMC_FinalStates(std::vector< std::string > &list)
Function for retrieving list of final states for cascade decays.
Definition: Cascades.cpp:41
void GA_missingFinalStates(std::vector< std::string > &result)
Identification of final states that are not yet tabulated.
Definition: GamYields.cpp:61
void cascadeMC_fetchSpectra(std::map< std::string, daFunk::Funk > &spectra, std::string finalState, const std::vector< std::string > &ini, const std::vector< std::string > &fin, const std::map< std::string, std::map< std::string, SimpleHist > > &h, const std::map< std::string, int > &eventCounts)
Convenience function for getting a daFunk::Funk object of a given spectrum.
Definition: Cascades.cpp:625
void cascadeMC_EventCount(std::map< std::string, int > &counts)
Event counter for cascade decays.
Definition: Cascades.cpp:175
void cascadeMC_Histograms(std::map< std::string, std::map< std::string, SimpleHist > > &result)
Function responsible for histogramming, and evaluating end conditions for event loop.
Definition: Cascades.cpp:363
double E(const double x, const double y)
Here is the call graph for this function:

◆ cascadeMC_GenerateChain()

void Gambit::DarkBit::cascadeMC_GenerateChain ( DarkBit::DecayChain::ChainContainer chain)

Function for generating decay chains.

Option cMC_maxChainLength<int>: Maximum chain length, -1 is infinite (default -1)

Option cMC_Emin<double>: Cutoff energy for cascade particles (default 0)

Definition at line 198 of file Cascades.cpp.

References cascadeMC_DecayTable(), cascadeMC_InitialState(), Gambit::LogTags::err, MC_FINALIZE, MC_INIT, MC_NEXT_STATE, Gambit::piped_errors, and Gambit::Piped_exceptions::request().

Referenced by main().

200  {
201  using namespace DecayChain;
202  using namespace Pipes::cascadeMC_GenerateChain;
203  static int cMC_maxChainLength;
204  static double cMC_Emin;
205  switch(*Loop::iteration)
206  {
207  case MC_INIT:
209  cMC_maxChainLength = runOptions->getValueOrDef<int> (-1, "cMC_maxChainLength");
211  cMC_Emin = runOptions->getValueOrDef<double> (-1, "cMC_Emin");
212  return;
213  case MC_NEXT_STATE:
214  case MC_FINALIZE:
215  return;
216  }
217  shared_ptr<ChainParticle> chn;
218  try
219  {
220  chn.reset(new ChainParticle( vec3(0), &(*Dep::cascadeMC_DecayTable),
222  chn->generateDecayChainMC(cMC_maxChainLength,cMC_Emin);
223  }
224  catch(Piped_exceptions::description err)
225  {
226  Loop::wrapup();
227  piped_errors.request(err);
228  }
229  chain=ChainContainer(chn);
230  }
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:547
void cascadeMC_InitialState(std::string &pID)
Function selecting initial state for decay chain.
Definition: Cascades.cpp:140
void cascadeMC_GenerateChain(DarkBit::DecayChain::ChainContainer &chain)
Function for generating decay chains.
Definition: Cascades.cpp:198
void cascadeMC_DecayTable(DarkBit::DecayChain::DecayTable &table)
Function setting up the decay table used in decay chains.
Definition: Cascades.cpp:58
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cascadeMC_Histograms()

void Gambit::DarkBit::cascadeMC_Histograms ( std::map< std::string, std::map< std::string, SimpleHist > > &  result)

Function responsible for histogramming, and evaluating end conditions for event loop.

Option cMC_numSpecSamples<int>: number of samples to draw from tabulated spectra (default 25)

Option cMC_endCheckFrequency: number of events to wait between successive checks of the convergence criteria (default 25)

Option cMC_gammaBGPower: power-law slope to assume for astrophysical background (default -2.5)

Option cMC_gammaRelError: max allowed relative error in bin with highest expected signal-to-background (default 0.20)

Option cMC_NhistBins<int>: Number of histogram bins (default 140)

Option cMC_binLow<double>: Histogram min energy in GeV (default 0.001)

Option cMC_binHigh<double>: Histogram max energy in GeV (default 10000)

Definition at line 363 of file Cascades.cpp.

References Gambit::DarkBit::SimpleHist::binCenter(), Gambit::DarkBit::SimpleHist::binVals, cascadeMC_FinalStates(), cascadeMC_InitialState(), cascadeMC_sampleSimYield(), DarkBit_warning(), Gambit::FlavBit::LoopFunctions::E(), Gambit::DarkBit::SimpleHist::getBinValues(), Gambit::DarkBit::SimpleHist::getRelError(), Gambit::Piped_exceptions::inquire(), LOCAL_INFO, MC_FINALIZE, MC_INIT, MC_NEXT_STATE, Gambit::DarkBit::SimpleHist::nBins, Gambit::piped_errors, and Gambit::Scanner::pow().

Referenced by cascadeMC_gammaSpectra(), and main().

365  {
366  using namespace DecayChain;
367  using namespace Pipes::cascadeMC_Histograms;
368 
369  // YAML options
370  static int cMC_numSpecSamples;
371  static int cMC_endCheckFrequency;
372  static double cMC_gammaBGPower;
373  static double cMC_gammaRelError;
374  static int cMC_NhistBins;
375  static double cMC_binLow;
376  static double cMC_binHigh;
377  // Histogram list shared between all threads
378  static std::map<std::string, std::map<std::string, SimpleHist> > histList;
379 
380  switch(*Loop::iteration)
381  {
382  case MC_INIT:
383  // Initialization
386  cMC_numSpecSamples = runOptions->getValueOrDef<int> (25, "cMC_numSpecSamples");
389  cMC_endCheckFrequency =
390  runOptions->getValueOrDef<int> (25, "cMC_endCheckFrequency");
393  cMC_gammaBGPower =
394  runOptions->getValueOrDef<double>(-2.5, "cMC_gammaBGPower");
397  cMC_gammaRelError =
398  runOptions->getValueOrDef<double>(0.20, "cMC_gammaRelError");
399 
400  // Note: use same binning for all particle species
402  cMC_NhistBins = runOptions->getValueOrDef<int> (140, "cMC_NhistBins");
404  cMC_binLow = runOptions->getValueOrDef<double>(0.001, "cMC_binLow");
406  cMC_binHigh = runOptions->getValueOrDef<double>(10000.0,"cMC_binHigh");
407  histList.clear();
408  return;
409  case MC_NEXT_STATE:
410  // Initialize histograms
411  for(std::vector<std::string>::const_iterator it =
413  it!=Dep::cascadeMC_FinalStates->end(); ++it)
414  {
415 #ifdef DARKBIT_DEBUG
416  std::cout << "Defining new histList entry!!!" << std::endl;
417  std::cout << "for: " << *Dep::cascadeMC_InitialState
418  << " " << *it << std::endl;
419 #endif
420  histList[*Dep::cascadeMC_InitialState][*it]=
421  SimpleHist(cMC_NhistBins,cMC_binLow,cMC_binHigh,true);
422  }
423  return;
424  case MC_FINALIZE:
425  // For performance, only return the actual result once finished
426  result = histList;
427  return;
428  }
429 
430  // Get list of endpoint states for this chain
431  vector<const ChainParticle*> endpoints;
432  (*Dep::cascadeMC_ChainEvent).chain->
433  collectEndpointStates(endpoints, false);
434  // Iterate over final states of interest
435  for(std::vector<std::string>::const_iterator pit =
437  pit!=Dep::cascadeMC_FinalStates->end(); ++pit)
438  {
439  // Iterate over all endpoint states of the decay chain. These can
440  // either be final state particles themselves or parents of final state
441  // particles. The reason for not using only final state particles is
442  // that certain endpoints (e.g. quark-antiquark pairs) cannot be
443  // handled as separate particles.
444  for(vector<const ChainParticle*>::const_iterator it =endpoints.begin();
445  it != endpoints.end(); it++)
446  {
447 #ifdef DARKBIT_DEBUG
448  std::cout << " working on endpoint: " << (*it)->getpID() << std::endl;
449  (*it)->printChain();
450 #endif
451 
452  // Weighting factor (correction for mismatch between decay width
453  // of available decay channels and total decay width)
454  double weight;
455  // Analyze single particle endpoints
456  bool ignored = true;
457  if((*it)->getnChildren() ==0)
458  {
459  weight = (*it)->getWeight();
460  // Check if the final state itself is the particle we are looking
461  // for.
462  if((*it)->getpID()==*pit)
463  {
464  double E = (*it)->E_Lab();
465 #pragma omp critical (cascadeMC_histList)
466  histList[*Dep::cascadeMC_InitialState][*pit].addEvent(E,weight);
467  ignored = false;
468  }
469  // Check if tabulated spectra exist for this final state
470  else if((*Dep::SimYieldTable).hasChannel( (*it)->getpID(), *pit ))
471  {
473  *Dep::SimYieldTable, *it, *pit, *Dep::TH_ProcessCatalog,
474  histList, *Dep::cascadeMC_InitialState, weight,
475  cMC_numSpecSamples
476  );
477  // Check if an error was raised
478  ignored = false;
479  if(piped_errors.inquire())
480  {
481  Loop::wrapup();
482  return;
483  }
484  }
485  }
486  // Analyze multiparticle endpoints (the endpoint particle is here the
487  // parent of final state particles).
488  else
489  {
490  weight = (*(*it))[0]->getWeight();
491  bool hasTabulated = false;
492  if((*it)->getnChildren() == 2)
493  {
494 #ifdef DARKBIT_DEBUG
495  std::cout << " check whether two-body final state is tabulated: "
496  << (*(*it))[0]->getpID() << " " << (*(*it))[1]->getpID() <<
497  std::endl;
498 #endif
499  // Check if tabulated spectra exist for this final state
500  if((*Dep::SimYieldTable).hasChannel(
501  (*(*it))[0]->getpID() , (*(*it))[1]->getpID(), *pit ))
502  {
503  hasTabulated = true;
504  cascadeMC_sampleSimYield(*Dep::SimYieldTable, *it, *pit,
505  *Dep::TH_ProcessCatalog, histList,
507  cMC_numSpecSamples
508  );
509  // Check if an error was raised
510  ignored = false;
511  if(piped_errors.inquire())
512  {
513  Loop::wrapup();
514  return;
515  }
516  }
517  }
518  if(!hasTabulated)
519  {
520  for(int i=0; i<((*it)->getnChildren()); i++)
521  {
522  const ChainParticle* child = (*(*it))[i];
523  // Check if the child particle is the particle we are looking
524  // for.
525  if(child->getpID()==*pit)
526  {
527  double E = child->E_Lab();
528 #pragma omp critical (cascadeMC_histList)
529  histList[
530  *Dep::cascadeMC_InitialState][*pit].addEvent(E,weight);
531  ignored = false;
532  }
533  // Check if tabulated spectra exist for this final state
534  else if((*Dep::SimYieldTable).hasChannel( child->getpID(),
535  *pit))
536  {
537  cascadeMC_sampleSimYield(*Dep::SimYieldTable, child, *pit,
538  *Dep::TH_ProcessCatalog, histList,
540  cMC_numSpecSamples
541  );
542  // Check if an error was raised
543  ignored = false;
544  if(piped_errors.inquire())
545  {
546  Loop::wrapup();
547  return;
548  }
549  }
550  }
551  }
552  }
553  if (ignored)
554  {
555  DarkBit_warning().raise(LOCAL_INFO,
556  "WARNING FCMC: Missing complete decay information for "
557  + (*it)->getpID() + ". This state is ignored.");
558  }
559  }
560  }
561  // Check if finished every cMC_endCheckFrequency events
562  if((*Loop::iteration % cMC_endCheckFrequency) == 0)
563  {
564  enum status{untouched,unfinished,finished};
565  status cond = untouched;
566  for(std::vector<std::string>::const_iterator it =
568  it != Dep::cascadeMC_FinalStates->end(); ++it)
569  {
570  // End conditions currently only implemented for gamma final state
571  if(*it=="gamma")
572  {
573  SimpleHist hist;
574 #pragma omp critical (cascadeMC_histList)
575  hist = histList[*Dep::cascadeMC_InitialState][*it];
576 #ifdef DARKBIT_DEBUG
577  std::cout << "Checking whether convergence is reached" << std::endl;
578  for ( int i = 0; i < hist.nBins; i++ )
579  std::cout << "Estimated error at " << hist.binCenter(i) << " GeV : " << hist.getRelError(i) << std::endl;
580 #endif
581  double sbRatioMax=-1.0;
582  int maxBin=0;
583  for(int i=0; i<hist.nBins; i++)
584  {
585  double E = hist.binCenter(i);
586  double background = pow(E,cMC_gammaBGPower);
587  double sbRatio = hist.binVals[i]/background;
588  if(sbRatio>sbRatioMax)
589  {
590  sbRatioMax = sbRatio;
591  maxBin=i;
592  }
593  }
594 #ifdef DARKBIT_DEBUG
595  std::cout << "Estimated maxBin: " << maxBin << std::endl;
596  std::cout << "Energy at maxBin: " << hist.binCenter(maxBin) << std::endl;
597  std::cout << "Estimated error at maxBin: " << hist.getRelError(maxBin) << std::endl;
598  std::cout << "Value at maxBin: " << hist.getBinValues()[maxBin];
599 #endif
600  // Check if end condition is fulfilled. If not, set cond to
601  // unfinished.
602  if(hist.getRelError(maxBin) > cMC_gammaRelError) cond = unfinished;
603 
604  // If end condition is fulfilled, set cond to finished, unless
605  // already set to unfinished by another condition.
606  else if(cond != unfinished) cond = finished;
607  }
608  }
609  // Break Monte Carlo loop if all end conditions are fulfilled.
610  if(cond==finished)
611  {
612 #ifdef DARKBIT_DEBUG
613  std::cout << "!! wrapping up !!" << std::endl;
614  std::cout << "Performed iterations: " << *Loop::iteration << std::endl;
615 #endif
616  Loop::wrapup();
617  }
618  }
619  }
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
void cascadeMC_FinalStates(std::vector< std::string > &list)
Function for retrieving list of final states for cascade decays.
Definition: Cascades.cpp:41
#define LOCAL_INFO
Definition: local_info.hpp:34
void cascadeMC_InitialState(std::string &pID)
Function selecting initial state for decay chain.
Definition: Cascades.cpp:140
bool inquire()
Check whether any exceptions were requested without handling them.
Definition: exceptions.cpp:607
void cascadeMC_sampleSimYield(const SimYieldTable &table, const DarkBit::DecayChain::ChainParticle *endpoint, std::string finalState, const TH_ProcessCatalog &catalog, std::map< std::string, std::map< std::string, SimpleHist > > &histList, std::string initialState, double weight, int cMC_numSpecSamples)
Function for sampling SimYieldTables (tabulated spectra).
Definition: Cascades.cpp:235
warning & DarkBit_warning()
double pow(const double &a)
Outputs a^i.
void cascadeMC_Histograms(std::map< std::string, std::map< std::string, SimpleHist > > &result)
Function responsible for histogramming, and evaluating end conditions for event loop.
Definition: Cascades.cpp:363
double E(const double x, const double y)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cascadeMC_InitialState()

void Gambit::DarkBit::cascadeMC_InitialState ( std::string &  pID)

Function selecting initial state for decay chain.

Definition at line 140 of file Cascades.cpp.

References GA_missingFinalStates(), LOCAL_INFO, MC_FINALIZE, MC_INIT, MC_NEXT_STATE, Gambit::piped_errors, and Gambit::Piped_exceptions::request().

Referenced by cascadeMC_EventCount(), cascadeMC_GenerateChain(), cascadeMC_Histograms(), and main().

141  {
142  using namespace DecayChain;
143  using namespace Pipes::cascadeMC_InitialState;
144  std::vector<std::string> chainList = *Dep::GA_missingFinalStates;
145  static int iteration;
146  switch(*Loop::iteration)
147  {
148  case MC_INIT:
149  iteration = -1;
150  return;
151  case MC_NEXT_STATE:
152  iteration++;
153  break;
154  case MC_FINALIZE:
155  return;
156  }
157  if(int(chainList.size()) <= iteration)
158  {
159  Loop::wrapup();
161  "Desync between cascadeMC_LoopManager and cascadeMC_InitialState");
162  }
163  else
164  pID = chainList[iteration];
165 #ifdef DARKBIT_DEBUG
166  std::cout << "cascadeMC_InitialState" << std::endl;
167  std::cout << " Iteration: " << *Loop::iteration << std::endl;
168  std::cout << " Number of states to simulate: "
169  << chainList.size() << std::endl;
170  std::cout << " Current state: " << pID << std::endl;
171 #endif
172  }
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
void GA_missingFinalStates(std::vector< std::string > &result)
Identification of final states that are not yet tabulated.
Definition: GamYields.cpp:61
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:547
#define LOCAL_INFO
Definition: local_info.hpp:34
void cascadeMC_InitialState(std::string &pID)
Function selecting initial state for decay chain.
Definition: Cascades.cpp:140
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cascadeMC_LoopManager()

void Gambit::DarkBit::cascadeMC_LoopManager ( )

Loop manager for cascade decays.

Option cMC_maxEvents<int>: Maximum number of cascade MC runs (default 20000)

Definition at line 80 of file Cascades.cpp.

References Gambit::Piped_invalid_point::check(), Gambit::Piped_exceptions::check(), DarkBit_error(), DarkBit_warning(), GA_missingFinalStates(), Gambit::Piped_exceptions::inquire(), LOCAL_INFO, MC_FINALIZE, MC_INIT, MC_NEXT_STATE, Gambit::piped_errors, Gambit::piped_invalid_point, and Gambit::piped_warnings.

Referenced by dumpSpectrum(), and main().

81  {
82  using namespace Pipes::cascadeMC_LoopManager;
83  std::vector<std::string> chainList = *Dep::GA_missingFinalStates;
84  int cMC_minEvents = 2; // runOptions->getValueOrDef<int>(2, "cMC_minEvents");
85  // Get YAML options
87  int cMC_maxEvents = runOptions->getValueOrDef<int>(20000, "cMC_maxEvents");
88 
89 
90  // Initialization run
91  Loop::executeIteration(MC_INIT);
92 
93  // Check whether there is anything to do
94  if ( chainList.size() == 0 )
95  {
96  return;
97  }
98 
99  // Iterate over initial state particles
100  for(std::vector<std::string>::const_iterator
101  cit =chainList.begin(); cit != chainList.end(); cit++)
102  {
103  int it;
104  int counter = 0;
105  bool finished = false;
106  // Set next initial state
107  Loop::executeIteration(MC_NEXT_STATE);
108  // Event generation loop
109  #pragma omp parallel private(it) shared(counter, finished)
110  {
111  while (!finished)
112  {
113  #pragma omp critical (cascadeMC_Counter)
114  {
115  counter++;
116  it = counter;
117  }
118  Loop::executeIteration(it);
119  #pragma omp critical (cascadeMC_Counter)
120  {
121  if((*Loop::done and ((counter >= cMC_minEvents) or piped_errors.inquire()))
122  or (counter >= cMC_maxEvents))
123  finished=true;
124  if (counter >= cMC_maxEvents)
125  DarkBit_warning().raise(LOCAL_INFO,
126  "WARNING FCMC: cMC_maxEvents reached without convergence.");
127  }
128  }
129  }
130  // Raise any exceptions
134  Loop::reset();
135  }
136  Loop::executeIteration(MC_FINALIZE);
137  }
error & DarkBit_error()
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
void GA_missingFinalStates(std::vector< std::string > &result)
Identification of final states that are not yet tabulated.
Definition: GamYields.cpp:61
#define LOCAL_INFO
Definition: local_info.hpp:34
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
bool inquire()
Check whether any exceptions were requested without handling them.
Definition: exceptions.cpp:607
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:544
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:564
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:481
warning & DarkBit_warning()
void cascadeMC_LoopManager()
Loop manager for cascade decays.
Definition: Cascades.cpp:80
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cascadeMC_sampleSimYield()

void Gambit::DarkBit::cascadeMC_sampleSimYield ( const SimYieldTable table,
const DarkBit::DecayChain::ChainParticle endpoint,
std::string  finalState,
const TH_ProcessCatalog catalog,
std::map< std::string, std::map< std::string, SimpleHist > > &  histList,
std::string  initialState,
double  weight,
int  cMC_numSpecSamples 
)

Function for sampling SimYieldTables (tabulated spectra).

This is a convenience function used in cascadeMC_Histograms, and does not have an associated capability.

Definition at line 235 of file Cascades.cpp.

References Gambit::DarkBit::SimpleHist::addBox(), beta, Gambit::DarkBit::SimYieldChannel::dNdE_bound, Gambit::Random::draw(), Gambit::DarkBit::DecayChain::ChainParticle::E_Lab(), Gambit::DarkBit::DecayChain::ChainParticle::E_parentFrame(), Gambit::DarkBit::SimYieldChannel::Ecm_max, Gambit::DarkBit::SimYieldChannel::Ecm_min, Gambit::DarkBit::DecayChain::ChainParticle::getBoost(), Gambit::DarkBit::SimYieldTable::getChannel(), Gambit::DarkBit::DecayChain::ChainParticle::getnChildren(), Gambit::DarkBit::DecayChain::ChainParticle::getParent(), Gambit::DarkBit::TH_ProcessCatalog::getParticleProperty(), Gambit::DarkBit::DecayChain::ChainParticle::getpID(), LOCAL_INFO, Gambit::Stats::logmin, M, Gambit::DarkBit::DecayChain::ChainParticle::m, Gambit::DarkBit::TH_ParticleProperty::mass, Gambit::DarkBit::SimpleHist::multiply(), Gambit::DarkBit::DecayChain::ChainParticle::p_Lab(), Gambit::piped_errors, and Gambit::Piped_exceptions::request().

Referenced by cascadeMC_Histograms().

243  {
244 #ifdef DARKBIT_DEBUG
245  std::cout << "SampleSimYield" << std::endl;
246 #endif
247  std::string p1,p2;
248  double gamma,beta;
249  double M;
250  switch(endpoint->getnChildren())
251  {
252  case 0:
253  {
254  p1 = endpoint->getpID();
255  p2 = "";
256  const DarkBit::DecayChain::ChainParticle* parent = endpoint->getParent();
257  if(parent == NULL)
258  {
259  endpoint->getBoost(gamma,beta);
260  M = endpoint->m;
261  }
262  else
263  {
264  parent->getBoost(gamma,beta);
265  M = endpoint->E_parentFrame();
266  }
267  break;
268  }
269  case 2:
270  {
271  p1=(*endpoint)[0]->getpID();
272  p2=(*endpoint)[1]->getpID();
273  endpoint->getBoost(gamma,beta);
274  M = endpoint->m;
275  break;
276  }
277  default:
279  "cascadeMC_sampleSimYield called with invalid endpoint state.");
280  return;
281  }
282  const SimYieldChannel &chn = table.getChannel(p1 , p2, finalState);
283  // Get Lorentz boost information
284 
285  const double gammaBeta = gamma*beta;
286  // Mass of final state squared
287  const double m = catalog.getParticleProperty(finalState).mass;
288  const double msq = m*m;
289  // Get histogram edges
290  double histEmin, histEmax;
291  histList[initialState][finalState].getEdges(histEmin, histEmax);
292 
293  // Calculate energies to sample between. A particle decaying
294  // isotropically in its rest frame will give a box spectrum. This is
295  // assumed and used here to add box contributions rather than points to
296  // the histograms. Limits are chosen such that we only sample energies
297  // that can contribute to histogram bins.
298  const double Ecmin = std::max( gamma*histEmin
299  - gammaBeta*sqrt(histEmin*histEmin-msq) , 0*chn.Ecm_min ); // CW: chn.Ecm_min refers to initial not final energies
300  const double Ecmax = std::min(std::min(
301  // Highest energy that can contribute to the histogram
302  gamma*histEmax + gammaBeta*sqrt(histEmax*histEmax-msq),
303  // Highest energy in SimYieldChannel
304  chn.Ecm_max ),
305  // Estimate for highest kinematically allowed CoM energy
306  0.5*(M*M + msq)/M );
307  if(Ecmin>=Ecmax) return;
308  const double logmin = log(Ecmin);
309  const double logmax = log(Ecmax);
310  const double dlogE=logmax-logmin;
311 
312 #ifdef DARKBIT_DEBUG
313  std::cout << "M = " << M << std::endl;
314  std::cout << "E_lab = " << endpoint->E_Lab() << std::endl;
315  std::cout << "p_lab = " << endpoint->p_Lab() << std::endl;
316  std::cout << "Lorentz factors gamma, beta: " << gamma << ", "
317  << beta << std::endl;
318  std::cout << "Initial state: " << initialState << std::endl;
319  std::cout << "Channel: " << p1 << " " << p2 << std::endl;
320  std::cout << "Final particles: " << finalState << std::endl;
321  std::cout << "Event weight: " << weight << std::endl;
322  std::cout << "histEmin/histEmax: " << histEmin << " " << histEmax
323  << std::endl;
324  std::cout << "chn.Ecm_min/max: " << chn.Ecm_min << " " << chn.Ecm_max
325  << std::endl;
326  std::cout << "Ecmin/max: " << Ecmin << " " << Ecmax << std::endl;
327  std::cout << "Final state mass^2: " << msq << std::endl;
328 #endif
329 
330  double specSum=0;
331  int Nsampl=0;
332  SimpleHist spectrum(histList[initialState][finalState].binLower);
333  while(Nsampl<cMC_numSpecSamples)
334  {
335  // Draw an energy in the CoM frame of the endpoint. Logarithmic
336  // sampling.
337  double E_CoM= exp(logmin+(logmax-logmin)*Random::draw());
338  double dN_dE = chn.dNdE_bound->eval(E_CoM, M);
339 
340  double weight2 = E_CoM*dlogE*dN_dE;
341  specSum += weight2;
342  weight2 *= weight;
343  double tmp1 = gamma*E_CoM;
344  double tmp2 = gammaBeta*sqrt(E_CoM*E_CoM-msq);
345  // Add box spectrum to histogram
346  spectrum.addBox(tmp1-tmp2,tmp1+tmp2,weight2);
347  Nsampl++;
348  }
349 #ifdef DARKBIT_DEBUG
350  std::cout << "Number of samples = " << Nsampl << std::endl;
351 #endif
352  if(Nsampl>0)
353  {
354  spectrum.multiply(1.0/Nsampl);
355  // Add bin contents of spectrum histogram to main histogram as weighted
356  // events
357  #pragma omp critical (cascadeMC_histList)
358  histList[initialState][finalState].addHistAsWeights_sameBin(spectrum);
359  }
360  }
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
START_MODEL beta
Definition: Axions.hpp:36
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:547
#define LOCAL_INFO
Definition: local_info.hpp:34
const double logmin
Minimum finite result returnable from log(double x);.
Definition: statistics.cpp:37
START_MODEL M
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CAST_lnL_general()

double Gambit::DarkBit::CAST_lnL_general ( std::vector< double s,
const std::vector< double bkg_counts,
const std::vector< int sig_counts 
)

Definition at line 1350 of file Axions.cpp.

References Gambit::DecayBit::SM_Z::mu.

Referenced by calc_lnL_CAST2007(), and calc_lnL_CAST2017().

1351  {
1352  double result = 0.0;
1353  int n_bins = s.size();
1354 
1355  for (int i = 0; i < n_bins; i++)
1356  {
1357  double mu = s[i] + bkg_counts[i];
1358  result += sig_counts[i]*gsl_sf_log(mu) - mu;
1359  }
1360 
1361  return result;
1362  }
const double mu
Definition: SM_Z.hpp:42
Here is the caller graph for this function:

◆ CMC_dummy()

void Gambit::DarkBit::CMC_dummy ( DarkBit::stringFunkMap result)

Definition at line 56 of file DarkBit_standalone_MSSM.cpp.

57  {
59  result = sfm;
60  }
std::map< str, daFunk::Funk > stringFunkMap
Definition: SimpleHist.hpp:101

◆ createDecays()

void Gambit::DarkBit::createDecays ( DecayTable outDecays)

Option inputFileName<std::string>: Input SLHA (required)

Definition at line 73 of file DarkBit_standalone_MSSM.cpp.

Referenced by main(), and QUICK_FUNCTION().

74  {
75  using namespace Pipes::createDecays;
77  std::string inputFileName = runOptions->getValue<std::string>("filename");
78  std::cout << "Loading decays from: " << inputFileName << std::endl;
79  outDecays = DecayTable(inputFileName, 0, true);
80  }
void createDecays(DecayTable &outDecays)
GAMBIT native decay table class.
Definition: decay_table.hpp:35
Here is the caller graph for this function:

◆ createSLHA1Names()

void Gambit::DarkBit::createSLHA1Names ( mass_es_pseudonyms &  names)

Definition at line 83 of file DarkBit_standalone_MSSM.cpp.

References Gambit::LogTags::debug, and MSSM_spectrum.

Referenced by main().

84  {
85  const double gauge_mixing_tol = 0.5;
86  const bool tol_invalidates_pt = true;
87  const bool debug = false;
88  names.refill(Pipes::createSLHA1Names::Dep::MSSM_spectrum->get_HE(), gauge_mixing_tol, tol_invalidates_pt, debug);
89  }
Here is the caller graph for this function:

◆ createSpectrum()

void Gambit::DarkBit::createSpectrum ( Spectrum outSpec)

Option inputFileName<std::string>: Input SLHA (required)

Definition at line 63 of file DarkBit_standalone_MSSM.cpp.

References createSpectrum.

64  {
65  using namespace Pipes::createSpectrum;
67  std::string inputFileName = runOptions->getValue<std::string>("filename");
68  std::cout << "Loading spectrum from: " << inputFileName << std::endl;
69  outSpec = spectrum_from_SLHA<MSSMSimpleSpec>(inputFileName, Spectrum::mc_info(), Spectrum::mr_info());
70  }
std::vector< YAML::sdd > mc_info
Typedefs for making it easier to manipulate mass cut and mass ratio cut info.
Definition: spectrum.hpp:119
std::vector< YAML::ssdd > mr_info
Definition: spectrum.hpp:120

◆ DarkBit_error()

◆ DarkBit_warning()

◆ DarkMatter_ID_DiracSingletDM()

void Gambit::DarkBit::DarkMatter_ID_DiracSingletDM ( std::string &  result)

Definition at line 202 of file DiracSingletDM.cpp.

202 { result = "F"; }

◆ DarkMatter_ID_MajoranaSingletDM()

void Gambit::DarkBit::DarkMatter_ID_MajoranaSingletDM ( std::string &  result)

Definition at line 205 of file MajoranaSingletDM.cpp.

205 { result = "X"; }

◆ DarkMatter_ID_MSSM()

void Gambit::DarkBit::DarkMatter_ID_MSSM ( std::string &  result)

Definition at line 742 of file MSSM.cpp.

References b, Gambit::Spectrum::get(), Gambit::invalid_point(), MSSM_spectrum, Gambit::Par::Pole_Mass, and Gambit::invalid_point_exception::raise().

Referenced by main().

743  {
744  using namespace Pipes::DarkMatter_ID_MSSM;
745 
747 
748  // Usual candidates are the lightest neutralino, chargino, up-type squark, down-type-squark, slepton, sneutrino or gluino
749  sdpair msqd = {"~d_1", spec.get(Par::Pole_Mass, "~d_1")};
750  sdpair msqu = {"~u_1", spec.get(Par::Pole_Mass, "~u_1")};
751  sdpair msl = {"~e-_1", spec.get(Par::Pole_Mass, "~e-_1")};
752  sdpair msneu = {"~nu_1", spec.get(Par::Pole_Mass, "~nu_1")};
753  sdpair mglui = {"~g", spec.get(Par::Pole_Mass, "~g")};
754  sdpair mchi0 = {"~chi0_1", std::abs(spec.get(Par::Pole_Mass, "~chi0_1"))};
755  sdpair mchip = {"~chi+_1", std::abs(spec.get(Par::Pole_Mass, "~chi+_1"))};
756 
757  auto min = [&](sdpair a, sdpair b) { return a.second < b.second ? a : b; };
758 
759  sdpair lsp = min(min(min(min(msqd, msqu), min(msl, msneu)), min(mchi0, mchip)),mglui);
760 
761  if (lsp == msqd or lsp == msqu or lsp == msl or lsp == mglui or lsp == mchip)
762  {
763  invalid_point().raise("Point invalidated for having charged LSP.");
764  }
765 
766  result = lsp.first;
767  }
std::pair< str, double > sdpair
Shorthand for a pair of string and double.
Definition: util_types.hpp:70
START_MODEL b
Definition: demo.hpp:270
virtual void raise(const std::string &)
Raise the exception, i.e. throw it. Exact override of base method.
Definition: exceptions.cpp:422
void DarkMatter_ID_MSSM(std::string &result)
Definition: MSSM.cpp:742
invalid_point_exception & invalid_point()
Invalid point exceptions.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DarkMatter_ID_ScalarSingletDM()

void Gambit::DarkBit::DarkMatter_ID_ScalarSingletDM ( std::string &  result)

Definition at line 196 of file ScalarSingletDM.cpp.

Referenced by main().

196 { result = "S"; }
Here is the caller graph for this function:

◆ DarkMatter_ID_VectorSingletDM()

void Gambit::DarkBit::DarkMatter_ID_VectorSingletDM ( std::string &  result)

Definition at line 192 of file VectorSingletDM.cpp.

192 { result = "V"; }

◆ DarkMatter_ID_WIMP()

void Gambit::DarkBit::DarkMatter_ID_WIMP ( std::string &  result)

Definition at line 202 of file DarkBit_standalone_WIMP.cpp.

203  {
204  result = "WIMP";
205  }

◆ DarkMatterConj_ID_DiracSingletDM()

void Gambit::DarkBit::DarkMatterConj_ID_DiracSingletDM ( std::string &  result)

Definition at line 203 of file DiracSingletDM.cpp.

203 { result = "F"; }

◆ DarkMatterConj_ID_MajoranaSingletDM()

void Gambit::DarkBit::DarkMatterConj_ID_MajoranaSingletDM ( std::string &  result)

Definition at line 206 of file MajoranaSingletDM.cpp.

206 { result = "X"; }

◆ DarkMatterConj_ID_MSSM()

void Gambit::DarkBit::DarkMatterConj_ID_MSSM ( std::string &  result)

Definition at line 769 of file MSSM.cpp.

References b, Gambit::Spectrum::get(), Gambit::invalid_point(), MSSM_spectrum, Gambit::Par::Pole_Mass, and Gambit::invalid_point_exception::raise().

Referenced by main().

770  {
771  using namespace Pipes::DarkMatterConj_ID_MSSM;
772 
774 
775  // Usual candidates are the lightest neutralino, chargino, up-type squark, down-type-squark, slepton, sneutrino or gluino
776  sdpair msqd = {"~d_1", spec.get(Par::Pole_Mass, "~d_1")};
777  sdpair msqu = {"~u_1", spec.get(Par::Pole_Mass, "~u_1")};
778  sdpair msl = {"~e+_1", spec.get(Par::Pole_Mass, "~e-_1")};
779  sdpair msneu = {"~nu_1", spec.get(Par::Pole_Mass, "~nu_1")};
780  sdpair mglui = {"~g", spec.get(Par::Pole_Mass, "~g")};
781  sdpair mchi0 = {"~chi0_1", std::abs(spec.get(Par::Pole_Mass, "~chi0_1"))};
782  sdpair mchip = {"~chi-_1", std::abs(spec.get(Par::Pole_Mass, "~chi-_1"))};
783 
784  auto min = [&](sdpair a, sdpair b) { return a.second < b.second ? a : b; };
785 
786  sdpair lsp = min(min(min(min(msqd, msqu), min(msl, msneu)), min(mchi0, mchip)),mglui);
787 
788  if (lsp == msqd or lsp == msqu or lsp == msl or lsp == mglui or lsp == mchip)
789  {
790  invalid_point().raise("Point invalidated for having charged LSP.");
791  }
792 
793  result = lsp.first;
794  }
std::pair< str, double > sdpair
Shorthand for a pair of string and double.
Definition: util_types.hpp:70
void DarkMatterConj_ID_MSSM(std::string &result)
Definition: MSSM.cpp:769
START_MODEL b
Definition: demo.hpp:270
virtual void raise(const std::string &)
Raise the exception, i.e. throw it. Exact override of base method.
Definition: exceptions.cpp:422
invalid_point_exception & invalid_point()
Invalid point exceptions.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DarkMatterConj_ID_ScalarSingletDM()

void Gambit::DarkBit::DarkMatterConj_ID_ScalarSingletDM ( std::string &  result)

Definition at line 197 of file ScalarSingletDM.cpp.

Referenced by main().

197 { result = "S"; }
Here is the caller graph for this function:

◆ DarkMatterConj_ID_VectorSingletDM()

void Gambit::DarkBit::DarkMatterConj_ID_VectorSingletDM ( std::string &  result)

Definition at line 193 of file VectorSingletDM.cpp.

193 { result = "V"; }

◆ DarkMatterConj_ID_WIMP()

void Gambit::DarkBit::DarkMatterConj_ID_WIMP ( std::string &  result)

Definition at line 208 of file DarkBit_standalone_WIMP.cpp.

Referenced by dumpSpectrum(), and main().

209  {
210  result = "WIMP";
211  }
Here is the caller graph for this function:

◆ DarkSUSY5_PointInit_LocalHalo_func()

void Gambit::DarkBit::DarkSUSY5_PointInit_LocalHalo_func ( bool &  result)

Function to set Local Halo Parameters in DarkSUSY (DS5 only)

Option v_earth<double>: Keplerian velocity of the Earth around the Sun in km/s (default 29.78)

Definition at line 711 of file SunNeutrinos.cpp.

References Gambit::LogTags::debug, Gambit::EOM, Gambit::logger(), Gambit::LocalMaxwellianHalo::rho0, Gambit::LocalMaxwellianHalo::v0, vesc, Gambit::LocalMaxwellianHalo::vesc, vrot, and Gambit::LocalMaxwellianHalo::vrot.

Referenced by main().

712  {
714 
715  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
716 
717  double rho0 = LocalHaloParameters.rho0;
718  double rho0_eff = (*Dep::RD_fraction)*rho0;
719  double vrot = LocalHaloParameters.vrot;
720  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
721  double vesc = LocalHaloParameters.vesc;
723  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
724 
725  BEreq::dshmcom->rho0 = rho0;
726  BEreq::dshmcom->v_sun = vrot;
727  BEreq::dshmcom->v_earth = v_earth;
728  BEreq::dshmcom->rhox = rho0_eff;
729 
730  BEreq::dshmframevelcom->v_obs = vrot;
731 
732  BEreq::dshmisodf->vd_3d = vd_3d;
733  BEreq::dshmisodf->vgalesc = vesc;
734 
735  BEreq::dshmnoclue->vobs = vrot;
736 
738  << "Updating DarkSUSY halo parameters:" << std::endl
739  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
740  << " rho0_eff [GeV/cm^3] = " << rho0_eff << std::endl
741  << " v_sun [km/s] = " << vrot<< std::endl
742  << " v_earth [km/s] = " << v_earth << std::endl
743  << " v_obs [km/s] = " << vrot << std::endl
744  << " vd_3d [km/s] = " << vd_3d << std::endl
745  << " v_esc [km/s] = " << vesc << EOM;
746 
747  result = true;
748 
749  return;
750  }
void DarkSUSY5_PointInit_LocalHalo_func(bool &result)
Function to set Local Halo Parameters in DarkSUSY (DS5 only)
START_MODEL vesc
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
START_MODEL vrot
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DarkSUSY_PointInit_LocalHalo_func()

void Gambit::DarkBit::DarkSUSY_PointInit_LocalHalo_func ( bool &  result)

Function to set Local Halo Parameters in DarkSUSY (DS 6)

Option v_earth<double>: Keplerian velocity of the Earth around the Sun in km/s (default 29.78)

Definition at line 753 of file SunNeutrinos.cpp.

References Gambit::LogTags::debug, Gambit::EOM, Gambit::logger(), Gambit::LocalMaxwellianHalo::rho0, Gambit::LocalMaxwellianHalo::v0, vesc, Gambit::LocalMaxwellianHalo::vesc, vrot, and Gambit::LocalMaxwellianHalo::vrot.

Referenced by main().

754  {
756 
757  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
758 
759  double rho0 = LocalHaloParameters.rho0;
760  double vrot = LocalHaloParameters.vrot;
761  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
762  double vesc = LocalHaloParameters.vesc;
764  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
765 
766  BEreq::dshmcom->rho0 = rho0;
767  BEreq::dshmcom->v_sun = vrot;
768  BEreq::dshmcom->v_earth = v_earth;
769 
770  BEreq::dshmframevelcom->v_obs = vrot;
771 
772  BEreq::dshmisodf->vd_3d = vd_3d;
773  BEreq::dshmisodf->vgalesc = vesc;
774 
775  BEreq::dshmnoclue->vobs = vrot;
776 
778  << "Updating DarkSUSY halo parameters:" << std::endl
779  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
780  << " v_sun [km/s] = " << vrot<< std::endl
781  << " v_earth [km/s] = " << v_earth << std::endl
782  << " v_obs [km/s] = " << vrot << std::endl
783  << " vd_3d [km/s] = " << vd_3d << std::endl
784  << " v_esc [km/s] = " << vesc << EOM;
785 
786  result = true;
787 
788  return;
789  }
void DarkSUSY_PointInit_LocalHalo_func(bool &result)
Function to set Local Halo Parameters in DarkSUSY (DS 6)
START_MODEL vesc
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
START_MODEL vrot
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DD_couplings_DarkSUSY_DS5()

void Gambit::DarkBit::DD_couplings_DarkSUSY_DS5 ( DM_nucleon_couplings result)

Get direct detection couplings from initialized DarkSUSY 5.

Option rescale_couplings<double>: Rescaling factor for WIMP-nucleon couplings (default 1.)

Definition at line 63 of file DirectDetection.cpp.

References Gambit::LogTags::debug, Gambit::EOM, and Gambit::logger().

Referenced by main().

64  {
65  using namespace Pipes::DD_couplings_DarkSUSY_DS5;
66 
67  double fG;
68 
69  // Set proton hadronic matrix elements
70  BEreq::ddcom->ftp(7) = *Param["fpu"];
71  BEreq::ddcom->ftp(8) = *Param["fpd"];
72  BEreq::ddcom->ftp(10) = *Param["fps"];
73 
74  fG = 2./27.*(1. - *Param["fpu"] - *Param["fpd"] - *Param["fps"]);
75  BEreq::ddcom->ftp(9) = fG;
76  BEreq::ddcom->ftp(11) = fG;
77  BEreq::ddcom->ftp(12) = fG;
78 
79  logger() << LogTags::debug << "DarkSUSY proton hadronic matrix elements set to:" << endl;
80  logger() << LogTags::debug << "ftp(7) = fpu = " << BEreq::ddcom->ftp(7);
81  logger() << LogTags::debug << "\tftp(8) = fpd = " << BEreq::ddcom->ftp(8);
82  logger() << LogTags::debug << "\tftp(10) = fps = " << BEreq::ddcom->ftp(10) << endl;
83  logger() << LogTags::debug << "ftp(9) = ftp(11) = ftp(12) = 2/27 fG = " <<
84  BEreq::ddcom->ftp(9) << EOM;
85 
86  // Set neutron hadronic matrix elements
87  BEreq::ddcom->ftn(7) = *Param["fnu"];
88  BEreq::ddcom->ftn(8) = *Param["fnd"];
89  BEreq::ddcom->ftn(10) = *Param["fns"];
90 
91  fG = 2./27.*(1. - *Param["fnu"] - *Param["fnd"] - *Param["fns"]);
92  BEreq::ddcom->ftn(9) = fG;
93  BEreq::ddcom->ftn(11) = fG;
94  BEreq::ddcom->ftn(12) = fG;
95 
96  logger() << LogTags::debug << "DarkSUSY neutron hadronic matrix elements set to:" << endl;
97  logger() << LogTags::debug << "ftn(7) = fnu = " << BEreq::ddcom->ftn(7);
98  logger() << LogTags::debug << "\tftn(8) = fnd = " << BEreq::ddcom->ftn(8);
99  logger() << LogTags::debug << "\tftn(10) = fns = " << BEreq::ddcom->ftn(10) << endl;
100  logger() << LogTags::debug << "ftn(9) = ftn(11) = ftn(12) = 2/27 fG = " <<
101  BEreq::ddcom->ftn(9) << EOM;
102 
103  // Set deltaq
104  BEreq::ddcom->delu = *Param["deltau"];
105  BEreq::ddcom->deld = *Param["deltad"];
106  BEreq::ddcom->dels = *Param["deltas"];
107  logger() << LogTags::debug << "DarkSUSY delta q set to:" << endl;
108  logger() << LogTags::debug << "delu = delta u = " << BEreq::ddcom->delu;
109  logger() << LogTags::debug << "\tdeld = delta d = " << BEreq::ddcom->deld;
110  logger() << LogTags::debug << "\tdels = delta s = " << BEreq::ddcom->dels << EOM;
111 
112 
113  // Loop corrections and pole removal.
114 
115  // Option loop<bool>: If true, include 1-loop effects discussed in
116  // Drees Nojiri Phys.Rev. D48 (1993) 3483-3501 (default: true)
117  BEreq::ddcom->dddn = runOptions->getValueOrDef<bool>(true,"loop");
118 
119  // Option pole<bool>: If true, include pole in nuclear scattering cross-section.
120  // If false, approximate squark propagator as 1/m_sq^2 (default: false)
121  BEreq::ddcom->ddpole = runOptions->getValueOrDef<bool>(false,"pole");
122 
123  // Some version notes:
124  // The default in DS5 is to for both these options to be false (tree-level cross-section with pole removed).
125  // The default in DS6 is to use Drees-Nojiri (1 loop) with the pole removed, which isn't an
126  // option in the official DS5.1.3. The version in GAMBIT is patched to implement this option though,
127  // so we use it as the default here.
128 
129  // Calling DarkSUSY subroutine dsddgpgn(gps,gns,gpa,gna)
130  // to set all four couplings.
131  std::vector<double> DDcouplings=BEreq::get_DD_couplings();
132  double factor =
134  runOptions->getValueOrDef<double>(1., "rescale_couplings");
135  result.gps = factor*DDcouplings[0];// *= factor;
136  result.gns = factor*DDcouplings[1];// *= factor;
137  result.gpa = factor*DDcouplings[2];// *= factor;
138  result.gna = factor*DDcouplings[3];// *= factor;
139  logger() << LogTags::debug << "DarkSUSY dsddgpgn gives:" << std::endl;
140  logger() << LogTags::debug << " gps = " << result.gps << std::endl;
141  logger() << LogTags::debug << " gns = " << result.gns << std::endl;
142  logger() << LogTags::debug << " gpa = " << result.gpa << std::endl;
143  logger() << LogTags::debug << " gna = " << result.gna << EOM;
144  }
void DD_couplings_DarkSUSY_DS5(DM_nucleon_couplings &result)
Get direct detection couplings from initialized DarkSUSY 5.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DD_couplings_DarkSUSY_MSSM()

void Gambit::DarkBit::DD_couplings_DarkSUSY_MSSM ( DM_nucleon_couplings result)

Get direct detection couplings from DarkSUSY 6 initialized with MSSM module.

Option rescale_couplings<double>: Rescaling factor for WIMP-nucleon couplings (default 1.)

Definition at line 148 of file DirectDetection.cpp.

References Gambit::LogTags::debug, Gambit::EOM, and Gambit::logger().

Referenced by main().

149  {
150  using namespace Pipes::DD_couplings_DarkSUSY_MSSM;
151 
152  double fG;
153 
154  // Set proton hadronic matrix elements
155  BEreq::ddcomlegacy->ftp(7) = *Param["fpu"];
156  BEreq::ddcomlegacy->ftp(8) = *Param["fpd"];
157  BEreq::ddcomlegacy->ftp(10) = *Param["fps"];
158 
159  fG = 2./27.*(1. - *Param["fpu"] - *Param["fpd"] - *Param["fps"]);
160  BEreq::ddcomlegacy->ftp(9) = fG;
161  BEreq::ddcomlegacy->ftp(11) = fG;
162  BEreq::ddcomlegacy->ftp(12) = fG;
163 
164  logger() << LogTags::debug << "DarkSUSY proton hadronic matrix elements set to:" << endl;
165  logger() << LogTags::debug << "ftp(7) = fpu = " << BEreq::ddcomlegacy->ftp(7);
166  logger() << LogTags::debug << "\tftp(8) = fpd = " << BEreq::ddcomlegacy->ftp(8);
167  logger() << LogTags::debug << "\tftp(10) = fps = " << BEreq::ddcomlegacy->ftp(10) << endl;
168  logger() << LogTags::debug << "ftp(9) = ftp(11) = ftp(12) = 2/27 fG = " <<
169  BEreq::ddcomlegacy->ftp(9) << EOM;
170 
171  // Set neutron hadronic matrix elements
172  BEreq::ddcomlegacy->ftn(7) = *Param["fnu"];
173  BEreq::ddcomlegacy->ftn(8) = *Param["fnd"];
174  BEreq::ddcomlegacy->ftn(10) = *Param["fns"];
175 
176  fG = 2./27.*(1. - *Param["fnu"] - *Param["fnd"] - *Param["fns"]);
177  BEreq::ddcomlegacy->ftn(9) = fG;
178  BEreq::ddcomlegacy->ftn(11) = fG;
179  BEreq::ddcomlegacy->ftn(12) = fG;
180 
181  logger() << LogTags::debug << "DarkSUSY neutron hadronic matrix elements set to:" << endl;
182  logger() << LogTags::debug << "ftn(7) = fnu = " << BEreq::ddcomlegacy->ftn(7);
183  logger() << LogTags::debug << "\tftn(8) = fnd = " << BEreq::ddcomlegacy->ftn(8);
184  logger() << LogTags::debug << "\tftn(10) = fns = " << BEreq::ddcomlegacy->ftn(10) << endl;
185  logger() << LogTags::debug << "ftn(9) = ftn(11) = ftn(12) = 2/27 fG = " <<
186  BEreq::ddcomlegacy->ftn(9) << EOM;
187 
188  // Set deltaq
189  BEreq::ddcomlegacy->delu = *Param["deltau"];
190  BEreq::ddcomlegacy->deld = *Param["deltad"];
191  BEreq::ddcomlegacy->dels = *Param["deltas"];
192  logger() << LogTags::debug << "DarkSUSY delta q set to:" << endl;
193  logger() << LogTags::debug << "delu = delta u = " << BEreq::ddcomlegacy->delu;
194  logger() << LogTags::debug << "\tdeld = delta d = " << BEreq::ddcomlegacy->deld;
195  logger() << LogTags::debug << "\tdels = delta s = " << BEreq::ddcomlegacy->dels << EOM;
196 
197 
198  // Loop corrections and pole removal.
199 
200  // Option loop<bool>: If true, include 1-loop effects discussed in
201  // Drees Nojiri Phys.Rev. D48 (1993) 3483-3501 (default: true)
202  BEreq::ddmssmcom->dddn = runOptions->getValueOrDef<bool>(true,"loop");
203 
204  // Option pole<bool>: If true, include pole in nuclear scattering cross-section.
205  // If false, approximate squark propagator as 1/m_sq^2 (default: false)
206  BEreq::ddmssmcom->ddpole = runOptions->getValueOrDef<bool>(false,"pole");
207 
208  // Some version notes:
209  // The default in DS5 is for both these options to be false (tree-level cross-section with pole removed).
210  // The default in DS6 is to use Drees-Nojiri (1 loop) with the pole removed, which isn't an
211  // option in the official DS5.1.3. The version in GAMBIT is patched to implement this option though,
212  // so we use it as the default here.
213 
214  // Calling DarkSUSY subroutine dsddgpgn(gps,gns,gpa,gna)
215  // to set all four couplings.
216  std::vector<double> DDcouplings=BEreq::get_DD_couplings();
217  double factor =
219  runOptions->getValueOrDef<double>(1., "rescale_couplings");
220  result.gps = factor*DDcouplings[0];// *= factor;
221  result.gns = factor*DDcouplings[1];// *= factor;
222  result.gpa = factor*DDcouplings[2];// *= factor;
223  result.gna = factor*DDcouplings[3];// *= factor;
224  logger() << LogTags::debug << "DarkSUSY dsddgpgn gives:" << std::endl;
225  logger() << LogTags::debug << " gps = " << result.gps << std::endl;
226  logger() << LogTags::debug << " gns = " << result.gns << std::endl;
227  logger() << LogTags::debug << " gpa = " << result.gpa << std::endl;
228  logger() << LogTags::debug << " gna = " << result.gna << EOM;
229  }
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
void DD_couplings_DarkSUSY_MSSM(DM_nucleon_couplings &result)
Get direct detection couplings from DarkSUSY 6 initialized with MSSM module.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DD_couplings_DiracSingletDM_Z2()

void Gambit::DarkBit::DD_couplings_DiracSingletDM_Z2 ( DM_nucleon_couplings_fermionic_HP &  result)

Direct detection couplings for the DiracSingletDM_Z2 model.

Definition at line 206 of file DiracSingletDM.cpp.

References Gambit::LogTags::debug, Gambit::Par::dimensionless, Gambit::EOM, Gambit::SubSpectrum::get(), Gambit::Spectrum::get(), Gambit::Spectrum::get_HE(), Gambit::logger(), Gambit::m_neutron, Gambit::m_proton, mh, Gambit::Par::Pole_Mass, and Gambit::Scanner::pow().

207  {
209  const Spectrum& spec = *Dep::DiracSingletDM_Z2_spectrum;
210  const SubSpectrum& he = spec.get_HE();
211  //double mass = spec.get(Par::Pole_Mass,"F");
212  double lambda = he.get(Par::dimensionless,"lF");
213  double cosXI = std::cos(he.get(Par::dimensionless,"xi"));
214  double sinXI = std::sin(he.get(Par::dimensionless,"xi"));
215  double mh = spec.get(Par::Pole_Mass,"h0_1");
216 
217  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
218  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
219  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
220 
221  // SI scalar and pseudoscalar couplings
222  result.gps = lambda*fp*m_proton*cosXI/pow(mh,2);
223  result.gns = lambda*fn*m_neutron*cosXI/pow(mh,2);
224  result.gp_q2 = lambda*fp*m_proton*sinXI/pow(mh,2);
225  result.gn_q2 = lambda*fn*m_neutron*sinXI/pow(mh,2);
226 
227  logger() << LogTags::debug << "Dirac DM DD couplings:" << std::endl;
228  logger() << " gps = " << result.gps << std::endl;
229  logger() << " gns = " << result.gns << std::endl;
230  logger() << " gp_q2 = " << result.gp_q2 << std::endl;
231  logger() << " gn_q2 = " << result.gn_q2 << EOM;
232 
233  } // function DD_couplings_DiracSingletDM_Z2
const double m_proton
const double m_neutron
void DD_couplings_DiracSingletDM_Z2(DM_nucleon_couplings_fermionic_HP &result)
Direct detection couplings for the DiracSingletDM_Z2 model.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
double lambda(double x, double y, double z)
Definition: MSSM_H.hpp:38
double pow(const double &a)
Outputs a^i.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum mh
Here is the call graph for this function:

◆ DD_couplings_MajoranaSingletDM_Z2()

void Gambit::DarkBit::DD_couplings_MajoranaSingletDM_Z2 ( DM_nucleon_couplings_fermionic_HP &  result)

Direct detection couplings for the MajoranaSingletDM_Z2 model.

Definition at line 209 of file MajoranaSingletDM.cpp.

References Gambit::LogTags::debug, Gambit::Par::dimensionless, Gambit::EOM, Gambit::SubSpectrum::get(), Gambit::Spectrum::get(), Gambit::Spectrum::get_HE(), Gambit::logger(), Gambit::m_neutron, Gambit::m_proton, MajoranaSingletDM_Z2_spectrum, mh, Gambit::Par::Pole_Mass, and Gambit::Scanner::pow().

210  {
213  const SubSpectrum& he = spec.get_HE();
214  //double mass = spec.get(Par::Pole_Mass,"X");
215  double lambda = he.get(Par::dimensionless,"lX");
216  double cosXI = std::cos(he.get(Par::dimensionless,"xi"));
217  double sinXI = std::sin(he.get(Par::dimensionless,"xi"));
218  double mh = spec.get(Par::Pole_Mass,"h0_1");
219 
220  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
221  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
222  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
223 
224  // SI scalar and pseudoscalar couplings
225  result.gps = lambda*fp*m_proton*cosXI/pow(mh,2);
226  result.gns = lambda*fn*m_neutron*cosXI/pow(mh,2);
227  result.gp_q2 = lambda*fp*m_proton*sinXI/pow(mh,2);
228  result.gn_q2 = lambda*fn*m_neutron*sinXI/pow(mh,2);
229 
230  logger() << LogTags::debug << "Majorana DM DD couplings:" << std::endl;
231  logger() << " gps = " << result.gps << std::endl;
232  logger() << " gns = " << result.gns << std::endl;
233  logger() << " gp_q2 = " << result.gp_q2 << std::endl;
234  logger() << " gn_q2 = " << result.gn_q2 << EOM;
235 
236  } // function DD_couplings_MajoranaSingletDM_Z2
Spectrum Spectrum Spectrum Spectrum Spectrum MajoranaSingletDM_Z2_spectrum
const double m_proton
const double m_neutron
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
double lambda(double x, double y, double z)
Definition: MSSM_H.hpp:38
double pow(const double &a)
Outputs a^i.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum mh
void DD_couplings_MajoranaSingletDM_Z2(DM_nucleon_couplings_fermionic_HP &result)
Direct detection couplings for the MajoranaSingletDM_Z2 model.
Here is the call graph for this function:

◆ DD_couplings_MicrOmegas()

void Gambit::DarkBit::DD_couplings_MicrOmegas ( DM_nucleon_couplings result)

Get direct detection couplings from initialized MicrOmegas.

Definition at line 234 of file DirectDetection.cpp.

References Gambit::byVal(), DarkBit_error(), Gambit::LogTags::debug, Gambit::EOM, LOCAL_INFO, and Gambit::logger().

Referenced by main().

235  {
236  using namespace Pipes::DD_couplings_MicrOmegas;
237 
238  // Set proton hadronic matrix elements.
239  BEreq::MOcommon->par[2] = *Param["fpd"];
240  BEreq::MOcommon->par[3] = *Param["fpu"];
241  BEreq::MOcommon->par[4] = *Param["fps"];
242 
243  logger() << LogTags::debug << "micrOMEGAs proton hadronic matrix elements set to:" << endl;
244  logger() << LogTags::debug << "ScalarFFPd = fpd = " << BEreq::MOcommon->par[2];
245  logger() << LogTags::debug << "\tScalarFFPu = fpu = " << BEreq::MOcommon->par[3];
246  logger() << LogTags::debug << "\tScalarFFPs = fps = " << BEreq::MOcommon->par[4] << EOM;
247 
248  // Set neutron hadronic matrix elements.
249  BEreq::MOcommon->par[11] = *Param["fnd"];
250  BEreq::MOcommon->par[12] = *Param["fnu"];
251  BEreq::MOcommon->par[13] = *Param["fns"];
252 
253  logger() << LogTags::debug << "micrOMEGAs neutron hadronic matrix elements set to:" << endl;
254  logger() << LogTags::debug << "ScalarFFNd = fnd = " << BEreq::MOcommon->par[11];
255  logger() << LogTags::debug << "\tScalarFFNu = fnu = " << BEreq::MOcommon->par[12];
256  logger() << LogTags::debug << "\tScalarFFNs = fns = " << BEreq::MOcommon->par[13] << EOM;
257 
258  //Set delta q.
259  BEreq::MOcommon->par[5] = *Param["deltad"];
260  BEreq::MOcommon->par[6] = *Param["deltau"];
261  BEreq::MOcommon->par[7] = *Param["deltas"];
262 
263  BEreq::MOcommon->par[14] = *Param["deltau"];
264  BEreq::MOcommon->par[15] = *Param["deltad"];
265  BEreq::MOcommon->par[16] = *Param["deltas"];
266 
267  logger() << LogTags::debug << "micrOMEGAs delta q set to:" << endl;
268  logger() << LogTags::debug << "pVectorFFPd = pVectorFFNu = delta d = "
269  << BEreq::MOcommon->par[5] << endl;
270  logger() << LogTags::debug << "pVectorFFPu = pVectorFFPd = delta u = "
271  << BEreq::MOcommon->par[6] << endl;
272  logger() << LogTags::debug << "pVectorFFPs = pVectorFFNs = delta s = "
273  << BEreq::MOcommon->par[7] << EOM;
274 
275  double p1[2], p2[2], p3[2], p4[2];
276  int error;
277  if (runOptions->getValueOrDef<bool>(true, "box"))
278  error = BEreq::nucleonAmplitudes(byVal(BEreq::FeScLoop.pointer()),
279  byVal(p1), byVal(p2), byVal(p3), byVal(p4));
280  else error = BEreq::nucleonAmplitudes(NULL,
281  byVal(p1), byVal(p2), byVal(p3), byVal(p4));
282  if(error!=0)
283  DarkBit_error().raise(LOCAL_INFO,
284  "micrOMEGAs nucleonAmplitudes function failed with "
285  "error code " + std::to_string(error) + ".");
286 
287  // Rescaling to agree with DarkSUSY convention:
288  result.gps = p1[0]*2;
289  result.gpa = p2[0]*2;
290  result.gns = p3[0]*2;
291  result.gna = p4[0]*2;
292 
293  logger() << LogTags::debug << "micrOMEGAs nucleonAmplitudes gives:" << endl;
294  logger() << LogTags::debug << " gps: " << result.gps << endl;
295  logger() << LogTags::debug << " gns: " << result.gns << endl;
296  logger() << LogTags::debug << " gpa: " << result.gpa << endl;
297  logger() << LogTags::debug << " gna: " << result.gna << EOM;
298  }
error & DarkBit_error()
#define LOCAL_INFO
Definition: local_info.hpp:34
void DD_couplings_MicrOmegas(DM_nucleon_couplings &result)
Get direct detection couplings from initialized MicrOmegas.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
T byVal(T t)
Redirection function to turn an lvalue into an rvalue, so that it is correctly passed by value when d...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DD_couplings_ScalarSingletDM_Z2()

void Gambit::DarkBit::DD_couplings_ScalarSingletDM_Z2 ( DM_nucleon_couplings result)

Direct detection couplings for Z2 scalar singlet DM.

Definition at line 224 of file ScalarSingletDM.cpp.

References get_ScalarSingletDM_DD_couplings(), and ScalarSingletDM_Z2_spectrum.

Referenced by main().

225  {
228  get_ScalarSingletDM_DD_couplings(spec, result, Param);
229  }
void DD_couplings_ScalarSingletDM_Z2(DM_nucleon_couplings &result)
Direct detection couplings for Z2 scalar singlet DM.
void get_ScalarSingletDM_DD_couplings(const Spectrum &spec, DM_nucleon_couplings &result, Models::safe_param_map< safe_ptr< const double > > &Param)
Common code for different scalar singlet direct detection coupling routines.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DD_couplings_ScalarSingletDM_Z3()

void Gambit::DarkBit::DD_couplings_ScalarSingletDM_Z3 ( DM_nucleon_couplings result)

Direct detection couplings for Z3 scalar singlet DM.

Definition at line 232 of file ScalarSingletDM.cpp.

References get_ScalarSingletDM_DD_couplings(), and ScalarSingletDM_Z3_spectrum.

233  {
236  get_ScalarSingletDM_DD_couplings(spec, result, Param);
237  }
void DD_couplings_ScalarSingletDM_Z3(DM_nucleon_couplings &result)
Direct detection couplings for Z3 scalar singlet DM.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum ScalarSingletDM_Z3_spectrum
void get_ScalarSingletDM_DD_couplings(const Spectrum &spec, DM_nucleon_couplings &result, Models::safe_param_map< safe_ptr< const double > > &Param)
Common code for different scalar singlet direct detection coupling routines.
Here is the call graph for this function:

◆ DD_couplings_VectorSingletDM_Z2()

void Gambit::DarkBit::DD_couplings_VectorSingletDM_Z2 ( DM_nucleon_couplings result)

Direct detection couplings for the VectorSingletDM_Z2 model.

Definition at line 196 of file VectorSingletDM.cpp.

References Gambit::LogTags::debug, Gambit::Par::dimensionless, Gambit::EOM, Gambit::SubSpectrum::get(), Gambit::Spectrum::get(), Gambit::Spectrum::get_HE(), Gambit::logger(), Gambit::m_neutron, Gambit::m_proton, mh, Gambit::Par::Pole_Mass, Gambit::Scanner::pow(), and VectorSingletDM_Z2_spectrum.

197  {
200  const SubSpectrum& he = spec.get_HE();
201  double mass = spec.get(Par::Pole_Mass,"V");
202  double lambda = he.get(Par::dimensionless,"lambda_hV");
203  double mh = spec.get(Par::Pole_Mass,"h0_1");
204 
205  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
206  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
207  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
208 
209  result.gps = lambda*fp*m_proton/pow(mh,2)/mass/2;
210  result.gns = lambda*fn*m_neutron/pow(mh,2)/mass/2;
211  result.gpa = 0; // Only SI cross-section
212  result.gna = 0;
213 
214  logger() << LogTags::debug << "Vector DM DD couplings:" << std::endl;
215  logger() << " gps = " << result.gps << std::endl;
216  logger() << " gns = " << result.gns << std::endl;
217  logger() << " gpa = " << result.gpa << std::endl;
218  logger() << " gna = " << result.gna << EOM;
219 
220  } // function DD_couplings_VectorSingletDM_Z2
const double m_proton
const double m_neutron
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
void DD_couplings_VectorSingletDM_Z2(DM_nucleon_couplings &result)
Direct detection couplings for the VectorSingletDM_Z2 model.
double lambda(double x, double y, double z)
Definition: MSSM_H.hpp:38
Spectrum Spectrum Spectrum Spectrum VectorSingletDM_Z2_spectrum
double pow(const double &a)
Outputs a^i.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum mh
Here is the call graph for this function:

◆ DD_couplings_WIMP()

void Gambit::DarkBit::DD_couplings_WIMP ( DM_nucleon_couplings result)

Option gps<double>: gps (default 0)

Option gns<double>: gns (default 0)

Option gpa<double>: gpa (default 0)

Option gna<double>: gna (default 0)

Definition at line 213 of file DarkBit_standalone_WIMP.cpp.

References DD_couplings_WIMP.

214  {
215  using namespace Pipes::DD_couplings_WIMP;
217  result.gps = runOptions->getValueOrDef<double>(0., "gps");
219  result.gns = runOptions->getValueOrDef<double>(0., "gns");
221  result.gpa = runOptions->getValueOrDef<double>(0., "gpa");
223  result.gna = runOptions->getValueOrDef<double>(0., "gna");
224  //std::cout << "DD_coupling says" << std::endl;
225  //std::cout << result.gps << std::endl;
226  }

◆ DEF_FUNKTRAIT()

Gambit::DarkBit::DEF_FUNKTRAIT ( RD_EFF_ANNRATE_FROM_PROCESSCATALOG_TRAIT  ) &

Infer Weff from process catalog.

Definition at line 554 of file RelicDensity.cpp.

References Gambit::DarkBit::TH_Process::channelList, DarkBit_error(), DarkMatter_ID, Gambit::DarkBit::TH_Process::genRateMisc, Gambit::gev2tocm3s1, LOCAL_INFO, daFunk::var(), and daFunk::zero().

556  {
557  using namespace Pipes::RD_eff_annrate_from_ProcessCatalog;
558 
559  std::string DMid= *Dep::DarkMatter_ID;
560  std::string DMbarid = *Dep::DarkMatterConj_ID;
561  TH_Process annProc = (*Dep::TH_ProcessCatalog).getProcess(DMid, DMbarid);
562  double mDM = (*Dep::TH_ProcessCatalog).getParticleProperty(DMid).mass;
563 
564  auto Weff = daFunk::zero("peff");
565  auto peff = daFunk::var("peff");
566  auto s = 4*(peff*peff + mDM*mDM);
567 
568  // Individual contributions to the invariant rate Weff. Note that no
569  // symmetry factor of 1/2 for non-identical initial state particles
570  // (non-self-conjugate DM) should appear here. This factor does explicitly
571  // enter, however, when calculating the relic density in RD_oh2_DS_general.
572  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
573  it != annProc.channelList.end(); ++it)
574  {
575  Weff = Weff +
576  it->genRate->set("v", 2*peff/sqrt(mDM*mDM+peff*peff))*s/gev2tocm3s1;
577  }
578  // Add genRateMisc to Weff
579  Weff = Weff + annProc.genRateMisc->set("v", 2*peff/sqrt(mDM*mDM+peff*peff))*s/gev2tocm3s1;
580  if ( Weff->getNArgs() != 1 )
581  DarkBit_error().raise(LOCAL_INFO,
582  "RD_eff_annrate_from_ProcessCatalog: Wrong number of arguments.\n"
583  "The probable cause are three-body final states, which are not supported for this function."
584  );
585  result = Weff->plain<RD_EFF_ANNRATE_FROM_PROCESSCATALOG_TRAIT>("peff");
586  } // function RD_eff_annrate_from_ProcessCatalog
error & DarkBit_error()
#define LOCAL_INFO
Definition: local_info.hpp:34
Funk var(std::string arg)
Definition: daFunk.hpp:921
Funk zero(Args... argss)
Definition: daFunk.hpp:604
const double gev2tocm3s1
Here is the call graph for this function:

◆ DM_process_from_ProcessCatalog()

void Gambit::DarkBit::DM_process_from_ProcessCatalog ( std::string &  result)

Information about the nature of the DM process in question (i.e.

decay or annihilation) to use the correct scaling for ID in terms of the DM density, phase space, etc.

Definition at line 114 of file DarkBit.cpp.

References DarkMatter_ID, and Gambit::DarkBit::TH_Process::find().

Referenced by dumpSpectrum(), and main().

115  {
117 
118  // Only need to check this once.
119  static bool first = true;
120  if (first)
121  {
122  // Can we find a process that is a decay?
123  // (This logic used so that this capability does not depend on DarkMatterConj_ID)
124  const TH_Process* p = (*Dep::TH_ProcessCatalog).find(*Dep::DarkMatter_ID);
125 
126  // If not, it's an annihilation.
127  if (p == NULL) result = "annihilation";
128  else result = "decay";
129  first = false;
130  }
131 
132  }
void DM_process_from_ProcessCatalog(std::string &result)
Information about the nature of the DM process in question (i.e.
Definition: DarkBit.cpp:114
Here is the call graph for this function:
Here is the caller graph for this function:

◆ dRdE()

double Gambit::DarkBit::dRdE ( double  E,
void params 
)

Definition at line 2227 of file Axions.cpp.

References cspline, and Gambit::Scanner::pow().

Referenced by calc_lnL_XENON1T_DM_Anomaly().

2228  {
2229  struct dRdE_params * par = (struct dRdE_params *)params;
2230  // Efficiency of the Xenon1T experiment from arXiv:2006.09721
2231  // Columns: Energy [keV] | Efficiency [dimensionless]
2232  static AxionInterpolator efficiency (GAMBIT_DIR "/DarkBit/data/XENON1T/efficiency.txt", InterpolationOptions1D::cspline);
2233  return std::exp(-0.5*pow((E - par->m)/par->sigma,2))*efficiency.interpolate(E);
2234  }
double pow(const double &a)
Outputs a^i.
double E(const double x, const double y)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ DSgamma3bdy()

double Gambit::DarkBit::DSgamma3bdy ( double(*)(int &, double &, double &)  IBfunc,
int  IBch,
double  Eg,
double  E1,
double  M_DM,
double  m_1,
double  m_2 
)

Fully initialize DarkSUSY to the current model point.

Only selected MSSM parameter spaces are implemented. Returns bool indicating if point initialization was successful, which is essentially always true for models that satisfy the dependency resolver.

Supported models: MSSM63atQWrapper around DarkSUSY kinematic functions

Definition at line 61 of file MSSM.cpp.

Referenced by TH_ProcessCatalog_DS5_MSSM(), and TH_ProcessCatalog_DS_MSSM().

64  {
65  double E2 = 2*M_DM - Eg - E1;
66  double p12 = E1*E1-m_1*m_1;
67  double p22 = E2*E2-m_2*m_2;
68  double p22min = Eg*Eg+p12-2*Eg*sqrt(p12);
69  double p22max = Eg*Eg+p12+2*Eg*sqrt(p12);
70  // Check if the process is kinematically allowed
71  if((E1 < m_1) || (E2 < m_2) || (p22<p22min) || (p22>p22max))
72  {
73  return 0;
74  }
75  double x = Eg/M_DM; // x = E_gamma/mx
76  double y = (m_2*m_2 + 4*M_DM * (M_DM - E2) ) / (4*M_DM*M_DM);
77  // Here, y = (p+k)^2/(4 mx^2), where p denotes the W+ momentum and k the
78  // photon momentum. (note that the expressions above and below only
79  // apply to the v->0 limit)
80 
81  double result = IBfunc(IBch,x,y);
82 
83  #ifdef DARKBIT_DEBUG
84  std::cout << " x, y = " << x << ", " << y << std::endl;
85  std::cout << " E, E1, E2 = " << Eg << ", " << E1 << ", "
86  << E2 << std::endl;
87  std::cout << " mDM, m1, m2 = " << M_DM << ", " << m_1 << ", "
88  << m_2 << std::endl;
89  std::cout << " IBfunc = " << result << std::endl;
90  #endif
91 
92  // M_DM^-2 is from the Jacobi determinant
93  return std::max(0., result) / (M_DM*M_DM);
94  }
Here is the caller graph for this function:

◆ dump_GammaSpectrum()

void Gambit::DarkBit::dump_GammaSpectrum ( double result)

Helper function to dump gamma-ray spectra.

Definition at line 475 of file SimpleLikelihoods.cpp.

References Gambit::EOM, Gambit::logger(), and Gambit::Scanner::pow().

Referenced by dumpSpectrum(), and main().

476  {
477  using namespace Pipes::dump_GammaSpectrum;
478  daFunk::Funk spectrum = (*Dep::GA_Yield)->set("v", 0.);
479  // Option filename<string>: Filename for gamma-ray spectrum dump
480  // (default: dNdE.dat)
481  std::string filename = runOptions->getValueOrDef<std::string>(
482  "dNdE.dat", "filename");
483  logger() << "FILENAME for gamma dump: " << filename << EOM;
484  std::ofstream myfile (filename);
485  if (myfile.is_open())
486  {
487  for (int i = 0; i<=1200; i++)
488  {
489  double energy = pow(10., i/200. - 4.);
490 
491  myfile << energy << " " << spectrum->bind("E")->eval(energy) << "\n";
492  }
493  myfile.close();
494  }
495  result = 0.;
496  }
void dump_GammaSpectrum(double &result)
Helper function to dump gamma-ray spectra.
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
double pow(const double &a)
Outputs a^i.
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
Here is the call graph for this function:
Here is the caller graph for this function:

◆ equation_Tosc()

double Gambit::DarkBit::equation_Tosc ( double  T,
void params 
)

Definition at line 1661 of file Axions.cpp.

References axion_mass_temp(), beta, gStar(), Gambit::m_planck_red, Gambit::pi, and Tchi.

Referenced by calc_AxionOscillationTemperature().

1662  {
1663  // T/MeV, ma0/eV, m_pl/10^12eV = m_pl/10^3 GeV
1664  const double m_pl = m_planck_red*1.0E-3;
1665  struct AxionEDT_params * p1 = (struct AxionEDT_params *)params;
1666  double ma0 = (p1->ma0);
1667  double beta = (p1->beta);
1668  double Tchi = (p1->Tchi);
1669 
1670  double result = 1.0 - gStar(T)*gsl_pow_2(T*T*pi/(ma0*m_pl*axion_mass_temp(T, beta, Tchi)))/10.0;
1671 
1672  return result;
1673  }
START_MODEL beta
Definition: Axions.hpp:36
double gStar(double T)
Various capabilities and functions to provide SM physics as well as QCD input for axions...
Definition: Axions.cpp:1075
const double m_planck_red
const double pi
START_MODEL Tchi
Definition: Axions.hpp:36
double axion_mass_temp(double T, double beta, double Tchi)
Definition: Axions.cpp:1648
Here is the call graph for this function:
Here is the caller graph for this function:

◆ equilibration_time_Sun()

void Gambit::DarkBit::equilibration_time_Sun ( double result)

Equilibration time for capture and annihilation of dark matter in the Sun (s)

Definition at line 176 of file SunNeutrinos.cpp.

References Gambit::DarkBit::TH_Process::channelList, DarkBit_error(), DarkMatter_ID, Gambit::DarkBit::TH_Process::find(), Gambit::DarkBit::TH_Process::genRateMisc, LOCAL_INFO, and Gambit::Scanner::pow().

Referenced by annihilation_rate_Sun(), and main().

177  {
178  using namespace Pipes::equilibration_time_Sun;
179 
180  double sigmav = 0;
181  double T_Sun_core = 1.35e-6; // Sun's core temperature (GeV)
182 
183  std::string DMid = *Dep::DarkMatter_ID;
184  std::string DMbarid = *Dep::DarkMatterConj_ID;
185 
186  // Make sure that we're not trying to work with decaying DM.
187  const TH_Process* p = Dep::TH_ProcessCatalog->find(DMid, DMbarid);
188  if (p == NULL) DarkBit_error().raise(LOCAL_INFO, "Sorry, decaying DM is not supported yet by the DarkBit neutrino routines.");
189  TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMbarid);
190 
191  // Add all the regular channels
192  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
193  it != annProc.channelList.end(); ++it)
194  {
195  if ( it->nFinalStates == 2 )
196  {
197  // (sv)(v = sqrt(2T/mDM)) for two-body final state
198  sigmav += it->genRate->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
199  }
200  }
201 
202  // Add invisible contributions
203  sigmav += annProc.genRateMisc->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
204 
205  double ca = sigmav/6.6e28 * pow(*Dep::mwimp/20.0, 1.5);
206  // Scale the annihilation rate down by a factor of two if the DM is not self-conjugate
207  if (not (*Dep::TH_ProcessCatalog).getProcess(*Dep::DarkMatter_ID, *Dep::DarkMatterConj_ID).isSelfConj) ca *= 0.5;
208  result = pow(*Dep::capture_rate_Sun * ca, -0.5);
209 
210  // std::cout << "v = " << sqrt(2.0*T_Sun_core/(*Dep::mwimp)) << " and sigmav inside equilibration_time_Sun = " << sigmav << std::endl;
211  // std::cout << "capture_rate_Sun inside equilibration_time_Sun = " << *Dep::capture_rate_Sun << std::endl;
212  }
error & DarkBit_error()
void equilibration_time_Sun(double &result)
Equilibration time for capture and annihilation of dark matter in the Sun (s)
#define LOCAL_INFO
Definition: local_info.hpp:34
double pow(const double &a)
Outputs a^i.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ erg_integrand()

double Gambit::DarkBit::erg_integrand ( double  erg,
void params 
)

Definition at line 710 of file Axions.cpp.

References Gambit::gev2cm, Gambit::pi, rad_integrand(), rel_prec, and rs.

Referenced by alt_erg_integrand().

711  {
712  const double eVm = gev2cm*1E7;
713  const double L = 9.26/eVm;
714  struct SolarModel_params3 * p3 = (struct SolarModel_params3 *)params;
715  SolarModel* sol = p3->sol;
716  double m_ax = p3->ma0;
717  double rs = p3->rs;
718 
719  double argument = 0.25*1.0E-3*L*m_ax*m_ax/erg;
720  double temp = gsl_pow_2(gsl_sf_sinc(argument/pi));
721  double exposure = p3->eff_exp->interpolate(erg);
722  //std::cout << "Energy: " << erg << ", expoure = " << exposure << "." << std::endl;
723  SolarModel_params2 p2 = {erg, rs, sol};
724 
725  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1E6);
726  double result, error;
727 
728  gsl_function F;
729  F.function = &rad_integrand;
730  F.params = &p2;
731 
732  // Max. and min. integration radius
733  double rmin = sol->r_lo, rmax = std::min(rs, sol->r_hi);
734 
735  gsl_integration_qag (&F, rmin, rmax, 1e-1*abs_prec, 1e-1*rel_prec, 1E6, method, w, &result, &error);
736  gsl_integration_workspace_free (w);
737 
738  return temp*exposure*result;
739  }
START_MODEL rs
const double pi
const int method
Definition: Axions.cpp:655
const double gev2cm
double rad_integrand(double rad, void *params)
Definition: Axions.cpp:685
const double rel_prec
Definition: Axions.cpp:654
const double abs_prec
Definition: Axions.cpp:654
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ExtractLocalMaxwellianHalo()

void Gambit::DarkBit::ExtractLocalMaxwellianHalo ( LocalMaxwellianHalo result)

Module function providing local density and velocity dispersion parameters.

Definition at line 177 of file DarkBit.cpp.

References Gambit::LocalMaxwellianHalo::rho0, v0, Gambit::LocalMaxwellianHalo::v0, vesc, Gambit::LocalMaxwellianHalo::vesc, vrot, and Gambit::LocalMaxwellianHalo::vrot.

Referenced by main().

178  {
179  using namespace Pipes::ExtractLocalMaxwellianHalo;
180  double rho0 = *Param["rho0"];
181  double v0 = *Param["v0"];
182  double vesc = *Param["vesc"];
183  double vrot = *Param["vrot"];
184  result.rho0 = rho0;
185  result.v0 = v0;
186  result.vesc = vesc;
187  result.vrot = vrot;
188  }
START_MODEL v0
START_MODEL vesc
START_MODEL vrot
void ExtractLocalMaxwellianHalo(LocalMaxwellianHalo &result)
Module function providing local density and velocity dispersion parameters.
Definition: DarkBit.cpp:177
Here is the caller graph for this function:

◆ GA_AnnYield_General()

void Gambit::DarkBit::GA_AnnYield_General ( daFunk::Funk result)

General routine to derive annihilation yield.

Depends on:

This function returns

k*dN/dE*(sv)/mDM**2 (E, v) [cm^3/s/GeV^3]

the energy spectrum of photons times sigma*v/m^2, as function of energy (GeV) and velocity (c), multiplied by k=1 for self-conjugate DM or k=1/2 for non-self conjugate. By default, only the v=0 component is calculated.

Option line_width<double>: Set relative line width used in gamma-ray spectra (default 0.03)

Definition at line 420 of file GamYields.cpp.

References cascadeMC_gammaSpectra, DarkBit_error(), DarkMatter_ID, Gambit::DarkBit::TH_Process::find(), getYield(), daFunk::ifelse(), Gambit::DarkBit::TH_Process::isSelfConj, LOCAL_INFO, daFunk::throwError(), and daFunk::var().

Referenced by dumpSpectrum(), and main().

421  {
422  using namespace Pipes::GA_AnnYield_General;
423 
424  std::string DMid= *Dep::DarkMatter_ID;
425  std::string DMbarid = *Dep::DarkMatterConj_ID;
426 
428  double line_width = runOptions->getValueOrDef<double>(0.03, "line_width");
429 
430  // Check that the ProcessCatalog process *is* an annihilation, not a decay.
431  const TH_Process* p = (*Dep::TH_ProcessCatalog).find(DMid, DMbarid);
432  if (p == NULL)
433  DarkBit_error().raise(LOCAL_INFO, "Annihilation process not found. For decaying DM please use GA_DecayYield_General.");
434 
435  // if (not *Dep::TH_ProcessCatalog.isAnnihilation)
436  // DarkBit_error().raise(LOCAL_INFO, "Process is not an annihilation. Please use GA_DecayYield_General.");
437 
438  // Get annihilation process from process catalog
439  TH_Process annProc = (*Dep::TH_ProcessCatalog).getProcess(DMid, DMbarid);
440 
441  // If process involves non-self-conjugate DM then we need to add a factor of 1/2
442  // to the final spectrum. This must be explicitly set in the process catalog.
443  double k = (annProc.isSelfConj) ? 1. : 0.5;
444 
445  // Get particle mass from process catalog
446  const double mass = (*Dep::TH_ProcessCatalog).getParticleProperty(DMid).mass;
447  const double Ecm = 2*mass;
448 
449  // Loop over all channels for that process
450  daFunk::Funk Yield = getYield(annProc, Ecm, mass, *Dep::TH_ProcessCatalog, *Dep::SimYieldTable,
451  line_width, *Dep::cascadeMC_gammaSpectra);
452 
453  result = k*daFunk::ifelse(1e-6 - daFunk::var("v"), Yield/(mass*mass),
454  daFunk::throwError("Spectrum currently only defined for v=0."));
455  }
error & DarkBit_error()
Funk throwError(std::string msg)
Definition: daFunk.hpp:1411
#define LOCAL_INFO
Definition: local_info.hpp:34
Funk var(std::string arg)
Definition: daFunk.hpp:921
void GA_AnnYield_General(daFunk::Funk &result)
General routine to derive annihilation yield.
Definition: GamYields.cpp:420
daFunk::Funk getYield(TH_Process process, double Ecm, double mass, TH_ProcessCatalog catalog, SimYieldTable table, double line_width, stringFunkMap cascadeMC_gammaSpectra)
Helper function returning yield from a given DM process.
Definition: GamYields.cpp:182
Funk ifelse(Funk f, Funk g, Funk h)
Definition: daFunk.hpp:1373
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GA_DecayYield_General()

void Gambit::DarkBit::GA_DecayYield_General ( daFunk::Funk result)

General routine to derive yield from decaying DM.

Depends on:

This function returns

dN/dE*(Gamma)/mDM (E, v) [cm^3/s/GeV^3]

the energy spectrum of photons times Gamma/m, as function of energy (GeV) and velocity (c). By default, only the v=0 component is calculated.

Option line_width<double>: Set relative line width used in gamma-ray spectra (default 0.03)

Definition at line 473 of file GamYields.cpp.

References cascadeMC_gammaSpectra, DarkBit_error(), DarkMatter_ID, Gambit::DarkBit::TH_Process::find(), getYield(), and LOCAL_INFO.

474  {
475  using namespace Pipes::GA_DecayYield_General;
476 
477  std::string DMid= *Dep::DarkMatter_ID;
478 
480  double line_width = runOptions->getValueOrDef<double>(0.03, "line_width");
481 
482  // // Check that the ProcessCatalog process *is* a decay, not an annihilation.
483  const TH_Process* p = (*Dep::TH_ProcessCatalog).find(DMid);
484  if (p == NULL)
485  DarkBit_error().raise(LOCAL_INFO, "Decay process not found. For annihilating DM please use GA_AnnYield_General.");
486 
487  // // Check that the ProcessCatalog process *is* a decay, not an annihilation.
488  // if (*Dep::TH_ProcessCatalog.isAnnihilation)
489  // DarkBit_error().raise(LOCAL_INFO, "Process is not a decay. Please use GA_AnnYield_General.");
490 
491  // Get annihilation process from process catalog
492  TH_Process decayProc = (*Dep::TH_ProcessCatalog).getProcess(DMid);
493 
494  // Get particle mass from process catalog
495  const double mass = (*Dep::TH_ProcessCatalog).getParticleProperty(DMid).mass;
496  double Ecm = mass;
497 
498  // Loop over all channels for that process
499  daFunk::Funk Yield = getYield(decayProc, Ecm, mass, *Dep::TH_ProcessCatalog, *Dep::SimYieldTable,
500  line_width, *Dep::cascadeMC_gammaSpectra);
501 
502  // Rescale the yield by the correct kinematic factor
503  result = Yield/(mass);
504  }
error & DarkBit_error()
#define LOCAL_INFO
Definition: local_info.hpp:34
void GA_DecayYield_General(daFunk::Funk &result)
General routine to derive yield from decaying DM.
Definition: GamYields.cpp:473
daFunk::Funk getYield(TH_Process process, double Ecm, double mass, TH_ProcessCatalog catalog, SimYieldTable table, double line_width, stringFunkMap cascadeMC_gammaSpectra)
Helper function returning yield from a given DM process.
Definition: GamYields.cpp:182
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
Here is the call graph for this function:

◆ GA_missingFinalStates()

void Gambit::DarkBit::GA_missingFinalStates ( std::vector< std::string > &  result)

Identification of final states that are not yet tabulated.

Structure

1) Go through process catalog and find all final states that require to be calculated in the cascade code. To this end, check whether two-body channels are tabulated for two-body final states, and whether one-particle spectra exist for one-particle final states.

2) Calculate via the cascade code the missing energy spectra.

3) Put together the full spectrum.

Option ignore_all<bool>: Ignore all missing final states (default false)

Option ignore_two_body<bool>: Ignore two-body missing final states (default false)

Option ignore_three_body<bool>: Ignore three-body missing final states (default false)

Definition at line 61 of file GamYields.cpp.

References Gambit::DarkBit::TH_Process::channelList, and DarkMatter_ID.

Referenced by cascadeMC_gammaSpectra(), cascadeMC_InitialState(), cascadeMC_LoopManager(), dumpSpectrum(), and main().

62  {
63  using namespace Pipes::GA_missingFinalStates;
64  std::set<std::string> missingFinalStates;
65  std::string DMid= *Dep::DarkMatter_ID;
66  std::string DMbarid = *Dep::DarkMatterConj_ID;
67 
69  if ( runOptions->getValueOrDef(false, "ignore_all") ) return;
70 
71  // What type of process are we dealing with?
72  std::string processtype = *Dep::DM_process;
73  TH_Process process = (*Dep::DM_process == "annihilation") ?
74  (*Dep::TH_ProcessCatalog).getProcess(DMid, DMbarid) : (*Dep::TH_ProcessCatalog).getProcess(DMid);
75 
76  // Add only gamma-ray spectra for two and three body final states
77  for (std::vector<TH_Channel>::iterator it = process.channelList.begin();
78  it != process.channelList.end(); ++it)
79  {
80  if ( it->nFinalStates == 2 )
81  {
83  if ( not runOptions->getValueOrDef(false, "ignore_two_body") )
84  {
85  #ifdef DARKBIT_DEBUG
86  std::cout << "Checking for missing two-body final states: "
87  << it->finalStateIDs[0] << " " << it->finalStateIDs[1] << std::endl;
88  #endif
89  if ( not Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], it->finalStateIDs[1], "gamma") )
90  {
91  if ( not Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], "gamma") )
92  missingFinalStates.insert(it->finalStateIDs[0]);
93  if ( not Dep::SimYieldTable->hasChannel(it->finalStateIDs[1], "gamma") )
94  missingFinalStates.insert(it->finalStateIDs[1]);
95  }
96  }
97  }
98  else if ( it->nFinalStates == 3 )
99  {
101  if ( not runOptions->getValueOrDef(false, "ignore_three_body") )
102  {
103  #ifdef DARKBIT_DEBUG
104  std::cout << "Checking for missing three-body final states: "
105  << it->finalStateIDs[0] << " " << it->finalStateIDs[1]
106  << " " << it->finalStateIDs[2] << std::endl;
107  #endif
108  if (not Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], "gamma") )
109  missingFinalStates.insert(it->finalStateIDs[0]);
110  if (not Dep::SimYieldTable->hasChannel(it->finalStateIDs[1], "gamma") )
111  missingFinalStates.insert(it->finalStateIDs[1]);
112  if (not Dep::SimYieldTable->hasChannel(it->finalStateIDs[2], "gamma") )
113  missingFinalStates.insert(it->finalStateIDs[2]);
114  }
115  }
116  }
117  // Remove particles we don't have decays for.
118  for (auto it = missingFinalStates.begin(); it != missingFinalStates.end();)
119  {
120  if ((*Dep::TH_ProcessCatalog).find(*it, "") == NULL)
121  {
122  #ifdef DARKBIT_DEBUG
123  std::cout << "Erasing (because no decays known): " << *it << std::endl;
124  #endif
125  missingFinalStates.erase(it++);
126  }
127  else
128  {
129  #ifdef DARKBIT_DEBUG
130  std::cout << "Keeping (because decay known): " << *it << std::endl;
131  #endif
132  ++it;
133  }
134  }
135 
136  #ifdef DARKBIT_DEBUG
137  std::cout << "Number of missing final states: " << missingFinalStates.size() << std::endl;
138  for (auto it = missingFinalStates.begin(); it != missingFinalStates.end(); it++)
139  {
140  std::cout << *it << std::endl;
141  }
142  #endif
143 
144  result.assign(missingFinalStates.begin(), missingFinalStates.end());
145  }
void GA_missingFinalStates(std::vector< std::string > &result)
Identification of final states that are not yet tabulated.
Definition: GamYields.cpp:61
Here is the caller graph for this function:

◆ GalacticHalo_Einasto()

void Gambit::DarkBit::GalacticHalo_Einasto ( GalacticHaloProperties result)

Module function to generate GalacticHaloProperties for Einasto profile.

Definition at line 165 of file DarkBit.cpp.

References alpha, Gambit::GalacticHaloProperties::DensityProfile, daFunk::func(), profile_Einasto(), r_sun, Gambit::GalacticHaloProperties::r_sun, rs, and daFunk::var().

Referenced by main().

166  {
167  using namespace Pipes::GalacticHalo_Einasto;
168  double rhos = *Param["rhos"];
169  double rs = *Param["rs"];
170  double r_sun = *Param["r_sun"];
171  double alpha = *Param["alpha"];
172  result.DensityProfile = daFunk::func(profile_Einasto, rhos, rs, alpha, daFunk::var("r"));
173  result.r_sun = r_sun;
174  }
double profile_Einasto(double rhos, double rs, double alpha, double r)
Einasto dark matter halo profile function.
Definition: DarkBit.cpp:147
START_MODEL r_sun
START_MODEL rs
Funk var(std::string arg)
Definition: daFunk.hpp:921
START_MODEL alpha
Funk func(double(*f)(funcargs...), Args... args)
Definition: daFunk.hpp:768
void GalacticHalo_Einasto(GalacticHaloProperties &result)
Module function to generate GalacticHaloProperties for Einasto profile.
Definition: DarkBit.cpp:165
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GalacticHalo_gNFW()

void Gambit::DarkBit::GalacticHalo_gNFW ( GalacticHaloProperties result)

Module function to generate GalacticHaloProperties for gNFW profile.

Definition at line 151 of file DarkBit.cpp.

References alpha, beta, Gambit::GalacticHaloProperties::DensityProfile, daFunk::func(), profile_gNFW(), r_sun, Gambit::GalacticHaloProperties::r_sun, rs, and daFunk::var().

152  {
153  using namespace Pipes::GalacticHalo_gNFW;
154  double rhos = *Param["rhos"];
155  double rs = *Param["rs"];
156  double r_sun = *Param["r_sun"];
157  double alpha = *Param["alpha"];
158  double beta = *Param["beta"];
159  double gamma = *Param["gamma"];
160  result.DensityProfile = daFunk::func(profile_gNFW, rhos, rs, alpha, beta, gamma, daFunk::var("r"));
161  result.r_sun = r_sun;
162  }
START_MODEL beta
Definition: Axions.hpp:36
START_MODEL r_sun
START_MODEL rs
Funk var(std::string arg)
Definition: daFunk.hpp:921
void GalacticHalo_gNFW(GalacticHaloProperties &result)
Module function to generate GalacticHaloProperties for gNFW profile.
Definition: DarkBit.cpp:151
START_MODEL alpha
Funk func(double(*f)(funcargs...), Args... args)
Definition: daFunk.hpp:768
double profile_gNFW(double rhos, double rs, double alpha, double beta, double gamma, double r)
Generalized NFW dark matter halo profile function.
Definition: DarkBit.cpp:143
Here is the call graph for this function:

◆ gaussian_nuisance_lnL()

double Gambit::DarkBit::gaussian_nuisance_lnL ( double  theo,
double  obs,
double  sigma 
)

Definition at line 1057 of file Axions.cpp.

References Gambit::Stats::gaussian_loglikelihood().

Referenced by QCDAxion_AxionPhotonConstant_Nuisance_lnL(), QCDAxion_TemperatureDependence_Nuisance_lnL(), and QCDAxion_ZeroTemperatureMass_Nuisance_lnL().

1057 { return Stats::gaussian_loglikelihood(theo, obs, 0, sigma, false); }
double gaussian_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is Gaussian...
Definition: statistics.cpp:133
const double sigma
Definition: SM_Z.hpp:43
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_ScalarSingletDM_DD_couplings()

void Gambit::DarkBit::get_ScalarSingletDM_DD_couplings ( const Spectrum spec,
DM_nucleon_couplings result,
Models::safe_param_map< safe_ptr< const double > > &  Param 
)

Common code for different scalar singlet direct detection coupling routines.

Definition at line 200 of file ScalarSingletDM.cpp.

References Gambit::LogTags::debug, Gambit::Par::dimensionless, Gambit::EOM, Gambit::SubSpectrum::get(), Gambit::Spectrum::get(), Gambit::Spectrum::get_HE(), Gambit::logger(), Gambit::m_neutron, Gambit::m_proton, mh, Gambit::Par::Pole_Mass, and Gambit::Scanner::pow().

Referenced by DD_couplings_ScalarSingletDM_Z2(), and DD_couplings_ScalarSingletDM_Z3().

201  {
202  const SubSpectrum& he = spec.get_HE();
203  double mass = spec.get(Par::Pole_Mass,"S");
204  double lambda = he.get(Par::dimensionless,"lambda_hS");
205  double mh = spec.get(Par::Pole_Mass,"h0_1");
206 
207  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
208  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
209  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
210 
211  result.gps = lambda*fp*m_proton/pow(mh,2)/mass/2;
212  result.gns = lambda*fn*m_neutron/pow(mh,2)/mass/2;
213  result.gpa = 0; // Only SI cross-section
214  result.gna = 0;
215 
216  logger() << LogTags::debug << "Singlet DM DD couplings:" << std::endl;
217  logger() << " gps = " << result.gps << std::endl;
218  logger() << " gns = " << result.gns << std::endl;
219  logger() << " gpa = " << result.gpa << std::endl;
220  logger() << " gna = " << result.gna << EOM;
221  }
const double m_proton
const double m_neutron
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:100
EXPORT_SYMBOLS Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
double lambda(double x, double y, double z)
Definition: MSSM_H.hpp:38
double pow(const double &a)
Outputs a^i.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum mh
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_semi_ann_MicrOmegas()

void Gambit::DarkBit::get_semi_ann_MicrOmegas ( double result)

Definition at line 1071 of file RelicDensity.cpp.

References Gambit::byVal().

1072  {
1073  using namespace Pipes::get_semi_ann_MicrOmegas;
1074 
1075  double Beps; // Beps=1e-5 recommended, Beps=1 switches coannihilation off
1076  Beps = runOptions->getValueOrDef<double>(1e-5, "Beps");
1077 
1078  double Xf = *Dep::Xf;
1079 
1080  char*n1 = (char *)"~SS";
1081  char*n2 = (char *)"~SS";
1082  char*n3 = (char *)"h";
1083  char*n4 = (char *)"~ss";
1084 
1085  result = BEreq::get_oneChannel(byVal(Xf),byVal(Beps),byVal(n1),byVal(n2),byVal(n3),byVal(n4));
1086 
1087  }
void get_semi_ann_MicrOmegas(double &result)
T byVal(T t)
Redirection function to turn an lvalue into an rvalue, so that it is correctly passed by value when d...
Here is the call graph for this function:

◆ getYield()

daFunk::Funk Gambit::DarkBit::getYield ( TH_Process  process,
double  Ecm,
double  mass,
TH_ProcessCatalog  catalog,
SimYieldTable  table,
double  line_width,
stringFunkMap  cascadeMC_gammaSpectra 
)

Helper function returning yield from a given DM process.

Definition at line 182 of file GamYields.cpp.

References boost_dNdE(), Gambit::DarkBit::TH_Process::channelList, DarkBit_warning(), Gambit::LogTags::debug, Gambit::FlavBit::LoopFunctions::E(), daFunk::func(), Gambit::DarkBit::DarkBit_utils::gamma3bdy_limits(), Gambit::DarkBit::DarkBit_utils::gamma3bdy_limits< 0 >(), Gambit::DarkBit::DarkBit_utils::gamma3bdy_limits< 1 >(), Gambit::DarkBit::TH_ProcessCatalog::getParticleProperty(), Gambit::DarkBit::SimYieldTable::hasChannel(), Gambit::DarkBit::TH_Process::isAnnihilation, LOCAL_INFO, daFunk::logspace(), generate_raster_scan_settings::m1, generate_raster_scan_settings::m2, Gambit::DarkBit::TH_ParticleProperty::mass, Gambit::Scanner::pow(), daFunk::var(), and daFunk::zero().

Referenced by GA_AnnYield_General(), and GA_DecayYield_General().

184  {
186 
187  // Set the daFunk object based on the type of process.
188  daFunk::Funk Yield;
189  if (process.isAnnihilation)
190  {
191  Yield = daFunk::zero("E", "v");
192  }
193  else
194  {
195  Yield = daFunk::zero("E");
196  }
197 
198  // Adding two-body channels
199  for (std::vector<TH_Channel>::iterator it = process.channelList.begin();
200  it != process.channelList.end(); ++it)
201  {
202  bool added = false; // If spectrum is not available from any source
203 
204  // Here only take care of two-body final states
205  if (it->nFinalStates != 2) continue;
206 
207  // Get final state masses
208  double m0 = catalog.getParticleProperty(it->finalStateIDs[0]).mass;
209  double m1 = catalog.getParticleProperty(it->finalStateIDs[1]).mass;
210 
211  // Ignore channels that are kinematically closed for v=0
212  if ( m0 + m1 > Ecm ) continue;
213 
214  // Ignore channels with 0 BR in v=0 limit (if "v" is a variable of genRate, i.e. not a decay).
215  if (it->genRate->hasArg("v") && it->genRate->bind("v")->eval(0.) <= 0.0) continue;
216  else if ( !(it->genRate->hasArgs()) && it->genRate->bind()->eval() <=0.0) continue;
217 
218  double E0 = 0.5*(Ecm*Ecm+m0*m0-m1*m1)/Ecm;
219  double E1 = Ecm-E0;
220 
221  // Check whether two-body final state is in SimYield table
222  if ( table.hasChannel(it->finalStateIDs[0], it->finalStateIDs[1], "gamma") )
223  {
224  Yield = Yield +
225  it->genRate*(table)(
226  it->finalStateIDs[0], it->finalStateIDs[1], "gamma", Ecm);
227  added = true;
228  }
229  // Deal with composite final states
230  else
231  {
232  daFunk::Funk spec0 = daFunk::zero("E");
233  daFunk::Funk spec1 = daFunk::zero("E");
234  added = true;
235 
236  // Final state particle one
237  // Tabulated spectrum available?
238  if ( table.hasChannel(it->finalStateIDs[0], "gamma") )
239  {
240  spec0 = (table)(it->finalStateIDs[0], "gamma")->set("Ecm",E0);
241  }
242  // Gamma-ray line?
243  else if ( it->finalStateIDs[0] == "gamma" )
244  {
245  daFunk::Funk E = daFunk::var("E");
246  spec0 = exp(-pow((E-E0)/line_width/E0,2)/2)/E0/sqrt(2*M_PI)/line_width;
247  }
248  // MC spectra available?
249  else if ( cascadeMC_gammaSpectra.count(it->finalStateIDs[0]) )
250  {
251  double gamma0 = E0/m0;
252  //std::cout << it->finalStateIDs[0] << " " << gamma0 << std::endl;
253  spec0 = boost_dNdE(cascadeMC_gammaSpectra.at(it->finalStateIDs[0]), gamma0, 0.0);
254  }
255  else added = false;
256 
257  // Final state particle two
258  if ( table.hasChannel(it->finalStateIDs[1], "gamma") )
259  {
260  spec1 = (table)(it->finalStateIDs[1], "gamma")->set("Ecm", E1);
261  }
262  else if ( it->finalStateIDs[1] == "gamma" )
263  {
264  daFunk::Funk E = daFunk::var("E");
265  spec1 = exp(-pow((E-E1)/line_width/E1,2)/2)/E1/sqrt(2*M_PI)/line_width;
266  }
267  else if ( cascadeMC_gammaSpectra.count(it->finalStateIDs[1]) )
268  {
269  double gamma1 = E1/m1;
270  //std::cout << it->finalStateIDs[1] << " " << gamma1 << std::endl;
271  spec1 = boost_dNdE(cascadeMC_gammaSpectra.at(it->finalStateIDs[1]), gamma1, 0.0);
272  }
273  else added = false;
274 
275  #ifdef DARKBIT_DEBUG
276  std::cout << it->finalStateIDs[0] << " " << it->finalStateIDs[1] << std::endl;
277  //std::cout << "gammas: " << gamma0 << ", " << gamma1 << std::endl;
278  daFunk::Funk chnSpec = (daFunk::zero("v", "E")
279  + spec0
280  + spec1)-> set("v", 0.);
281  auto x = daFunk::logspace(0, 3, 10);
282  std::vector<double> y = chnSpec->bind("E")->vect(x);
283  std::cout << it->finalStateIDs[0] << it->finalStateIDs[1] << ":\n";
284  std::cout << " E: [";
285  for (std::vector<double>::iterator it2 = x.begin(); it2 != x.end(); it2++)
286  std::cout << *it2 << ", ";
287  std::cout << "]\n";
288  std::cout << " dNdE: [";
289  for (std::vector<double>::iterator it2 = y.begin(); it2 != y.end(); it2++)
290  std::cout << *it2 << ", ";
291  std::cout << "]\n";
292  #endif
293 
294  if (!added)
295  {
296  DarkBit_warning().raise(LOCAL_INFO,
297  "GA_AnnYield_General: cannot find spectra for "
298  + it->finalStateIDs[0] + " " + it->finalStateIDs[1]);
299  }
300 
301  Yield = Yield + (spec0 + spec1) * it->genRate;
302  }
303  } // End adding two-body final states
304 
305  #ifdef DARKBIT_DEBUG
306  std::vector<std::string> test1 = initVector<std::string> ("h0_1_test","h0_2_test","h0_2_test","h0_1_test","WH_test", "A0_test", "h0_1_test", "W+");
307  std::vector<std::string> test2 = initVector<std::string> ("A0_test", "A0_test", "Z0_test", "Z0_test", "WH_test", "Z0_test", "h0_2_test", "W-");
308 
309  for(size_t i=0; i<test1.size();i++)
310  {
311  daFunk::Funk chnSpec = (table)(test1[i], test2[i], "gamma", Ecm);
312  std::vector<double> y = chnSpec->bind("E")->vect(x);
313  os << test1[i] << test2[i] << ":\n";
314  os << " E: [";
315  for (std::vector<double>::iterator it2 = x.begin(); it2 != x.end(); it2++)
316  os << *it2 << ", ";
317  os << "]\n";
318  os << " dNdE: [";
319  for (std::vector<double>::iterator it2 = y.begin(); it2 != y.end(); it2++)
320  os << *it2 << ", ";
321  os << "]\n";
322  }
323  #endif
324 
325  // Adding three-body final states
326  //
327  // NOTE: Three body processes are added even if they are closed for at v=0
328  for (std::vector<TH_Channel>::iterator it = process.channelList.begin();
329  it != process.channelList.end(); ++it)
330  {
331  bool added = true;
332 
333  // Here only take care of three-body final states
334  if (it->nFinalStates != 3) continue;
335 
336  // Implement tabulated three-body final states
337  /*
338  if ( it->nFinalStates == 3
339  and table->hasChannel(it->finalStateIDs[0], "gamma")
340  and table->hasChannel(it->finalStateIDs[1], "gamma")
341  and table->hasChannel(it->finalStateIDs[2], "gamma")
342  )
343  {
344  daFunk::Funk dNdE1dE2 = it->genRate->set("v",0.);
345  daFunk::Funk spec0 =
346  (table)(it->finalStateIDs[0], "gamma");
347  daFunk::Funk spec1 =
348  (table)(it->finalStateIDs[1], "gamma");
349  daFunk::Funk spec2 =
350  (table)(it->finalStateIDs[2], "gamma");
351  Yield = Yield + convspec(spec0, spec1, spec2, dNdE1dE2);
352  }
353  */
354 
355  if ( it->finalStateIDs[0] == "gamma" )
356  {
357  if ( it->finalStateIDs[1] == "gamma" or it->finalStateIDs[2] == "gamma")
358  {
359  DarkBit_warning().raise(LOCAL_INFO, "Second and/or third primary gamma rays in three-body final states ignored.");
360  }
361  double m1 = catalog.getParticleProperty(it->finalStateIDs[1]).mass;
362  double m2 = catalog.getParticleProperty(it->finalStateIDs[2]).mass;
364  mass, m1, m2);
366  mass, m1, m2);
367  daFunk::Funk dsigmavde = it->genRate->gsl_integration(
368  "E1", E1_low, E1_high);
369 
370  #ifdef DARKBIT_DEBUG
371  daFunk::Funk chnSpec = (daFunk::zero("v", "E") + dsigmavde)-> set("v", 0.);
372  std::vector<double> y = chnSpec->bind("E")->vect(x);
373  os << it->finalStateIDs[0] << it->finalStateIDs[1] << it->finalStateIDs[2] << ":\n";
374  os << " E: [";
375  for (std::vector<double>::iterator it2 = x.begin(); it2 != x.end(); it2++)
376  os << *it2 << ", ";
377  os << "]\n";
378  os << " dNdE: [";
379  for (std::vector<double>::iterator it2 = y.begin(); it2 != y.end(); it2++)
380  os << *it2 << ", ";
381  os << "]\n";
382  #endif
383 
384  Yield = Yield + dsigmavde;
385  }
386  else added = false;
387 
388  if (!added)
389  {
390  DarkBit_warning().raise(LOCAL_INFO,
391  "GA_AnnYield_General: ignoring final state "
392  + it->finalStateIDs[0] + " " + it->finalStateIDs[1] + " " + it->finalStateIDs[2]);
393  }
394  }
395  #ifdef DARKBIT_DEBUG
396  if(debug) os.close();
397  #endif
398 
399  return Yield;
400 
401  }
template double gamma3bdy_limits< 0 >(double, double, double, double)
#define LOCAL_INFO
Definition: local_info.hpp:34
Funk var(std::string arg)
Definition: daFunk.hpp:921
daFunk::Funk boost_dNdE(daFunk::Funk dNdE, double gamma, double mass)
Boosts an energy spectrum of isotropic particles into another frame (and isotropizes again)...
Definition: GamYields.cpp:154
std::vector< double > logspace(double x0, double x1, unsigned int n)
Definition: daFunk.hpp:186
Funk func(double(*f)(funcargs...), Args... args)
Definition: daFunk.hpp:768
template double gamma3bdy_limits< 1 >(double, double, double, double)
double gamma3bdy_limits(double Eg, double M_DM, double m1, double m2)
Calculate kinematical limits for three-body final states.
warning & DarkBit_warning()
double pow(const double &a)
Outputs a^i.
Funk zero(Args... argss)
Definition: daFunk.hpp:604
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
double E(const double x, const double y)
Here is the call graph for this function: