gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
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  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 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 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)
 
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 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...
 
void GA_AnnYield_General (daFunk::Funk &result)
 General routine to derive annihilation yield. 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 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 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 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 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 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 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.

Function Documentation

◆ ALPS1_lnL_general()

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

Definition at line 1230 of file Axions.cpp.

Referenced by calc_lnL_ALPS1().

1231  {
1232  // Propagate uncertainty from efficiency in chi^2-likelihood.
1233  return -0.5*gsl_pow_2(s-mu)/(gsl_pow_2(0.05*s/0.82)+gsl_pow_2(sigma));
1234  }
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 1183 of file Axions.cpp.

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

Referenced by calc_ALPS1_signal_gas(), and calc_ALPS1_signal_vac().

1184  {
1185  const double eVm = gev2cm*1E7;
1186  // Photon energy in eV.
1187  const double erg = 2.0*pi*eVm/532.0E-9;
1188  const double len = 4.2;
1189  // We include the uncertainty of the detection efficiency eff = 0.82(5) in the likelihood.
1190  const double eff = 0.82;
1191 
1192  double result = 0.0;
1193 
1194  // CAVE: Approximations/conversion only valid/possible for m_ax << 2.33 eV (532 nm).
1195  if (m_ax < 1.0)
1196  {
1197  // Effective photon mass and momentum transfer.
1198  double m_ph_sq = 2.0*erg*erg*nm1;
1199  double q = 0.5*(m_ax*m_ax+m_ph_sq)/(eVm*erg);
1200  double factor = gsl_pow_4((gagg*1E17)*gsl_sf_sinc(0.5*q*len/pi));
1201 
1202  // Prefactor: 1096 W * 1 h * (10^-17/eV * 4.98 T * 4.2 m)^4 / 16.
1203  result = 0.00282962979*eff*factor*(power/1096.0)/erg;
1204  }
1205 
1206  return result;
1207  }
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 732 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.

733  {
734  const double eVm = gev2cm*1E7;
735  const double L = 9.26/eVm;
736  struct SolarModel_params4 * p4 = (struct SolarModel_params4 *)params;
737  double m_ax = p4->ma0;
738 
739  double argument = 0.25*1.0E-3*L*m_ax*m_ax/erg;
740  double temp = gsl_pow_2(gsl_sf_sinc(argument/pi));
741  double exposure = p4->eff_exp->interpolate(erg);
742  double gaee_flux = p4->gaee_flux->interpolate(erg);
743 
744  return temp*exposure*gaee_flux;
745  }
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 210 of file SunNeutrinos.cpp.

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

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

211  {
212  using namespace Pipes::annihilation_rate_Sun;
213  double tt_sun = 1.5e17 / *Dep::equilibration_time_Sun;
214  result = *Dep::capture_rate_Sun * 0.5 * pow(tanh(tt_sun),2.0);
215  }
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 1638 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().

1639  {
1640  double res = 1.0;
1641  if (T > Tchi) { res = pow(T/Tchi,-0.5*beta); }
1642  return res;
1643  }
START_MODEL Tchi
Definition: Axions.hpp:35
START_MODEL beta
Definition: Axions.hpp:35
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 150 of file GamYields.cpp.

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

Referenced by GA_AnnYield_General().

151  {
152  if ( gamma < 1.0 + .02 ) // Ignore less than 2% boosts
153  {
154  if (gamma < 1.0)
155  DarkBit_error().raise(LOCAL_INFO,
156  "boost_dNdE: Requested Lorentz boost with gamma < 1");
157  return dNdE;
158  }
159  double betaGamma = sqrt(gamma*gamma-1);
160  daFunk::Funk E = daFunk::var("E");
161  daFunk::Funk lnE = daFunk::var("lnE");
162  daFunk::Funk Ep = daFunk::var("Ep");
163  daFunk::Funk halfBox_int = betaGamma*sqrt(E*E-mass*mass);
164  daFunk::Funk halfBox_bound = betaGamma*sqrt(Ep*Ep-mass*mass);
165  daFunk::Funk integrand = dNdE/(2*halfBox_int);
166  return integrand->gsl_integration("E", Ep*gamma-halfBox_bound, Ep*gamma+halfBox_bound)
167  ->set_epsabs(0)->set_limit(100)->set_epsrel(1e-3)->set_use_log_fallback(true)->set("Ep", daFunk::var("E"));
168  //
169  // Note: integration over lnE causes problems in the WIMP example (3) as the singularity is dropped.
170  // return (integrand*E)->set("E", exp(lnE))->gsl_integration("lnE", log(Ep*gamma-halfBox_bound), log(Ep*gamma+halfBox_bound))
171  // ->set_epsabs(0)->set_epsrel(1e-3)->set("Ep", daFunk::var("E"));
172  }
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 1220 of file Axions.cpp.

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

1221  {
1222  using namespace Pipes::calc_ALPS1_signal_gas;
1223  double m_ax = *Param["ma0"];
1224  double gagg = 1.0E-9*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1225 
1226  result = ALPS1_signal_general(1044.0, 5.0E-8, m_ax, gagg);
1227  }
void calc_ALPS1_signal_gas(double &result)
Definition: Axions.cpp:1220
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:1183
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 1210 of file Axions.cpp.

References ALPS1_signal_general().

1211  {
1212  using namespace Pipes::calc_ALPS1_signal_vac;
1213  double m_ax = *Param["ma0"];
1214  double gagg = 1.0E-9*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1215 
1216  result = ALPS1_signal_general(1096.0, 0.0, m_ax, gagg);
1217  }
void calc_ALPS1_signal_vac(double &result)
Definition: Axions.cpp:1210
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:1183
Here is the call graph for this function:

◆ calc_AxionOscillationTemperature()

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

Definition at line 1666 of file Axions.cpp.

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

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

References Gambit::Scanner::pow().

1266  {
1267  using namespace Pipes::calc_CAST2007_signal_vac;
1268  double m_ax = *Param["ma0"];
1269  double gagg = 1.0E-9*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1270  double gaee = std::fabs(*Param["gaee"]);
1271 
1272  // Initialise the Solar model calculator and get the reference counts for a given mass.
1273  // Get Solar model we are working with; set default value here
1274  static std::string solar_model_gagg = runOptions->getValueOrDef<std::string> ("AGSS09met", "solar_model_gagg");
1275  static std::string solar_model_gaee = runOptions->getValueOrDef<std::string> ("AGSS09met_old", "solar_model_gaee");
1276  static CAST_SolarModel_Interpolator lg_ref_counts (solar_model_gagg, solar_model_gaee, "CAST2007");
1277  std::vector<double> lg_ref_counts_gagg = lg_ref_counts.evaluate_gagg_contrib(m_ax);
1278  std::vector<double> lg_ref_counts_gaee = lg_ref_counts.evaluate_gaee_contrib(m_ax);
1279  static int n_bins = lg_ref_counts_gagg.size();
1280 
1281  std::vector<double> counts;
1282  double dummy;
1283  for (int i = 0; i < n_bins; i++)
1284  {
1285  dummy = gsl_pow_2(gagg*1E19)*pow(10,lg_ref_counts_gagg[i]) + gsl_pow_2(gaee*1E13)*pow(10,lg_ref_counts_gaee[i]);
1286  counts.push_back(gsl_pow_2(gagg*1E19)*dummy);
1287  }
1288 
1289  result = counts;
1290  }
void calc_CAST2007_signal_vac(std::vector< double > &result)
Definition: Axions.cpp:1265
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 1293 of file Axions.cpp.

References Gambit::Scanner::pow().

1294  {
1295  using namespace Pipes::calc_CAST2017_signal_vac;
1296  double m_ax = *Param["ma0"];
1297  double gagg = 1.0E-9*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1298  double gaee = std::fabs(*Param["gaee"]);
1299  std::vector<std::vector<double>> res;
1300 
1301  // Initialise the Solar model calculator and get the reference counts for a given mass.
1302  // Get Solar model we are working with; set default value here
1303  static std::string solar_model_gagg = runOptions->getValueOrDef<std::string> ("AGSS09met", "solar_model_gagg");
1304  static std::string solar_model_gaee = runOptions->getValueOrDef<std::string> ("AGSS09met_old", "solar_model_gaee");
1305 
1306  const int n_exps = 12;
1307  const int n_bins = 10;
1308  const std::string exp_names [n_exps] = {"A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"};
1309  static std::vector<CAST_SolarModel_Interpolator> lg_ref_counts;
1310  static bool lg_ref_counts_not_calculated = true;
1311  if (lg_ref_counts_not_calculated)
1312  {
1313  for (int e = 0; e < n_exps; e++)
1314  {
1315  CAST_SolarModel_Interpolator dummy (solar_model_gagg, solar_model_gaee, "CAST2017_"+exp_names[e]);
1316  lg_ref_counts.push_back(std::move(dummy));
1317  }
1318  }
1319  lg_ref_counts_not_calculated = false;
1320 
1321  for (int e = 0; e < n_exps; e++)
1322  {
1323  std::vector<double> lg_ref_counts_gagg = lg_ref_counts[e].evaluate_gagg_contrib(m_ax);
1324  std::vector<double> lg_ref_counts_gaee = lg_ref_counts[e].evaluate_gaee_contrib(m_ax);
1325 
1326  std::vector<double> counts;
1327  double dummy;
1328  for (int bin = 0; bin < n_bins; bin++)
1329  {
1330  dummy = gsl_pow_2(gagg*1E19)*pow(10,lg_ref_counts_gagg[bin]) + gsl_pow_2(gaee*1E13)*pow(10,lg_ref_counts_gaee[bin]);
1331  counts.push_back(gsl_pow_2(gagg*1E19)*dummy);
1332  }
1333 
1334  res.push_back(counts);
1335  }
1336 
1337  result = res;
1338  }
void calc_CAST2017_signal_vac(std::vector< std::vector< double >> &result)
Definition: Axions.cpp:1293
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 1438 of file Axions.cpp.

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

1439  {
1440  using namespace Pipes::calc_Haloscope_signal;
1441  double gagg = 1.0E-9*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
1442  // Get the DM fraction in axions and the local DM density.
1443  double fraction = *Dep::RD_fraction;
1444  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
1445  double rho0 = LocalHaloParameters.rho0;
1446 
1447  // Signal relative to a reference coupling and local DM density.
1448  double s = gsl_pow_2(gagg/1.0E-24) * fraction * (rho0/0.45);
1449 
1450  result = s;
1451  }
void calc_Haloscope_signal(double &result)
Definition: Axions.cpp:1438
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 1237 of file Axions.cpp.

References ALPS1_lnL_general().

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

◆ calc_lnL_CAST2007()

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

Definition at line 1356 of file Axions.cpp.

References CAST_lnL_general().

1357  {
1358  using namespace Pipes::calc_lnL_CAST2007;
1359  std::vector<double> sig_vac = *Dep::CAST2007_signal_vac;
1360 
1361  // CAST CCD 2004 vacuum data (based on hep-ex/0702006).
1362  const int n_bins = 20;
1363  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};
1364  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,
1365  1.247346181, 1.247346181, 2.286801272, 1.247346181};
1366 
1367  // Only calculate norm once.
1368  static double norm = 0.0;
1369  static bool norm_not_calculated = true;
1370  if (norm_not_calculated)
1371  {
1372  for (int i = 0; i < n_bins; i++) { norm += gsl_sf_lnfact(dat_vac[i]); }
1373  }
1374  norm_not_calculated = false;
1375 
1376  result = CAST_lnL_general(sig_vac, bkg_vac, dat_vac) - norm;
1377  }
void calc_lnL_CAST2007(double &result)
Definition: Axions.cpp:1356
double CAST_lnL_general(std::vector< double > s, const std::vector< double > bkg_counts, const std::vector< int > sig_counts)
Definition: Axions.cpp:1341
Here is the call graph for this function:

◆ calc_lnL_CAST2017()

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

Definition at line 1380 of file Axions.cpp.

References CAST_lnL_general().

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

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

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

◆ calc_lnL_Haloscope_ADMX2()

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

Definition at line 1496 of file Axions.cpp.

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

◆ calc_lnL_Haloscope_RBF()

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

Definition at line 1547 of file Axions.cpp.

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

1548  {
1549  using namespace Pipes::calc_lnL_Haloscope_RBF;
1550  // Rescale the axion mass to ueV.
1551  double m = 1.0E+6*(*Param["ma0"]);
1552  double l = 0.0;
1553  // Results from Phys. Rev. D 40, 3153 (1989)
1554  // The bins and results below (from Table I) appear to be inconsitent with Figure 14.
1555  //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,
1556  // 11.298638, 11.583999, 12.845377, 13.234130, 15.301962, 16.2655809};
1557  // N_sigma values as defined in the paper.
1558  //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};
1559  // Proportionality factors ("sigma") inferred from Table I in in units of GeV^-1/cm^3.
1560  //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,
1561  // 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};
1562  // The results below are derived from Fig. 14 (assuming N_sigma = 4 for all values).
1563  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};
1564  const double N_sigma = 4.0;
1565  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};
1566 
1567  // Likelihood shape parameters based on information from PRD 40 (1989) 3153.
1568  // 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.
1569  if ( ((m > bins.front()) && (m < 10.067396)) || ((m > 11.213613) && (m < bins.back())))
1570  {
1571  // Use the standard search algorthim to identify the bin index and use the appropriate values for the likelihood.
1572  auto index = upper_bound(bins.begin(), bins.end(), m);
1573  double sigma = eff_sigma[index-bins.begin()-1];
1574  // Uncomment and comment out lines below to swap between implementations using Table I and Figure 14, respectively.
1575  //double offset = sigma*N_sigma[index-bins.begin()-1];
1576  double offset = sigma*N_sigma;
1577  // The "signal" needs to be rescaled to 0.3 GeV/cm^3, which is their reference value.
1578  double s = (0.45/0.3)*(*Dep::Haloscope_signal);
1579  if (s > offset) {
1580  l = -0.5 * gsl_pow_2( (s - offset)/sigma );
1581  }
1582  }
1583 
1584  result = l;
1585  }
const double sigma
Definition: SM_Z.hpp:43
void calc_lnL_Haloscope_RBF(double &result)
Definition: Axions.cpp:1547

◆ calc_lnL_Haloscope_UF()

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

Definition at line 1523 of file Axions.cpp.

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

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

◆ calc_lnL_HESS_GCMF()

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

Definition at line 2032 of file Axions.cpp.

References daFunk::interp().

2033  {
2034  using namespace Pipes::calc_lnL_HESS_GCMF;
2035  double m_ax = *Param["ma0"];
2036  double gagg = 1.0E-9*std::fabs(*Param["gagg"]); // gagg needs to be in eV^-1.
2037 
2038  // Compute the domensionless parameters Epsilon and Gamma from the axion mass and axion-photon coupling (see 1311.3148).
2039  const double c_epsilon = 0.071546787;
2040  const double c_gamma = 0.015274036*370.0/sqrt(37.0);
2041  double epsilon = log10(m_ax*c_epsilon) + 5.0;
2042  double gamma = log10(gagg*c_gamma) + 20.0;
2043 
2044  // Initialise the interpolation and extrapolation routies for the H.E.S.S. results.
2045  static HESS_Interpolator interp (GAMBIT_DIR "/DarkBit/data/HESS_GCMF_Table.dat");
2046 
2047  result = interp.lnL(epsilon, gamma);
2048  }
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:2032
Here is the call graph for this function:

◆ calc_lnL_RParameter()

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

Definition at line 1887 of file Axions.cpp.

1888  {
1889  using namespace Pipes::calc_lnL_RParameter;
1890  double Rtheo = *Dep::RParameter;
1891 
1892  // Observed R values from astro-ph/0403600.
1893  const double Robs = 1.39;
1894  const double RobsErr = 0.03;
1895  // Value for He-abundance Y from 1503.08146: <Y> = 0.2515(17).
1896  const double YobsErrContrib = 7.3306*0.0017;
1897 
1898  result = -0.5*gsl_pow_2(Rtheo - Robs)/(RobsErr*RobsErr+YobsErrContrib*YobsErrContrib);
1899  }
void calc_lnL_RParameter(double &result)
Definition: Axions.cpp:1887

◆ calc_lnL_SN1987A()

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

Definition at line 1998 of file Axions.cpp.

1999  {
2000  using namespace Pipes::calc_lnL_SN1987A;
2001  double f_10s = *Dep::PhotonFluence_SN1987A_Conversion;
2002 
2003  // Standard devations of the null observation.
2004  const double sigma_10s = 0.2;
2005 
2006  double ratio = f_10s/sigma_10s;
2007 
2008  result = -0.5*ratio*ratio;
2009  }
void calc_lnL_SN1987A(double &result)
Definition: Axions.cpp:1998

◆ calc_lnL_WDVar_G117B15A()

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

Definition at line 1906 of file Axions.cpp.

References Gambit::LogTags::err.

1907  {
1908  using namespace Pipes::calc_lnL_WDVar_G117B15A;
1909  // Rescale coupling to be used in their model prediction.
1910  double x = (1.0E+14 * std::fabs(*Param["gaee"]))/2.8;
1911  // We only have predictions up to x = 30. Limits should get stronger for x > 30, so
1912  // it is conservative to use the prediction for x = 30 for x > 30.
1913  x = std::min(x,30.0);
1914 
1915  // Values for the model prediction provided by the authors.
1916  const std::vector<double> xvals = {0.0, 1.0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 17.5, 20.1, 22.5, 25.0, 27.5, 30.0};
1917  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};
1918  const double err = 0.09;
1919 
1920  // Use interpolation for the model predction, but only initialise once.
1921  static AxionInterpolator period_change (xvals, dPidts, "cspline");
1922 
1923  double pred = period_change.interpolate(x);
1924  result = -0.5 * gsl_pow_2(4.19 - pred) / (0.73*0.73 + err*err);
1925  }
void calc_lnL_WDVar_G117B15A(double &result)
Definition: Axions.cpp:1906

◆ calc_lnL_WDVar_L192()

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

Definition at line 1972 of file Axions.cpp.

References Gambit::LogTags::err.

1973  {
1974  using namespace Pipes::calc_lnL_WDVar_L192;
1975  // Rescale coupling to be used in their model prediction.
1976  double x = (1.0E+14 * std::fabs(*Param["gaee"]))/2.8;
1977  // We only have predictions up to x = 30. Limits should get stronger for x > 30, so
1978  // it is conservative to use the prediction for x = 30 for x > 30.
1979  x = std::min(x,23.0);
1980 
1981  // Values for the model prediction provided by the authors.
1982  const std::vector<double> xvals = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0};
1983  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};
1984  const double err = 0.85;
1985 
1986  static AxionInterpolator period_change (xvals, dPidts, "cspline");
1987 
1988  double pred = period_change.interpolate(x);
1989  result = -0.5 * gsl_pow_2(3.0 - pred) / (0.6*0.6 + err*err);
1990  }
void calc_lnL_WDVar_L192(double &result)
Definition: Axions.cpp:1972

◆ calc_lnL_WDVar_PG1351489()

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

Definition at line 1949 of file Axions.cpp.

References Gambit::LogTags::err.

1950  {
1951  using namespace Pipes::calc_lnL_WDVar_PG1351489;
1952  // Rescale coupling to be used in their model prediction.
1953  double x = (1.0E+14 * std::fabs(*Param["gaee"]))/2.8;
1954  // We only have predictions up to x = 20. Limits should get stronger for x > 20, so
1955  // it is conservative to use the prodiction for x = 20 for x > 20.
1956  x = std::min(x,20.0);
1957 
1958  // Values for the model prediction provided by the authors.
1959  const std::vector<double> xvals = {0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0};
1960  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};
1961  const double err = 0.5;
1962 
1963  static AxionInterpolator period_change (xvals, dPidts, "cspline");
1964 
1965  // We only have predictions up to x = 20. Limits should get stronger for x > 20, so
1966  // it is conservative to use the prodiction for x = 20 for x > 20.
1967  double pred = period_change.interpolate(x);
1968  result = -0.5 * gsl_pow_2(2.0 - pred) / (0.9*0.9 + err*err);
1969  }
void calc_lnL_WDVar_PG1351489(double &result)
Definition: Axions.cpp:1949

◆ calc_lnL_WDVar_R548()

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

Definition at line 1928 of file Axions.cpp.

References Gambit::LogTags::err.

1929  {
1930  using namespace Pipes::calc_lnL_WDVar_R548;
1931  // Rescale coupling to be used in their model prediction.
1932  double x = (1.0E+14 * std::fabs(*Param["gaee"]))/2.8;
1933  // We only have predictions up to x = 30. Limits should get stronger for x > 30, so
1934  // it is conservative to use the prodiction for x = 30 for x > 30.
1935  x = std::min(x,30.0);
1936 
1937  // Values for the model prediction provided by the authors.
1938  const std::vector<double> xvals = {0.0, 1.0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 17.5, 20.0, 22.5, 25.0, 27.5, 30.0};
1939  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};
1940  const double err = 0.09;
1941 
1942  static AxionInterpolator period_change (xvals, dPidts, "cspline");
1943 
1944  double pred = period_change.interpolate(x);
1945  result = -0.5 * gsl_pow_2(3.3 - pred) / (1.1*1.1 + err*err);
1946  }
void calc_lnL_WDVar_R548(double &result)
Definition: Axions.cpp:1928

◆ calc_PhotonFluence_SN1987A_Conversion()

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

Definition at line 2017 of file Axions.cpp.

References Gambit::Scanner::pow().

2018  {
2020  double m = (1.0E+10*(*Param["ma0"]))/5.433430;
2021  double g = (1.0E+12*std::fabs(*Param["gagg"]))/5.339450;
2022 
2023  result = 0.570589*gsl_pow_4(g);
2024  if (m > 1.0) { result = result*pow(m, -4.021046); }
2025  }
double pow(const double &a)
Outputs a^i.
void calc_PhotonFluence_SN1987A_Conversion(double &result)
Definition: Axions.cpp:2017
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 1863 of file Axions.cpp.

References Gambit::FlavBit::LoopFunctions::E(), and Gambit::Scanner::pow().

1864  {
1865  using namespace Pipes::calc_RParameter;
1866  double gaee2 = gsl_pow_2(1.0E+13*std::fabs(*Param["gaee"]));
1867  double gagg = 1.0E+10*std::fabs(std::fabs(*Param["gagg"]));
1868  double lgma0 = log10(*Param["ma0"]);
1869  // Value for He-abundance Y from 1503.08146: <Y> = 0.2515(17).
1870  const double Y = 0.2515;
1871  // Use interpolation for the finite-mass correction.
1872  static AxionInterpolator correction (GAMBIT_DIR "/DarkBit/data/Axions_RParameterCorrection.dat", "linear");
1873  // Initialise an effective axion-photon coupling, valid for low masses.
1874  double geff = gagg;
1875  // Apply correction for higher mass values...
1876  if ((lgma0 > correction.lower()) && (lgma0 < correction.upper())) { geff *= pow(10, 0.5*correction.interpolate(lgma0)); }
1877  // ... or set to zero if mass is too high.
1878  if (lgma0 >= correction.upper()) { geff = 0.0; }
1879  // Expressions only valid for gaee2 < 35.18 but limits should become stronger for gaee2 > 35.18 (but perhaps not gaee2 >> 35.18).
1880  // Conservative approach: Constrain gaee2 > 35.18 at the level of gaee2 = 35.18.
1881  if (gaee2 > 35.18) { gaee2 = 35.18; }
1882 
1883  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;
1884  }
void calc_RParameter(double &result)
Capabilities relating to astrophysical observations (R-parameter, H.E.S.S. telescope search...
Definition: Axions.cpp:1863
double pow(const double &a)
Outputs a^i.
double E(const double x, const double y)
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:99
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(), combine_hdf5::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:527
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:587
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:527
#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:587
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:544
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:473
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:35
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:527
#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 1341 of file Axions.cpp.

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

Referenced by calc_lnL_CAST2007(), and calc_lnL_CAST2017().

1342  {
1343  double result = 0.0;
1344  int n_bins = s.size();
1345 
1346  for (int i = 0; i < n_bins; i++)
1347  {
1348  double mu = s[i] + bkg_counts[i];
1349  result += sig_counts[i]*gsl_sf_log(mu) - mu;
1350  }
1351 
1352  return result;
1353  }
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 743 of file MSSM.cpp.

Referenced by main().

744  {
745  using namespace Pipes::DarkMatter_ID_MSSM;
746  // TODO: need ask Dep::MSSM_spectrum in future; might have sneutralinos and gravitinos later.
747  result = "~chi0_1";
748  }
void DarkMatter_ID_MSSM(std::string &result)
Definition: MSSM.cpp:743
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 199 of file DarkBit_standalone_WIMP.cpp.

200  {
201  result = "WIMP";
202  }

◆ 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 705 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().

706  {
708 
709  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
710 
711  double rho0 = LocalHaloParameters.rho0;
712  double rho0_eff = (*Dep::RD_fraction)*rho0;
713  double vrot = LocalHaloParameters.vrot;
714  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
715  double vesc = LocalHaloParameters.vesc;
717  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
718 
719  BEreq::dshmcom->rho0 = rho0;
720  BEreq::dshmcom->v_sun = vrot;
721  BEreq::dshmcom->v_earth = v_earth;
722  BEreq::dshmcom->rhox = rho0_eff;
723 
724  BEreq::dshmframevelcom->v_obs = vrot;
725 
726  BEreq::dshmisodf->vd_3d = vd_3d;
727  BEreq::dshmisodf->vgalesc = vesc;
728 
729  BEreq::dshmnoclue->vobs = vrot;
730 
732  << "Updating DarkSUSY halo parameters:" << std::endl
733  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
734  << " rho0_eff [GeV/cm^3] = " << rho0_eff << std::endl
735  << " v_sun [km/s] = " << vrot<< std::endl
736  << " v_earth [km/s] = " << v_earth << std::endl
737  << " v_obs [km/s] = " << vrot << std::endl
738  << " vd_3d [km/s] = " << vd_3d << std::endl
739  << " v_esc [km/s] = " << vesc << EOM;
740 
741  result = true;
742 
743  return;
744  }
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:99
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 747 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().

748  {
750 
751  LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
752 
753  double rho0 = LocalHaloParameters.rho0;
754  double vrot = LocalHaloParameters.vrot;
755  double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
756  double vesc = LocalHaloParameters.vesc;
758  double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");
759 
760  BEreq::dshmcom->rho0 = rho0;
761  BEreq::dshmcom->v_sun = vrot;
762  BEreq::dshmcom->v_earth = v_earth;
763 
764  BEreq::dshmframevelcom->v_obs = vrot;
765 
766  BEreq::dshmisodf->vd_3d = vd_3d;
767  BEreq::dshmisodf->vgalesc = vesc;
768 
769  BEreq::dshmnoclue->vobs = vrot;
770 
772  << "Updating DarkSUSY halo parameters:" << std::endl
773  << " rho0 [GeV/cm^3] = " << rho0 << std::endl
774  << " v_sun [km/s] = " << vrot<< std::endl
775  << " v_earth [km/s] = " << v_earth << std::endl
776  << " v_obs [km/s] = " << vrot << std::endl
777  << " vd_3d [km/s] = " << vd_3d << std::endl
778  << " v_esc [km/s] = " << vesc << EOM;
779 
780  result = true;
781 
782  return;
783  }
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:99
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:99
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:99
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 205 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().

206  {
208  const Spectrum& spec = *Dep::DiracSingletDM_Z2_spectrum;
209  const SubSpectrum& he = spec.get_HE();
210  //double mass = spec.get(Par::Pole_Mass,"F");
211  double lambda = he.get(Par::dimensionless,"lF");
212  double cosXI = std::cos(he.get(Par::dimensionless,"xi"));
213  double sinXI = std::sin(he.get(Par::dimensionless,"xi"));
214  double mh = spec.get(Par::Pole_Mass,"h0_1");
215 
216  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
217  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
218  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
219 
220  // SI scalar and pseudoscalar couplings
221  result.gps = lambda*fp*m_proton*cosXI/pow(mh,2);
222  result.gns = lambda*fn*m_neutron*cosXI/pow(mh,2);
223  result.gp_q2 = lambda*fp*m_proton*sinXI/pow(mh,2);
224  result.gn_q2 = lambda*fn*m_neutron*sinXI/pow(mh,2);
225 
226  logger() << LogTags::debug << "Dirac DM DD couplings:" << std::endl;
227  logger() << " gps = " << result.gps << std::endl;
228  logger() << " gns = " << result.gns << std::endl;
229  logger() << " gp_q2 = " << result.gp_q2 << std::endl;
230  logger() << " gn_q2 = " << result.gn_q2 << EOM;
231 
232  } // 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:99
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
Spectrum Spectrum
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 208 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().

209  {
212  const SubSpectrum& he = spec.get_HE();
213  //double mass = spec.get(Par::Pole_Mass,"X");
214  double lambda = he.get(Par::dimensionless,"lX");
215  double cosXI = std::cos(he.get(Par::dimensionless,"xi"));
216  double sinXI = std::sin(he.get(Par::dimensionless,"xi"));
217  double mh = spec.get(Par::Pole_Mass,"h0_1");
218 
219  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
220  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
221  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
222 
223  // SI scalar and pseudoscalar couplings
224  result.gps = lambda*fp*m_proton*cosXI/pow(mh,2);
225  result.gns = lambda*fn*m_neutron*cosXI/pow(mh,2);
226  result.gp_q2 = lambda*fp*m_proton*sinXI/pow(mh,2);
227  result.gn_q2 = lambda*fn*m_neutron*sinXI/pow(mh,2);
228 
229  logger() << LogTags::debug << "Majorana DM DD couplings:" << std::endl;
230  logger() << " gps = " << result.gps << std::endl;
231  logger() << " gns = " << result.gns << std::endl;
232  logger() << " gp_q2 = " << result.gp_q2 << std::endl;
233  logger() << " gn_q2 = " << result.gn_q2 << EOM;
234 
235  } // 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:99
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
Spectrum Spectrum
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:99
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 223 of file ScalarSingletDM.cpp.

References get_ScalarSingletDM_DD_couplings(), and ScalarSingletDM_Z2_spectrum.

Referenced by main().

224  {
227  get_ScalarSingletDM_DD_couplings(spec, result, Param);
228  }
void DD_couplings_ScalarSingletDM_Z2(DM_nucleon_couplings &result)
Direct detection couplings for Z2 scalar singlet DM.
Spectrum 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:
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 231 of file ScalarSingletDM.cpp.

References get_ScalarSingletDM_DD_couplings(), and ScalarSingletDM_Z3_spectrum.

232  {
235  get_ScalarSingletDM_DD_couplings(spec, result, Param);
236  }
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
Spectrum 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 195 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.

196  {
199  const SubSpectrum& he = spec.get_HE();
200  double mass = spec.get(Par::Pole_Mass,"V");
201  double lambda = he.get(Par::dimensionless,"lambda_hV");
202  double mh = spec.get(Par::Pole_Mass,"h0_1");
203 
204  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
205  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
206  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
207 
208  result.gps = lambda*fp*m_proton/pow(mh,2)/mass/2;
209  result.gns = lambda*fn*m_neutron/pow(mh,2)/mass/2;
210  result.gpa = 0; // Only SI cross-section
211  result.gna = 0;
212 
213  logger() << LogTags::debug << "Vector DM DD couplings:" << std::endl;
214  logger() << " gps = " << result.gps << std::endl;
215  logger() << " gns = " << result.gns << std::endl;
216  logger() << " gpa = " << result.gpa << std::endl;
217  logger() << " gna = " << result.gna << EOM;
218 
219  } // 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:99
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
Spectrum Spectrum
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 204 of file DarkBit_standalone_WIMP.cpp.

Referenced by main().

205  {
206  using namespace Pipes::DD_couplings_WIMP;
208  result.gps = runOptions->getValueOrDef<double>(0., "gps");
210  result.gns = runOptions->getValueOrDef<double>(0., "gns");
212  result.gpa = runOptions->getValueOrDef<double>(0., "gpa");
214  result.gna = runOptions->getValueOrDef<double>(0., "gna");
215  //std::cout << "DD_coupling says" << std::endl;
216  //std::cout << result.gps << std::endl;
217  }
void DD_couplings_WIMP(DM_nucleon_couplings &result)
Here is the caller graph for this function:

◆ DEF_FUNKTRAIT()

Gambit::DarkBit::DEF_FUNKTRAIT ( RD_EFF_ANNRATE_FROM_PROCESSCATALOG_TRAIT  ) &

Infer Weff from process catalog.

Definition at line 553 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().

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

◆ 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 412 of file SimpleLikelihoods.cpp.

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

Referenced by dumpSpectrum(), and main().

413  {
414  using namespace Pipes::dump_GammaSpectrum;
415  daFunk::Funk spectrum = (*Dep::GA_AnnYield)->set("v", 0.);
416  // Option filename<string>: Filename for gamma-ray spectrum dump
417  // (default: dNdE.dat)
418  std::string filename = runOptions->getValueOrDef<std::string>(
419  "dNdE.dat", "filename");
420  logger() << "FILENAME for gamma dump: " << filename << EOM;
421  std::ofstream myfile (filename);
422  if (myfile.is_open())
423  {
424  for (int i = 0; i<=1200; i++)
425  {
426  double energy = pow(10., i/200. - 4.);
427 
428  myfile << energy << " " << spectrum->bind("E")->eval(energy) << "\n";
429  }
430  myfile.close();
431  }
432  result = 0.;
433  }
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:99
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 1651 of file Axions.cpp.

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

Referenced by calc_AxionOscillationTemperature().

1652  {
1653  // T/MeV, ma0/eV, m_pl/10^12eV = m_pl/10^3 GeV
1654  const double m_pl = m_planck_red*1.0E-3;
1655  struct AxionEDT_params * p1 = (struct AxionEDT_params *)params;
1656  double ma0 = (p1->ma0);
1657  double beta = (p1->beta);
1658  double Tchi = (p1->Tchi);
1659 
1660  double result = 1.0 - gStar(T)*gsl_pow_2(T*T*pi/(ma0*m_pl*axion_mass_temp(T, beta, Tchi)))/10.0;
1661 
1662  return result;
1663  }
START_MODEL Tchi
Definition: Axions.hpp:35
START_MODEL beta
Definition: Axions.hpp:35
double gStar(double T)
Various capabilities and functions to provide SM physics as well as QCD input for axions...
Definition: Axions.cpp:1066
const double m_planck_red
const double pi
double axion_mass_temp(double T, double beta, double Tchi)
Definition: Axions.cpp:1638
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, DarkMatter_ID, Gambit::DarkBit::TH_Process::genRateMisc, 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  TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMid);
185 
186  // Add all the regular channels
187  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
188  it != annProc.channelList.end(); ++it)
189  {
190  if ( it->nFinalStates == 2 )
191  {
192  // (sv)(v = sqrt(2T/mDM)) for two-body final state
193  sigmav += it->genRate->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
194  }
195  }
196 
197  // Add invisible contributions
198  sigmav += annProc.genRateMisc->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
199 
200  double ca = sigmav/6.6e28 * pow(*Dep::mwimp/20.0, 1.5);
201  // Scale the annihilation rate down by a factor of two if the DM is not self-conjugate
202  if (not (*Dep::TH_ProcessCatalog).getProcess(*Dep::DarkMatter_ID, *Dep::DarkMatter_ID).isSelfConj) ca *= 0.5;
203  result = pow(*Dep::capture_rate_Sun * ca, -0.5);
204 
205  // std::cout << "v = " << sqrt(2.0*T_Sun_core/(*Dep::mwimp)) << " and sigmav inside equilibration_time_Sun = " << sigmav << std::endl;
206  // std::cout << "capture_rate_Sun inside equilibration_time_Sun = " << *Dep::capture_rate_Sun << std::endl;
207  }
void equilibration_time_Sun(double &result)
Equilibration time for capture and annihilation of dark matter in the Sun (s)
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 701 of file Axions.cpp.

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

Referenced by alt_erg_integrand().

702  {
703  const double eVm = gev2cm*1E7;
704  const double L = 9.26/eVm;
705  struct SolarModel_params3 * p3 = (struct SolarModel_params3 *)params;
706  SolarModel* sol = p3->sol;
707  double m_ax = p3->ma0;
708  double rs = p3->rs;
709 
710  double argument = 0.25*1.0E-3*L*m_ax*m_ax/erg;
711  double temp = gsl_pow_2(gsl_sf_sinc(argument/pi));
712  double exposure = p3->eff_exp->interpolate(erg);
713  //std::cout << "Energy: " << erg << ", expoure = " << exposure << "." << std::endl;
714  SolarModel_params2 p2 = {erg, rs, sol};
715 
716  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1E6);
717  double result, error;
718 
719  gsl_function F;
720  F.function = &rad_integrand;
721  F.params = &p2;
722 
723  // Max. and min. integration radius
724  double rmin = sol->r_lo, rmax = std::min(rs, sol->r_hi);
725 
726  gsl_integration_qag (&F, rmin, rmax, 1e-1*abs_prec, 1e-1*rel_prec, 1E6, method, w, &result, &error);
727  gsl_integration_workspace_free (w);
728 
729  return temp*exposure*result;
730  }
START_MODEL rs
const double pi
const int method
Definition: Axions.cpp:646
const double gev2cm
double rad_integrand(double rad, void *params)
Definition: Axions.cpp:676
const double rel_prec
Definition: Axions.cpp:645
const double abs_prec
Definition: Axions.cpp:645
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 152 of file DarkBit.cpp.

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

Referenced by main().

153  {
154  using namespace Pipes::ExtractLocalMaxwellianHalo;
155  double rho0 = *Param["rho0"];
156  double v0 = *Param["v0"];
157  double vesc = *Param["vesc"];
158  double vrot = *Param["vrot"];
159  result.rho0 = rho0;
160  result.v0 = v0;
161  result.vesc = vesc;
162  result.vrot = vrot;
163  }
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:152
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 191 of file GamYields.cpp.

References boost_dNdE(), cascadeMC_gammaSpectra, DarkBit_warning(), DarkMatter_ID, 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 >(), daFunk::ifelse(), LOCAL_INFO, daFunk::logspace(), Gambit::Scanner::pow(), daFunk::throwError(), daFunk::var(), and daFunk::zero().

Referenced by dumpSpectrum(), and main().

192  {
193  using namespace Pipes::GA_AnnYield_General;
195 
196  std::string DMid= *Dep::DarkMatter_ID;
197 
199  double line_width = runOptions->getValueOrDef<double>(0.03, "line_width");
200 
201  // Get annihilation process from process catalog
202  TH_Process annProc = (*Dep::TH_ProcessCatalog).getProcess(DMid, DMid);
203 
204  // If process involves non-self-conjugate DM then we need to add a factor of 1/2 to the final spectrum. This must be explicitly set in the process catalogue.
205  double k = (annProc.isSelfConj) ? 1. : 0.5;
206 
207  // Get particle mass from process catalog
208  const double mass = (*Dep::TH_ProcessCatalog).getParticleProperty(DMid).mass;
209  const double Ecm = 2*mass;
210 
211  // Loop over all channels for that process
212  daFunk::Funk Yield = daFunk::zero("E", "v");
213 
214  // Adding two-body channels
215  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
216  it != annProc.channelList.end(); ++it)
217  {
218  bool added = false; // If spectrum is not available from any source
219 
220  // Here only take care of two-body final states
221  if (it->nFinalStates != 2) continue;
222 
223  // Get final state masses
224  double m0 = (*Dep::TH_ProcessCatalog).getParticleProperty(
225  it->finalStateIDs[0]).mass;
226  double m1 = (*Dep::TH_ProcessCatalog).getParticleProperty(
227  it->finalStateIDs[1]).mass;
228 
229  // Ignore channels that are kinematically closed for v=0
230  if ( m0 + m1 > Ecm ) continue;
231 
232  // Ignore channels with 0 BR in v=0 limit
233  if (it->genRate->bind("v")->eval(0.) <= 0.) continue;
234 
235  double E0 = 0.5*(Ecm*Ecm+m0*m0-m1*m1)/Ecm;
236  double E1 = Ecm-E0;
237 
238  // Check whether two-body final state is in SimYield table
239  if ( Dep::SimYieldTable->hasChannel(
240  it->finalStateIDs[0], it->finalStateIDs[1], "gamma") )
241  {
242  Yield = Yield +
243  it->genRate*(*Dep::SimYieldTable)(
244  it->finalStateIDs[0], it->finalStateIDs[1], "gamma", Ecm);
245  added = true;
246  }
247  // Deal with composite final states
248  else
249  {
250  daFunk::Funk spec0 = daFunk::zero("E");
251  daFunk::Funk spec1 = daFunk::zero("E");
252  added = true;
253 
254  // Final state particle one
255  // Tabulated spectrum available?
256  if ( Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], "gamma") )
257  {
258  spec0 = (*Dep::SimYieldTable)(it->finalStateIDs[0], "gamma")->set("Ecm",E0);
259  }
260  // Gamma-ray line?
261  else if ( it->finalStateIDs[0] == "gamma" )
262  {
263  daFunk::Funk E = daFunk::var("E");
264  spec0 = exp(-pow((E-E0)/line_width/E0,2)/2)/E0/sqrt(2*M_PI)/line_width;
265  }
266  // MC spectra available?
267  else if ( Dep::cascadeMC_gammaSpectra->count(it->finalStateIDs[0]) )
268  {
269  double gamma0 = E0/m0;
270  //std::cout << it->finalStateIDs[0] << " " << gamma0 << std::endl;
271  spec0 = boost_dNdE(Dep::cascadeMC_gammaSpectra->at(it->finalStateIDs[0]), gamma0, 0.0);
272  }
273  else added = false;
274 
275  // Final state particle two
276  if ( Dep::SimYieldTable->hasChannel(it->finalStateIDs[1], "gamma") )
277  {
278  spec1 = (*Dep::SimYieldTable)(it->finalStateIDs[1], "gamma")->set("Ecm", E1);
279  }
280  else if ( it->finalStateIDs[1] == "gamma" )
281  {
282  daFunk::Funk E = daFunk::var("E");
283  spec1 = exp(-pow((E-E1)/line_width/E1,2)/2)/E1/sqrt(2*M_PI)/line_width;
284  }
285  else if ( Dep::cascadeMC_gammaSpectra->count(it->finalStateIDs[1]) )
286  {
287  double gamma1 = E1/m1;
288  //std::cout << it->finalStateIDs[1] << " " << gamma1 << std::endl;
289  spec1 = boost_dNdE(Dep::cascadeMC_gammaSpectra->at(it->finalStateIDs[1]), gamma1, 0.0);
290  }
291  else added = false;
292 
293  #ifdef DARKBIT_DEBUG
294  std::cout << it->finalStateIDs[0] << " " << it->finalStateIDs[1] << std::endl;
295  //std::cout << "gammas: " << gamma0 << ", " << gamma1 << std::endl;
296  daFunk::Funk chnSpec = (daFunk::zero("v", "E")
297  + spec0
298  + spec1)-> set("v", 0.);
299  auto x = daFunk::logspace(0, 3, 10);
300  std::vector<double> y = chnSpec->bind("E")->vect(x);
301  std::cout << it->finalStateIDs[0] << it->finalStateIDs[1] << ":\n";
302  std::cout << " E: [";
303  for (std::vector<double>::iterator it2 = x.begin(); it2 != x.end(); it2++)
304  std::cout << *it2 << ", ";
305  std::cout << "]\n";
306  std::cout << " dNdE: [";
307  for (std::vector<double>::iterator it2 = y.begin(); it2 != y.end(); it2++)
308  std::cout << *it2 << ", ";
309  std::cout << "]\n";
310  #endif
311 
312  if (!added)
313  {
314  DarkBit_warning().raise(LOCAL_INFO,
315  "GA_AnnYield_General: cannot find spectra for "
316  + it->finalStateIDs[0] + " " + it->finalStateIDs[1]);
317  }
318 
319  Yield = Yield + (spec0 + spec1) * it->genRate;
320  }
321  } // End adding two-body final states
322 
323  #ifdef DARKBIT_DEBUG
324  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+");
325  std::vector<std::string> test2 = initVector<std::string> ("A0_test", "A0_test", "Z0_test", "Z0_test", "WH_test", "Z0_test", "h0_2_test", "W-");
326 
327  for(size_t i=0; i<test1.size();i++)
328  {
329  daFunk::Funk chnSpec = (*Dep::SimYieldTable)(test1[i], test2[i], "gamma", Ecm);
330  std::vector<double> y = chnSpec->bind("E")->vect(x);
331  os << test1[i] << test2[i] << ":\n";
332  os << " E: [";
333  for (std::vector<double>::iterator it2 = x.begin(); it2 != x.end(); it2++)
334  os << *it2 << ", ";
335  os << "]\n";
336  os << " dNdE: [";
337  for (std::vector<double>::iterator it2 = y.begin(); it2 != y.end(); it2++)
338  os << *it2 << ", ";
339  os << "]\n";
340  }
341  #endif
342 
343  // Adding three-body final states
344  //
345  // NOTE: Three body processes are added even if they are closed for at v=0
346  for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
347  it != annProc.channelList.end(); ++it)
348  {
349  bool added = true;
350 
351  // Here only take care of three-body final states
352  if (it->nFinalStates != 3) continue;
353 
354  // Implement tabulated three-body final states
355  /*
356  if ( it->nFinalStates == 3
357  and Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], "gamma")
358  and Dep::SimYieldTable->hasChannel(it->finalStateIDs[1], "gamma")
359  and Dep::SimYieldTable->hasChannel(it->finalStateIDs[2], "gamma")
360  )
361  {
362  daFunk::Funk dNdE1dE2 = it->genRate->set("v",0.);
363  daFunk::Funk spec0 =
364  (*Dep::SimYieldTable)(it->finalStateIDs[0], "gamma");
365  daFunk::Funk spec1 =
366  (*Dep::SimYieldTable)(it->finalStateIDs[1], "gamma");
367  daFunk::Funk spec2 =
368  (*Dep::SimYieldTable)(it->finalStateIDs[2], "gamma");
369  Yield = Yield + convspec(spec0, spec1, spec2, dNdE1dE2);
370  }
371  */
372 
373  if ( it->finalStateIDs[0] == "gamma" )
374  {
375  if ( it->finalStateIDs[1] == "gamma" or it->finalStateIDs[2] == "gamma")
376  {
377  DarkBit_warning().raise(LOCAL_INFO, "Second and/or third primary gamma rays in three-body final states ignored.");
378  }
379  double m1 = (*Dep::TH_ProcessCatalog).getParticleProperty(
380  it->finalStateIDs[1]).mass;
381  double m2 = (*Dep::TH_ProcessCatalog).getParticleProperty(
382  it->finalStateIDs[2]).mass;
384  mass, m1, m2);
386  mass, m1, m2);
387  daFunk::Funk dsigmavde = it->genRate->gsl_integration(
388  "E1", E1_low, E1_high);
389 
390  #ifdef DARKBIT_DEBUG
391  daFunk::Funk chnSpec = (daFunk::zero("v", "E") + dsigmavde)-> set("v", 0.);
392  std::vector<double> y = chnSpec->bind("E")->vect(x);
393  os << it->finalStateIDs[0] << it->finalStateIDs[1] << it->finalStateIDs[2] << ":\n";
394  os << " E: [";
395  for (std::vector<double>::iterator it2 = x.begin(); it2 != x.end(); it2++)
396  os << *it2 << ", ";
397  os << "]\n";
398  os << " dNdE: [";
399  for (std::vector<double>::iterator it2 = y.begin(); it2 != y.end(); it2++)
400  os << *it2 << ", ";
401  os << "]\n";
402  #endif
403 
404  Yield = Yield + dsigmavde;
405  }
406  else added = false;
407 
408  if (!added)
409  {
410  DarkBit_warning().raise(LOCAL_INFO,
411  "GA_AnnYield_General: ignoring final state "
412  + it->finalStateIDs[0] + " " + it->finalStateIDs[1] + " " + it->finalStateIDs[2]);
413  }
414  }
415  #ifdef DARKBIT_DEBUG
416  if(debug) os.close();
417  #endif
418 
419  result = k*daFunk::ifelse(1e-6 - daFunk::var("v"), Yield/(mass*mass),
420  daFunk::throwError("Spectrum currently only defined for v=0."));
421  }
Funk throwError(std::string msg)
Definition: daFunk.hpp:1411
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:150
void GA_AnnYield_General(daFunk::Funk &result)
General routine to derive annihilation yield.
Definition: GamYields.cpp:191
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
Funk ifelse(Funk f, Funk g, Funk h)
Definition: daFunk.hpp:1373
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:

◆ 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 
68  if ( runOptions->getValueOrDef(false, "ignore_all") ) return;
69 
70  TH_Process process = (*Dep::TH_ProcessCatalog).getProcess(DMid, DMid);
71 
72  // Add only gamma-ray spectra for two and three body final states
73  for (std::vector<TH_Channel>::iterator it = process.channelList.begin();
74  it != process.channelList.end(); ++it)
75  {
76  if ( it->nFinalStates == 2 )
77  {
79  if ( not runOptions->getValueOrDef(false, "ignore_two_body") )
80  {
81  #ifdef DARKBIT_DEBUG
82  std::cout << "Checking for missing two-body final states: "
83  << it->finalStateIDs[0] << " " << it->finalStateIDs[1] << std::endl;
84  #endif
85  if ( not Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], it->finalStateIDs[1], "gamma") )
86  {
87  if ( not Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], "gamma") )
88  missingFinalStates.insert(it->finalStateIDs[0]);
89  if ( not Dep::SimYieldTable->hasChannel(it->finalStateIDs[1], "gamma") )
90  missingFinalStates.insert(it->finalStateIDs[1]);
91  }
92  }
93  }
94  else if ( it->nFinalStates == 3 )
95  {
97  if ( not runOptions->getValueOrDef(false, "ignore_three_body") )
98  {
99  #ifdef DARKBIT_DEBUG
100  std::cout << "Checking for missing three-body final states: "
101  << it->finalStateIDs[0] << " " << it->finalStateIDs[1]
102  << " " << it->finalStateIDs[2] << std::endl;
103  #endif
104  if (not Dep::SimYieldTable->hasChannel(it->finalStateIDs[0], "gamma") )
105  missingFinalStates.insert(it->finalStateIDs[0]);
106  if (not Dep::SimYieldTable->hasChannel(it->finalStateIDs[1], "gamma") )
107  missingFinalStates.insert(it->finalStateIDs[1]);
108  if (not Dep::SimYieldTable->hasChannel(it->finalStateIDs[2], "gamma") )
109  missingFinalStates.insert(it->finalStateIDs[2]);
110  }
111  }
112  }
113  // Remove particles we don't have decays for.
114  for (auto it = missingFinalStates.begin(); it != missingFinalStates.end();)
115  {
116  if ((*Dep::TH_ProcessCatalog).find(*it, "") == NULL)
117  {
118  #ifdef DARKBIT_DEBUG
119  std::cout << "Erasing (because no decays known): " << *it << std::endl;
120  #endif
121  missingFinalStates.erase(it++);
122  }
123  else
124  {
125  #ifdef DARKBIT_DEBUG
126  std::cout << "Keeping (because decay known): " << *it << std::endl;
127  #endif
128  ++it;
129  }
130  }
131 
132  #ifdef DARKBIT_DEBUG
133  std::cout << "Number of missing final states: " << missingFinalStates.size() << std::endl;
134  for (auto it = missingFinalStates.begin(); it != missingFinalStates.end(); it++)
135  {
136  std::cout << *it << std::endl;
137  }
138  #endif
139 
140  result.assign(missingFinalStates.begin(), missingFinalStates.end());
141  }
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 140 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().

141  {
142  using namespace Pipes::GalacticHalo_Einasto;
143  double rhos = *Param["rhos"];
144  double rs = *Param["rs"];
145  double r_sun = *Param["r_sun"];
146  double alpha = *Param["alpha"];
147  result.DensityProfile = daFunk::func(profile_Einasto, rhos, rs, alpha, daFunk::var("r"));
148  result.r_sun = r_sun;
149  }
double profile_Einasto(double rhos, double rs, double alpha, double r)
Einasto dark matter halo profile function.
Definition: DarkBit.cpp:122
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:140
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 126 of file DarkBit.cpp.

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

127  {
128  using namespace Pipes::GalacticHalo_gNFW;
129  double rhos = *Param["rhos"];
130  double rs = *Param["rs"];
131  double r_sun = *Param["r_sun"];
132  double alpha = *Param["alpha"];
133  double beta = *Param["beta"];
134  double gamma = *Param["gamma"];
135  result.DensityProfile = daFunk::func(profile_gNFW, rhos, rs, alpha, beta, gamma, daFunk::var("r"));
136  result.r_sun = r_sun;
137  }
START_MODEL beta
Definition: Axions.hpp:35
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:126
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:118
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 1048 of file Axions.cpp.

References Gambit::Stats::gaussian_loglikelihood().

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

1048 { 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 199 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().

200  {
201  const SubSpectrum& he = spec.get_HE();
202  double mass = spec.get(Par::Pole_Mass,"S");
203  double lambda = he.get(Par::dimensionless,"lambda_hS");
204  double mh = spec.get(Par::Pole_Mass,"h0_1");
205 
206  // Expressions taken from Cline et al. (2013, PRD 88:055025, arXiv:1306.4710)
207  double fp = 2./9. + 7./9.*(*Param["fpu"] + *Param["fpd"] + *Param["fps"]);
208  double fn = 2./9. + 7./9.*(*Param["fnu"] + *Param["fnd"] + *Param["fns"]);
209 
210  result.gps = lambda*fp*m_proton/pow(mh,2)/mass/2;
211  result.gns = lambda*fn*m_neutron/pow(mh,2)/mass/2;
212  result.gpa = 0; // Only SI cross-section
213  result.gna = 0;
214 
215  logger() << LogTags::debug << "Singlet DM DD couplings:" << std::endl;
216  logger() << " gps = " << result.gps << std::endl;
217  logger() << " gns = " << result.gns << std::endl;
218  logger() << " gpa = " << result.gpa << std::endl;
219  logger() << " gna = " << result.gna << EOM;
220  }
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:99
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 1069 of file RelicDensity.cpp.

References Gambit::byVal().

1070  {
1071  using namespace Pipes::get_semi_ann_MicrOmegas;
1072 
1073  double Beps; // Beps=1e-5 recommended, Beps=1 switches coannihilation off
1074  Beps = runOptions->getValueOrDef<double>(1e-5, "Beps");
1075 
1076  double Xf = *Dep::Xf;
1077 
1078  char*n1 = (char *)"~SS";
1079  char*n2 = (char *)"~SS";
1080  char*n3 = (char *)"h";
1081  char*n4 = (char *)"~ss";
1082 
1083  result = BEreq::get_oneChannel(byVal(Xf),byVal(Beps),byVal(n1),byVal(n2),byVal(n3),byVal(n4));
1084 
1085  }
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:

◆ gStar()

double Gambit::DarkBit::gStar ( double  T)

Various capabilities and functions to provide SM physics as well as QCD input for axions.

Supported models: QCDAxion

Definition at line 1066 of file Axions.cpp.

Referenced by equation_Tosc(), and hubble_rad_dom().

1067  {
1068  // Needs log10(T/GeV) for interpolation.
1069  double lgT = log10(T) - 3.0;
1070  // Interpolated effective relatvistic d.o.f. based on 0910.1066, deviations < 0.5%
1071  // Tabulated data: x = log10(T/GeV), y = gStar
1072  static AxionInterpolator gR (GAMBIT_DIR "/DarkBit/data/gR_WantzShellard.dat", "cspline");
1073  double res;
1074  if (lgT > 3.0) {
1075  res = gR.interpolate (2.99);
1076  } else if (lgT < -5.0) {
1077  res = gR.interpolate (-4.99);
1078  } else {
1079  res = gR.interpolate (lgT);
1080  }
1081 
1082  return res;
1083  }
Here is the caller graph for this function:

◆ gStar_S()

double Gambit::DarkBit::gStar_S ( double  T)

Definition at line 1086 of file Axions.cpp.

Referenced by RD_oh2_Axions().

1087  {
1088  // Need log10(T/GeV) for interpolation.
1089  double lgT = log10(T) - 3.0;
1090  // Interpolated effective relatvistic d.o.f. based on 0910.1066, deviations < 0.5%
1091  // Tabulated data: x = log10(T/GeV), y = gStar
1092  static AxionInterpolator gS (GAMBIT_DIR "/DarkBit/data/gS_WantzShellard.dat", "cspline");
1093  double res;
1094  if (lgT > 3.0) {
1095  res = gS.interpolate (2.99);
1096  } else if (lgT < -5.0) {
1097  res = gS.interpolate (-4.99);
1098  } else {
1099  res = gS.interpolate (lgT);
1100  }
1101 
1102  return res;
1103  }
Here is the caller graph for this function:

◆ hubble_rad_dom()

double Gambit::DarkBit::hubble_rad_dom ( double  T)

Definition at line 1630 of file Axions.cpp.

References gStar(), and Gambit::m_planck_red.

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

1631  {
1632  // H(T)/eV, T/MeV, m_pl/10^12eV = m_pl/10^3 GeV
1633  const double m_pl = m_planck_red*1.0E-3;
1634  return 0.331153*sqrt(gStar(T))*T*T/m_pl;
1635  }
double gStar(double T)
Various capabilities and functions to provide SM physics as well as QCD input for axions...
Definition: Axions.cpp:1066
const double m_planck_red
Here is the call graph for this function:
Here is the caller graph for this function:

◆ IC22_bg()

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

Definition at line 616 of file SunNeutrinos.cpp.

616  {
617  result = Pipes::IC22_bg ::Dep::IC22_data->bg; }

◆ IC22_bgloglike()

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

Definition at line 622 of file SunNeutrinos.cpp.

Referenced by IC_loglike().

622  {
623  result = Pipes::IC22_bgloglike::Dep::IC22_data->bgloglike; }
Here is the caller graph for this function:

◆ IC22_full()

void Gambit::DarkBit::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.

We define the final log-likelihood as delta = sum over analyses of (lnL_model - lnL_BG), conservatively forbidding delta > 0 in order to always just use the neutrino likelihood as a limit. This ignores small low-E excesses caused by impending breakdown of approximations used in IceCube response data and the nulike likelihood at very low E. This implies conditioning on all but one parameter (e.g. the cross-section), such that including any particular IC analysis adds just one additional degree of freedom to the fit.22-string IceCube sample: predicted signal and background counts, observed counts and likelihoods.

Option nulike_speed<int>: Speed setting for nulike backend (default 3)

Definition at line 485 of file SunNeutrinos.cpp.

References annihilation_rate_Sun(), Gambit::DarkBit::nudata::bg, Gambit::DarkBit::nudata::bgloglike, Gambit::byVal(), Gambit::Piped_invalid_point::check(), Gambit::Piped_exceptions::check(), DarkBit_error(), DarkBit_warning(), Gambit::DarkBit::nudata::loglike, Gambit::DarkBit::nudata::nobs, Gambit::piped_errors, Gambit::piped_invalid_point, Gambit::piped_warnings, Gambit::DarkBit::nudata::pvalue, and Gambit::DarkBit::nudata::signal.

486  {
487  using namespace Pipes::IC22_full;
488  static bool first = true;
489  const double bgloglike = -808.4581;
490  double sigpred, bgpred, lnLike, pval;
491  int totobs;
492  char experiment[300] = "IC-22";
493  void* context = NULL;
494  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
496  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
497  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
498  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
499  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
503  result.signal = sigpred;
504  result.bg = bgpred;
505  result.nobs = totobs;
506  result.loglike = lnLike;
507  result.pvalue = pval;
508  if (first)
509  {
510  result.bgloglike = bgloglike;
511  first = false;
512  }
513  }
error & DarkBit_error()
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
void annihilation_rate_Sun(double &result)
Annihilation rate of dark matter in the Sun (s^-1)
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
void IC22_full(nudata &result)
Likelihood calculators for different IceCube event samples These functions all include the likelihood...
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:544
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:473
warning & DarkBit_warning()
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:

◆ IC22_loglike()

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

Definition at line 620 of file SunNeutrinos.cpp.

Referenced by IC_loglike().

620  {
621  result = Pipes::IC22_loglike::Dep::IC22_data->loglike; }
Here is the caller graph for this function:

◆ IC22_nobs()

void Gambit::DarkBit::IC22_nobs ( int result)

Definition at line 618 of file SunNeutrinos.cpp.

618  {
619  result = Pipes::IC22_nobs ::Dep::IC22_data->nobs; }

◆ IC22_pvalue()

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

Definition at line 624 of file SunNeutrinos.cpp.

624  {
625  result = Pipes::IC22_pvalue ::Dep::IC22_data->pvalue; }

◆ IC22_signal()

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

22-string extractor module functions

Definition at line 614 of file SunNeutrinos.cpp.

614  {
615  result = Pipes::IC22_signal ::Dep::IC22_data->signal; }

◆ IC79_loglike()

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

Composite IceCube 79-string likelihood function.

Definition at line 677 of file SunNeutrinos.cpp.

References IC79SL_bgloglike(), IC79SL_loglike(), IC79WH_bgloglike(), IC79WH_loglike(), IC79WL_bgloglike(), and IC79WL_loglike().

Referenced by main().

678  {
679  using namespace Pipes::IC79_loglike;
683  if (result > 0.0) result = 0.0;
684  }
void IC79SL_bgloglike(double &result)
void IC79WL_bgloglike(double &result)
void IC79WL_loglike(double &result)
void IC79WH_loglike(double &result)
void IC79SL_loglike(double &result)
void IC79_loglike(double &result)
Composite IceCube 79-string likelihood function.
void IC79WH_bgloglike(double &result)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ IC79SL_bg()

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

Definition at line 664 of file SunNeutrinos.cpp.

664  {
665  result = Pipes::IC79SL_bg ::Dep::IC79SL_data->bg; }

◆ IC79SL_bgloglike()

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

Definition at line 670 of file SunNeutrinos.cpp.

Referenced by IC79_loglike(), IC_loglike(), and main().

670  {
671  result = Pipes::IC79SL_bgloglike::Dep::IC79SL_data->bgloglike; }
Here is the caller graph for this function:

◆ IC79SL_full()

void Gambit::DarkBit::IC79SL_full ( nudata result)

79-string IceCube SL sample: predicted signal and background counts, observed counts and likelihoods.

Option nulike_speed<int>: Speed setting for nulike backend (default 3)

Definition at line 581 of file SunNeutrinos.cpp.

References annihilation_rate_Sun(), Gambit::DarkBit::nudata::bg, Gambit::DarkBit::nudata::bgloglike, Gambit::byVal(), Gambit::Piped_invalid_point::check(), Gambit::Piped_exceptions::check(), DarkBit_error(), DarkBit_warning(), Gambit::DarkBit::nudata::loglike, Gambit::DarkBit::nudata::nobs, Gambit::piped_errors, Gambit::piped_invalid_point, Gambit::piped_warnings, Gambit::DarkBit::nudata::pvalue, and Gambit::DarkBit::nudata::signal.

Referenced by main().

582  {
583  static bool first = true;
584  using namespace Pipes::IC79SL_full;
585  const double bgloglike = -5015.6474;
586  double sigpred, bgpred, lnLike, pval;
587  int totobs;
588  char experiment[300] = "IC-79 SL";
589  void* context = NULL;
590  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
592  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
593  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
594  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
595  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
599  result.signal = sigpred;
600  result.bg = bgpred;
601  result.nobs = totobs;
602  result.loglike = lnLike;
603  result.pvalue = pval;
604  if (first)
605  {
606  result.bgloglike = bgloglike;
607  first = false;
608  }
609  }
error & DarkBit_error()
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
void annihilation_rate_Sun(double &result)
Annihilation rate of dark matter in the Sun (s^-1)
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
void IC79SL_full(nudata &result)
79-string IceCube SL sample: predicted signal and background counts, observed counts and likelihoods...
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:544
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:473
warning & DarkBit_warning()
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:

◆ IC79SL_loglike()

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

Definition at line 668 of file SunNeutrinos.cpp.

Referenced by IC79_loglike(), IC_loglike(), and main().

668  {
669  result = Pipes::IC79SL_loglike::Dep::IC79SL_data->loglike; }
Here is the caller graph for this function:

◆ IC79SL_nobs()

void Gambit::DarkBit::IC79SL_nobs ( int result)

Definition at line 666 of file SunNeutrinos.cpp.

666  {
667  result = Pipes::IC79SL_nobs ::Dep::IC79SL_data->nobs; }

◆ IC79SL_pvalue()

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

Definition at line 672 of file SunNeutrinos.cpp.

672  {
673  result = Pipes::IC79SL_pvalue ::Dep::IC79SL_data->pvalue; }

◆ IC79SL_signal()

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

79-string SL extractor module functions

Definition at line 662 of file SunNeutrinos.cpp.

662  {
663  result = Pipes::IC79SL_signal ::Dep::IC79SL_data->signal; }

◆ IC79WH_bg()

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

Definition at line 632 of file SunNeutrinos.cpp.

632  {
633  result = Pipes::IC79WH_bg ::Dep::IC79WH_data->bg; }

◆ IC79WH_bgloglike()

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

Definition at line 638 of file SunNeutrinos.cpp.

Referenced by IC79_loglike(), IC_loglike(), and main().

638  {
639  result = Pipes::IC79WH_bgloglike::Dep::IC79WH_data->bgloglike; }
Here is the caller graph for this function:

◆ IC79WH_full()

void Gambit::DarkBit::IC79WH_full ( nudata result)

79-string IceCube WH sample: predicted signal and background counts, observed counts and likelihoods.

Option nulike_speed<int>: Speed setting for nulike backend (default 3)

Definition at line 517 of file SunNeutrinos.cpp.

References annihilation_rate_Sun(), Gambit::DarkBit::nudata::bg, Gambit::DarkBit::nudata::bgloglike, Gambit::byVal(), Gambit::Piped_invalid_point::check(), Gambit::Piped_exceptions::check(), DarkBit_error(), DarkBit_warning(), Gambit::DarkBit::nudata::loglike, Gambit::DarkBit::nudata::nobs, Gambit::piped_errors, Gambit::piped_invalid_point, Gambit::piped_warnings, Gambit::DarkBit::nudata::pvalue, and Gambit::DarkBit::nudata::signal.

Referenced by main().

518  {
519  static bool first = true;
520  using namespace Pipes::IC79WH_full;
521  const double bgloglike = -11874.8689;
522  double sigpred, bgpred, lnLike, pval;
523  int totobs;
524  char experiment[300] = "IC-79 WH";
525  void* context = NULL;
526  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
528  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
529  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
530  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
531  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
535  result.signal = sigpred;
536  result.bg = bgpred;
537  result.nobs = totobs;
538  result.loglike = lnLike;
539  result.pvalue = pval;
540  if (first)
541  {
542  result.bgloglike = bgloglike;
543  first = false;
544  }
545  }
error & DarkBit_error()
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
void annihilation_rate_Sun(double &result)
Annihilation rate of dark matter in the Sun (s^-1)
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
void IC79WH_full(nudata &result)
79-string IceCube WH sample: predicted signal and background counts, observed counts and likelihoods...
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:544
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:473
warning & DarkBit_warning()
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:

◆ IC79WH_loglike()

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

Definition at line 636 of file SunNeutrinos.cpp.

Referenced by IC79_loglike(), IC_loglike(), and main().

636  {
637  result = Pipes::IC79WH_loglike::Dep::IC79WH_data->loglike; }
Here is the caller graph for this function:

◆ IC79WH_nobs()

void Gambit::DarkBit::IC79WH_nobs ( int result)

Definition at line 634 of file SunNeutrinos.cpp.

634  {
635  result = Pipes::IC79WH_nobs ::Dep::IC79WH_data->nobs; }

◆ IC79WH_pvalue()

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

Definition at line 640 of file SunNeutrinos.cpp.

640  {
641  result = Pipes::IC79WH_pvalue ::Dep::IC79WH_data->pvalue; }

◆ IC79WH_signal()

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

79-string WH extractor module functions

Definition at line 630 of file SunNeutrinos.cpp.

630  {
631  result = Pipes::IC79WH_signal ::Dep::IC79WH_data->signal; }

◆ IC79WL_bg()

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

Definition at line 648 of file SunNeutrinos.cpp.

648  {
649  result = Pipes::IC79WL_bg ::Dep::IC79WL_data->bg; }

◆ IC79WL_bgloglike()

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

Definition at line 654 of file SunNeutrinos.cpp.

Referenced by IC79_loglike(), IC_loglike(), and main().

654  {
655  result = Pipes::IC79WL_bgloglike::Dep::IC79WL_data->bgloglike; }
Here is the caller graph for this function:

◆ IC79WL_full()

void Gambit::DarkBit::IC79WL_full ( nudata result)

79-string IceCube WL sample: predicted signal and background counts, observed counts and likelihoods.

Option nulike_speed<int>: Speed setting for nulike backend (default 3)

Definition at line 549 of file SunNeutrinos.cpp.

References annihilation_rate_Sun(), Gambit::DarkBit::nudata::bg, Gambit::DarkBit::nudata::bgloglike, Gambit::byVal(), Gambit::Piped_invalid_point::check(), Gambit::Piped_exceptions::check(), DarkBit_error(), DarkBit_warning(), Gambit::DarkBit::nudata::loglike, Gambit::DarkBit::nudata::nobs, Gambit::piped_errors, Gambit::piped_invalid_point, Gambit::piped_warnings, Gambit::DarkBit::nudata::pvalue, and Gambit::DarkBit::nudata::signal.

Referenced by main().

550  {
551  static bool first = true;
552  using namespace Pipes::IC79WL_full;
553  const double bgloglike = -1813.4503;
554  double sigpred, bgpred, lnLike, pval;
555  int totobs;
556  char experiment[300] = "IC-79 WL";
557  void* context = NULL;
558  double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
560  int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
561  BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
562  byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
563  theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
567  result.signal = sigpred;
568  result.bg = bgpred;
569  result.nobs = totobs;
570  result.loglike = lnLike;
571  result.pvalue = pval;
572  if (first)
573  {
574  result.bgloglike = bgloglike;
575  first = false;
576  }
577  }
error & DarkBit_error()
Piped_exceptions piped_errors
Global instance of Piped_exceptions class for errors.
void IC79WL_full(nudata &result)
79-string IceCube WL sample: predicted signal and background counts, observed counts and likelihoods...
Piped_exceptions piped_warnings
Global instance of Piped_exceptions class for warnings.
void annihilation_rate_Sun(double &result)
Annihilation rate of dark matter in the Sun (s^-1)
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:524
void check(exception &excep)
Check whether any exceptions were requested, and raise them.
Definition: exceptions.cpp:544
void check()
Check whether an exception was requested, and throw it if necessary.
Definition: exceptions.cpp:473
warning & DarkBit_warning()
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:

◆ IC79WL_loglike()

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

Definition at line 652 of file SunNeutrinos.cpp.

Referenced by IC79_loglike(), IC_loglike(), and main().

652  {
653  result = Pipes::IC79WL_loglike::Dep::IC79WL_data->loglike; }
Here is the caller graph for this function:

◆ IC79WL_nobs()

void Gambit::DarkBit::IC79WL_nobs ( int result)

Definition at line 650 of file SunNeutrinos.cpp.

650  {
651  result = Pipes::IC79WL_nobs ::Dep::IC79WL_data->nobs; }

◆ IC79WL_pvalue()

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

Definition at line 656 of file SunNeutrinos.cpp.

656  {
657  result = Pipes::IC79WL_pvalue ::Dep::IC79WL_data->pvalue; }

◆ IC79WL_signal()

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

79-string WL extractor module functions

Definition at line 646 of file SunNeutrinos.cpp.

646  {
647  result = Pipes::IC79WL_signal ::Dep::IC79WL_data->signal; }

◆ IC_loglike()

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

Complete composite IceCube likelihood function.

Definition at line 687 of file SunNeutrinos.cpp.

References IC22_bgloglike(), IC22_loglike(), IC79SL_bgloglike(), IC79SL_loglike(), IC79WH_bgloglike(), IC79WH_loglike(), IC79WL_bgloglike(), and IC79WL_loglike().

688  {
689  using namespace Pipes::IC_loglike;
694  if (result > 0.0) result = 0.0;
695 #ifdef DARKBIT_DEBUG
696  cout << "IC likelihood: " << result << endl;
697  cout << "IC79SL contribution: " << *Dep::IC79SL_loglike - *Dep::IC79SL_bgloglike << endl;
698  cout << "IC79WL contribution: " << *Dep::IC79WL_loglike - *Dep::IC79WL_bgloglike << endl;
699  cout << "IC79WH contribution: " << *Dep::IC79WH_loglike - *Dep::IC79WH_bgloglike << endl;
700  cout << "IC22 contribution: " << *Dep::IC22_loglike - *Dep::IC22_bgloglike << endl;
701 #endif
702  }
void IC22_bgloglike(double &result)
void IC79SL_bgloglike(double &result)
void IC79WL_bgloglike(double &result)
void IC79WL_loglike(double &result)
void IC79WH_loglike(double &result)
void IC79SL_loglike(double &result)
void IC22_loglike(double &result)
void IC_loglike(double &result)
Complete composite IceCube likelihood function.
void IC79WH_bgloglike(double &result)
Here is the call graph for this function:

◆ intersect_parabola_line()

double Gambit::DarkBit::intersect_parabola_line ( double  a,
double  b,
double  sign,
const double  pparams[] 
)

Definition at line 345 of file Axions.cpp.

References Gambit::alpha_EM, Gambit::atomic_mass_unit, b, DarkBit_error(), Gambit::GreAT::data, Gambit::LogTags::debug, Gambit::FlavBit::LoopFunctions::E(), Gambit::EOM, Gambit::eV2g, Gambit::gev2cm, Gambit::LogTags::info, daFunk::interp(), Gambit::K2eV, LOCAL_INFO, Gambit::logger(), Gambit::m_electron, parabola(), Gambit::pi, r, and Gambit::swap().

346  {
347  const double x1 = -3.673776;
348  const double y1 = 0.4;
349  double x2 = a - pparams[1];
350  double temp1 = a*a + 4.0*b*pparams[0] - 2.0*a*pparams[1] + pparams[1]*pparams[1] - 4.0*pparams[0]*pparams[2];
351  x2 = x2 - sign*sqrt(temp1);
352  x2 = x2/(2.0*pparams[0]);
353  double y2 = parabola(x2, pparams);
354  temp1 = x1 - x2;
355  double temp2 = y1 - y2;
356  return sqrt(temp1*temp1 + temp2*temp2);
357  }
START_MODEL b
Definition: demo.hpp:235
double parabola(double x, const double params[])
H.E.S.S.-likelihood-related interpolation routines.
Definition: Axions.cpp:342
Here is the call graph for this function:

◆ lnL_CTAGC_gamLike()

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

Definition at line 189 of file SimpleLikelihoods.cpp.

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

Referenced by main().

190  {
191  using namespace Pipes::lnL_CTAGC_gamLike;
192 
193  double fraction = *Dep::RD_fraction;
194  result = 0;
195 
196  // from 25 GeV to 10 TeV
197  std::vector<double> x = daFunk::logspace(1.39, 4.00, 100);
198  x = daFunk::augmentSingl(x, (*Dep::GA_AnnYield)->set("v",0));
199  std::vector<double> y = ((*Dep::GA_AnnYield)/8./M_PI*fraction*fraction)->
200  set("v", 0)->bind("E")->vect(x);
201 
202  result = BEreq::lnL(5, x, y);
203 
204  logger() << LogTags::debug << "GamLike CTA GC likelihood is lnL = " << result << EOM;
205  }
std::vector< double > logspace(double x0, double x1, unsigned int n)
Definition: daFunk.hpp:186
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
std::vector< double > augmentSingl(const std::vector< double > &xgrid, Funk f, int N=100, double sigma=3)
Definition: daFunk.hpp:1714
void lnL_CTAGC_gamLike(double &result)
Here is the call graph for this function: