gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
Axions.cpp
Go to the documentation of this file.
1 
21 #include <algorithm>
22 #include <cmath>
23 #include <vector>
24 #include <string>
25 #include <iostream>
26 #include <sstream>
27 
28 #include <gsl/gsl_math.h>
29 #include <gsl/gsl_sf.h>
30 #include <gsl/gsl_sf_trig.h>
31 #include <gsl/gsl_spline.h>
32 #include <gsl/gsl_interp2d.h>
33 #include <gsl/gsl_spline2d.h>
34 #include <gsl/gsl_histogram.h>
35 #include <gsl/gsl_roots.h>
36 #include <gsl/gsl_matrix.h>
37 #include <gsl/gsl_errno.h>
38 #include <gsl/gsl_odeiv2.h>
39 
47 
48 //#define AXION_DEBUG_MODE
49 //#define AXION_OMP_DEBUG_MODE
50 
51 namespace Gambit
52 {
53  namespace DarkBit
54  {
56  // //
57  // General Functions and Classes for Axions //
58  // //
60 
64  // Auxillary functions and classes for interpolation //
67 
71  // AxionInterpolator class: Provides a general 1-D interpolation container based on the gsl library.
72  // Can be declared static for efficiency & easy one-time initialisation of interpolating functions.
73  class AxionInterpolator
74  {
75  public:
76  // Overloaded class creators for the AxionInterpolator class using the init function below.
77  AxionInterpolator();
78  AxionInterpolator(const std::vector<double> x, const std::vector<double> y, std::string type);
79  AxionInterpolator(const std::vector<double> x, const std::vector<double> y);
80  AxionInterpolator(std::string file, std::string type);
81  AxionInterpolator(std::string file);
82  AxionInterpolator& operator=(AxionInterpolator&&);
83  // Destructor.
84  ~AxionInterpolator();
85  // Delete copy constructor and assignment operator to avoid shallow copies.
86  AxionInterpolator(const AxionInterpolator&) = delete;
87  AxionInterpolator operator=(const AxionInterpolator&) = delete;
88  // Routine to access interpolated values.
89  double interpolate(double x);
90  // Routine to access upper and lower boundaries of available data.
91  double lower();
92  double upper();
93  private:
94  // Initialiser for the AxionInterpolator class.
95  void init(std::string file, std::string type);
96  void init(const std::vector<double> x, const std::vector<double> y, std::string type);
97  // The gsl objects for the interpolating functions.
98  gsl_interp_accel *acc;
99  gsl_spline *spline;
100  // Upper and lower boundaries available for the interpolating function.
101  double lo;
102  double up;
103  };
104 
105  // Default constructor.
106  AxionInterpolator::AxionInterpolator() {}
107 
108  // Initialiser for the AxionInterpolator class.
109  void AxionInterpolator::init(const std::vector<double> x, const std::vector<double> y, std::string type)
110  {
111  int pts = x.size();
112  // Get first and last value of the "x" component.
113  lo = x.front();
114  up = x.back();
115  acc = gsl_interp_accel_alloc ();
116  if (type == "cspline")
117  {
118  spline = gsl_spline_alloc (gsl_interp_cspline, pts);
119  }
120  else if (type == "linear")
121  {
122  spline = gsl_spline_alloc (gsl_interp_linear, pts);
123  }
124  else
125  {
126  DarkBit_error().raise(LOCAL_INFO, "ERROR! Interpolation type '"+type+"' not known to class AxionInterpolator.\n Available types: 'linear' and 'cspline'.");
127  }
128 
129  gsl_spline_init (spline, &x[0], &y[0], pts);
130  }
131 
132  // Overloaded class creators for the AxionInterpolator class using the init function above.
133  AxionInterpolator::AxionInterpolator(const std::vector<double> x, const std::vector<double> y, std::string type) { init(x, y, type); }
134  AxionInterpolator::AxionInterpolator(const std::vector<double> x, const std::vector<double> y) { init(x, y, "linear"); }
135 
136  // Initialiser for the AxionInterpolator class.
137  void AxionInterpolator::init(std::string file, std::string type)
138  {
139  // Check if file exists.
140  if (not(Utils::file_exists(file)))
141  {
142  DarkBit_error().raise(LOCAL_INFO, "ERROR! File '"+file+"' not found!");
143  } else {
144  logger() << LogTags::debug << "Reading data from file '"+file+"' and interpolating it with '"+type+"' method." << EOM;
145  }
146  // Read numerical values from data file.
147  ASCIItableReader tab (file);
148  tab.setcolnames("x", "y");
149 
150  init(tab["x"],tab["y"],type);
151  }
152 
153  // Overloaded class creators for the AxionInterpolator class using the init function above.
154  AxionInterpolator::AxionInterpolator(std::string file, std::string type) { init(file, type); }
155  AxionInterpolator::AxionInterpolator(std::string file) { init(file, "linear"); }
156 
157  // Move assignment operator
158  AxionInterpolator& AxionInterpolator::operator=(AxionInterpolator&& interp)
159  {
160  if(this != &interp)
161  {
162  std::swap(acc,interp.acc);
163  std::swap(spline,interp.spline);
164  std::swap(lo,interp.lo);
165  std::swap(up,interp.up);
166  }
167  return *this;
168  }
169 
170  // Destructor
171  AxionInterpolator::~AxionInterpolator()
172  {
173  gsl_spline_free (spline);
174  gsl_interp_accel_free (acc);
175  }
176 
177  // Routine to access interpolated values.
178  double AxionInterpolator::interpolate(double x) { return gsl_spline_eval(spline, x, acc); }
179 
180  // Routines to return upper and lower boundaries of interpolating function
181  double AxionInterpolator::lower() { return lo; }
182  double AxionInterpolator::upper() { return up; }
183 
184 
188  // AxionInterpolator2D class: Provides a 2-D interpolation container based on the gsl library.
189  // Can be declared static for efficiency & easy one-time initialisation of interpolating functions.
190  class AxionInterpolator2D
191  {
192  public:
193  // Overloaded class creators for the AxionInterpolator class using the init function below.
194  AxionInterpolator2D();
195  AxionInterpolator2D(std::string file, std::string type);
196  AxionInterpolator2D(std::string file);
197  AxionInterpolator2D& operator=(AxionInterpolator2D&&);
198  // Destructor.
199  ~AxionInterpolator2D();
200  // Delete copy constructor and assignment operator to avoid shallow copies.
201  AxionInterpolator2D(const AxionInterpolator2D&) = delete;
202  AxionInterpolator2D operator=(const AxionInterpolator2D&) = delete;
203  // Routine to access interpolated values.
204  double interpolate(double x, double y);
205  // Routine to check if a point is inside the interpolating box.
206  bool is_inside_box(double x, double y);
207  private:
208  // Initialiser for the AxionInterpolator2D class.
209  void init(std::string file, std::string type);
210  // The gsl objects for the interpolating functions that need to be available to the class routines.
211  gsl_interp_accel *x_acc;
212  gsl_interp_accel *y_acc;
213  gsl_spline2d *spline;
214  double* z;
215  // Upper and lower "x" and "y" values available to the interpolating function.
216  double x_lo, y_lo, x_up, y_up;
217  };
218 
219  // Move assignment operator
220  AxionInterpolator2D& AxionInterpolator2D::operator=(AxionInterpolator2D&& interp)
221  {
222  if(this != &interp)
223  {
224  std::swap(x_acc,interp.x_acc);
225  std::swap(y_acc,interp.y_acc);
226  std::swap(spline,interp.spline);
227  std::swap(x_lo,interp.x_lo);
228  std::swap(x_up,interp.x_up);
229  std::swap(y_lo,interp.y_lo);
230  std::swap(y_up,interp.y_up);
231  std::swap(z,interp.z);
232  }
233  return *this;
234  }
235 
236  // Destructor
237  AxionInterpolator2D::~AxionInterpolator2D()
238  {
239  gsl_spline2d_free (spline);
240  gsl_interp_accel_free (x_acc);
241  gsl_interp_accel_free (y_acc);
242  free(z);
243  }
244 
245  // Initialiser for the AxionInterpolator class.
246  void AxionInterpolator2D::init(std::string file, std::string type)
247  {
248  // Check if file exists.
249  if (not(Utils::file_exists(file)))
250  {
251  DarkBit_error().raise(LOCAL_INFO, "ERROR! File '"+file+"' not found!");
252  } else {
253  logger() << LogTags::debug << "Reading data from file '"+file+"' and interpolating it with '"+type+"' method." << EOM;
254  }
255  // Read numerical values from data file.
256  ASCIItableReader tab (file);
257  tab.setcolnames("x", "y", "z");
258  // Initialise gsl interpolation routine.
259  // Get unique entries of "x" and "y" for the grid and grid size.
260  std::vector<double> x_vec = tab["x"];
261  sort(x_vec.begin(), x_vec.end());
262  x_vec.erase(unique(x_vec.begin(), x_vec.end()), x_vec.end());
263  int nx = x_vec.size();
264  std::vector<double> y_vec = tab["y"];
265  sort(y_vec.begin(), y_vec.end());
266  y_vec.erase(unique(y_vec.begin(), y_vec.end()), y_vec.end());
267  int ny = y_vec.size();
268  int n_grid_pts = tab["z"].size();
269 
270  if (nx*ny != n_grid_pts)
271  {
272  DarkBit_error().raise(LOCAL_INFO, "ERROR! The number of grid points ("+std::to_string(n_grid_pts)+") for AxionInterpolator2D does not equal the number of unique 'x' and 'y' values ("+std::to_string(nx)+" and "+std::to_string(ny)+")!\n Check formatting of the file: '"+file+"'.");
273  }
274 
275  const double* x = &x_vec[0];
276  const double* y = &y_vec[0];
277  // Allocate memory for "z" values array in gsl format
278  z = (double*) malloc(nx * ny * sizeof(double));
279 
280  if (type == "bicubic")
281  {
282  spline = gsl_spline2d_alloc(gsl_interp2d_bicubic, nx, ny);
283  }
284  else if (type == "bilinear")
285  {
286  spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, nx, ny);
287  }
288  else
289  {
290  DarkBit_error().raise(LOCAL_INFO, "ERROR! Interpolation type '"+type+"' not known to class AxionInterpolator2D.\n Available types: 'bilinear' and 'bicubic'.");
291  }
292 
293  x_acc = gsl_interp_accel_alloc();
294  y_acc = gsl_interp_accel_alloc();
295 
296  // Determine first and last "x" and "y" values and grid step size.
297  x_lo = x_vec.front();
298  x_up = x_vec.back();
299  y_lo = y_vec.front();
300  y_up = y_vec.back();
301  double x_delta = (x_up-x_lo) / (nx-1);
302  double y_delta = (y_up-y_lo) / (ny-1);
303 
304  // Intialise grid.
305  for (int i = 0; i < n_grid_pts; i++)
306  {
307  // Determine appropriate indices for the grid points.
308  double temp = (tab["x"][i]-x_lo) / x_delta;
309  int ind_x = (int) (temp+0.5);
310  temp = (tab["y"][i]-y_lo) / y_delta;
311  int ind_y = (int) (temp+0.5);
312 
313  //std::cout << ind_x << "/" << nx-1 << " " << tab["x"][i] << " vs " << x[ind_x] << " " << ind_y << "/" << ny-1 << " " << tab["y"][i] << " vs " << y[ind_y] << std::endl;
314 
315  gsl_spline2d_set(spline, z, ind_x, ind_y, tab["z"][i]);
316  }
317  gsl_spline2d_init (spline, x, y, z, nx, ny);
318  }
319 
320  // Default creator with dummy entries for the objects w/ memory allocation
321  AxionInterpolator2D::AxionInterpolator2D()
322  {
323  x_acc = gsl_interp_accel_alloc();
324  y_acc = gsl_interp_accel_alloc();
325  spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, 2, 2);
326  z = (double*) malloc(2 * 2 * sizeof(double));
327  }
328  // Overloaded class creators for the AxionInterpolator class using the init function above.
329  AxionInterpolator2D::AxionInterpolator2D(std::string file, std::string type) { init(file, type); }
330  AxionInterpolator2D::AxionInterpolator2D(std::string file) { init(file, "bilinear"); }
331 
332  // Routine to access interpolated values.
333  double AxionInterpolator2D::interpolate(double x, double y) { return gsl_spline2d_eval(spline, x, y, x_acc, y_acc); }
334 
335  // Routine to check if a point is inside the interpolating box.
336  bool AxionInterpolator2D::is_inside_box(double x, double y) { return ((x >= x_lo) && (x <= x_up) && (y >= y_lo) && (y <= y_up)); }
337 
341  // Auxillary function for a parabola (needed for H.E.S.S. likelihood approximation).
342  double parabola(double x, const double params[]) { return params[0]*x*x + params[1]*x + params[2]; }
343 
344  // Auxillary function to return the appropriate intersection between a parabola and a line (needed for H.E.S.S. likelihood).
345  double intersect_parabola_line(double a, double b, double sign, const double pparams[])
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  }
358 
359  // HESS_Interpolator class: Provides a customised interpolation container for the H.E.S.S. likelihood.
360  class HESS_Interpolator
361  {
362  public:
363  // Class creator.
364  HESS_Interpolator(std::string file);
365  // Class destructor
366  ~HESS_Interpolator();
367  // Delete copy constructor and assignment operator to avoid shallow copies
368  HESS_Interpolator(const HESS_Interpolator&) = delete;
369  HESS_Interpolator operator=(const HESS_Interpolator&) = delete;
370  // Container for the tabulated data.
371  ASCIItableReader interp_lnL;
372  // Routine to return interpolated log-likelihood values.
373  double lnL(double epsilon, double gamma);
374  private:
375  gsl_interp_accel *acc[17];
376  gsl_spline *spline[17];
377  };
378 
379  // Class creator. Needs path to tabulated H.E.S.S. data.
380  HESS_Interpolator::HESS_Interpolator(std::string file)
381  {
382  // Initialise upper part of the likelihood interpolation (i.e. higher axion-photon coupling).
383  interp_lnL = ASCIItableReader(file);
384  interp_lnL.setcolnames("lnL16", "lnL15", "lnL14", "lnL13", "lnL12", "lnL11", "lnL10", "lnL9", "lnL8", "lnL7", "lnL6", "lnL5", "lnL4", "lnL3", "lnL2", "lnL1", "lnL0");
385  for (int i = 16; i >= 0; i--)
386  {
387  int pts = interp_lnL["lnL"+std::to_string(i)].size();
388  acc[i] = gsl_interp_accel_alloc ();
389  spline[i] = gsl_spline_alloc (gsl_interp_cspline, pts);
390  const double* epsvals = &interp_lnL["lnL"+std::to_string(i)][0];
391  if (pts==8) {
392  const double lnLvals [8] = {0., -2.30259, -2.99573, -4.60517, -4.60517, -2.99573, -2.30259, 0.};
393  gsl_spline_init (spline[i], epsvals, lnLvals, pts);
394  } else {
395  const double lnLvals [7] = {0., -2.30259, -2.99573, -4.60517, -2.99573, -2.30259, 0.};
396  gsl_spline_init (spline[i], epsvals, lnLvals, pts);
397  }
398  }
399  }
400 
401  // Destructor
402  HESS_Interpolator::~HESS_Interpolator()
403  {
404  for (auto spl : spline)
405  gsl_spline_free (spl);
406  for (auto ac : acc)
407  gsl_interp_accel_free (ac);
408  }
409 
410  // Rotuine to interpolate the H.E.S.S. log-likelihood values.
411  double HESS_Interpolator::lnL(double epsilon, double gamma)
412  {
413  // Parameters for the parabolae.
414  const double ppars00 [3] = {0.553040458173831, 3.9888540782199913, 6.9972958867687565};
415  const double ppars90 [3] = {1.2852894785722664, 9.42311266504736, 17.49643049277964};
416  const double ppars95 [3] = {1.4501115909461886, 10.647792383304218, 19.811978366687622};
417  double result = 0.0;
418 
419  // Check if we are inside the constrained region.
420  if ((gamma > parabola(epsilon, ppars00)) && (gamma < 1.2) && (epsilon > -4.64) && (epsilon < -2.57))
421  {
422  // Check if we are in the upper part (higher coupling; interpolation using linear and cubic splines).
423  if (gamma > 0.4)
424  {
425  // Cubic interpolation in Epsilon.
426  int index_lo = floor((gamma-0.4)/0.05);
427  int index_hi = index_lo + 1;
428  double z_lo = 0.0, z_hi = 0.0;
429  // Only use interpolating function where needed.
430  if ( (epsilon > interp_lnL["lnL"+std::to_string(index_lo)].front()) && (epsilon < interp_lnL["lnL"+std::to_string(index_lo)].back()) )
431  {
432  z_lo = gsl_spline_eval (spline[index_lo], epsilon, acc[index_lo]);
433  }
434 
435  if ( (epsilon > interp_lnL["lnL"+std::to_string(index_hi)].front()) && (epsilon < interp_lnL["lnL"+std::to_string(index_hi)].back()) )
436  {
437  z_hi = gsl_spline_eval (spline[index_hi], epsilon, acc[index_hi]);
438  }
439 
440  // Linear interpolation in Gamma.
441  double a = static_cast<double>(index_hi) - (gamma-0.4)/0.05;
442  double b = 1.0 - a;
443 
444  result = a*z_lo + b*z_hi;
445 
446  // If not in the upper part, we must be in the lower part.
447  } else {
448  const double loglikevals [4] = {-4.60517, -2.99573, -2.30259, 0.0};
449  // Gamma values belonging to the likelihood values along symmetry line (in terms of distance to 0.4).
450  double gammavals [4] = {0.0, 0.134006, 0.174898, 0.592678};
451  double distance = 0.4 - gamma;
452  // Check if point is on a vertical line with the 99% C.L. point
453  if (fabs(epsilon + 3.673776) > 1e-6)
454  {
455  // Calculate distance of point
456  double a = -3.673776 - epsilon;
457  double b = (-3.673776*gamma - 0.4*epsilon)/a;
458  double temp1 = distance;
459  distance = sqrt(a*a + temp1*temp1);
460  a = temp1/a;
461  temp1 = GSL_SIGN(-3.673776 - epsilon);
462  double temp2 = intersect_parabola_line(a, b, temp1, ppars00);
463  // CAVE: There used to be problems with distance > 1.0; these should be fixed now. Otherwise: replace by min(dist,1).
464  distance = distance/temp2;
465  gammavals[3] = 1.0;
466  gammavals[2] = intersect_parabola_line(a, b, temp1, ppars90)/temp2;
467  gammavals[1] = intersect_parabola_line(a, b, temp1, ppars95)/temp2;
468  }
469  gsl_interp_accel *acc = gsl_interp_accel_alloc ();
470  gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, 4);
471  gsl_spline_init (spline, gammavals, loglikevals, 4);
472  result = gsl_spline_eval (spline, distance, acc);
473  gsl_spline_free (spline);
474  gsl_interp_accel_free (acc);
475  }
476  }
477  // CAVE: There used to be a bug with log-likelihood > 0.0; this is fixed now, but still safeguard the result against roundoff errors.
478  return std::min(result,0.0);
479  }
480 
481 
483  // Solar model and integration routines to calculate the expected Helioscope signals //
485 
486  // SolarModel class: Provides a container to store a (tabulated) Solar model and functions to return its properties.
487  class SolarModel
488  {
489  public:
490  SolarModel();
491  SolarModel(std::string file);
492  ~SolarModel();
493  SolarModel& operator=(SolarModel&&);
494  // Delete copy constructor and assignment operator to avoid shallow copies
495  SolarModel(const SolarModel&) = delete;
496  SolarModel operator=(const SolarModel&) = delete;
497  // Min. and max. radius of the solar model file (distance r from the centre of the Sun in units of the solar radius)
498  double r_lo, r_hi;
499  // Routine to return the screening parameter kappa^2 (kappa^-1 = Debye-Hueckel radius).
500  double kappa_squared(double r);
501  // Routine to return the temperature in keV.
502  double temperature_in_keV(double r);
503  // Routine to return the plasma frequency squared.
504  double omega_pl_squared(double r);
505  private:
507  gsl_interp_accel *accel[3];
508  gsl_spline *linear_interp[3];
509  };
510 
511  SolarModel::SolarModel() {}
512  SolarModel::SolarModel(std::string file)
513  {
514  data = ASCIItableReader(file);
515  int pts = data.getnrow();
516  // Terminate GAMBIT is number of columns is wrong; i.e. the wrong solar model file format.
517  if (data.getncol() != 35)
518  {
519  DarkBit_error().raise(LOCAL_INFO, "ERROR! Solar model file '"+file+"' not compatible with GAMBIT!\n"
520  " See [arXiv:1810.07192] or example file in 'DarkBit/data/' for the correct format.");
521  }
522  data.setcolnames("mass", "radius", "temperature", "rho", "Pressure", "Luminosity", "X_H1", "X_He4", "X_He3", "X_C12", "X_C13", "X_N14", "X_N15", "X_O16", "X_O17", "X_O18", "X_Ne", "X_Na", "X_Mg", "X_Al", "X_Si", "X_P", "X_S", "X_Cl", "X_Ar",
523  "X_K", "X_Ca", "X_Sc", "X_Ti", "X_V", "X_Cr", "X_Mn", "X_Fe", "X_Co", "X_Ni");
524 
525  // Extract the radius from the files (in units of the solar radius).
526  const double* radius = &data["radius"][0];
527  r_lo = data["radius"][0];
528  r_hi = data["radius"][pts-1];
529 
530  // Extract the temperature from the files (has to be converted into keV) & calculate the screening scale kappa_s_squared.
531  // Initialise necessary variables for the screening scale calculation.
532  std::vector<double> temperature;
533  std::vector<double> kappa_s_sq;
534  std::vector<double> w_pl_sq;
535  // Multiplicative factor: (4 pi alpha_EM / atomic_mass_unit) x (1 g/cm^3) in units of keV^3
536  const double factor = 4.0*pi*alpha_EM*gsl_pow_3(1.0E+6*gev2cm)/((1.0E+9*eV2g)*atomic_mass_unit);
537  // Atomic weight of species i (exact weight if isotope is known OR estimate from average solar abundance from data if available OR estimate from natural terrestrial abundance).
538  const double A_vals [29] = {1.007825, 4.002603, 3.016029, 12.000000, 13.003355, 14.003074, 15.000109, 15.994915, 16.999132, 17.999160,
539  20.1312812, 22.989769, 24.3055, 26.9815385, 28.085, 30.973762, 32.0675, 35.4515, 36.275403, 39.0983, 40.078, 44.955908, 47.867, 50.9415, 51.9961, 54.938044, 55.845, 58.933194, 58.6934};
540  // Ionisation of species i assuming full ionisation.
541  const double Z_vals [29] = {1.0, 2.0, 2.0, 6.0, 6.0, 7.0, 7.0, 8.0, 8.0, 8.0,
542  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};
543 
544  #ifdef AXION_DEBUG_MODE
545  std::cout << "DEBUGGING INFO for solar models:\nradius/Rsol T/K kappa_s^2/keV^2 omega_pl^2/keV^2" << std::endl;
546  #endif
547 
548  // Linearly extrapolate the data in the solar model file to r = 0 if necessary.
549  if (r_lo > 0)
550  {
551  double r0 = data["radius"][0], r1 = data["radius"][1];
552  double t_intercept = (1.0E-3*K2eV)*(r0*data["temperature"][1]-r1*data["temperature"][0])/(r0-r1);
553  temperature.push_back(t_intercept);
554  double rho_intercept = (r0*data["rho"][1]-r1*data["rho"][0])/(r0-r1);
555  double sum = 0.0;
556  double ne = 0.0;
557  for (int j = 0; j < 29; j++)
558  {
559  double x_intercept = (r0* data[j+6][1]-r1* data[j+6][0])/(r0-r1);
560  double temp = x_intercept*Z_vals[j]/A_vals[j];
561  ne += temp;
562  sum += temp*(1.0 + Z_vals[j]);
563  }
564  double kss = factor*sum*rho_intercept/t_intercept;
565  kappa_s_sq.push_back(kss);
566  double wpls = factor*ne*rho_intercept/(1.0E+6*m_electron);
567  w_pl_sq.push_back(wpls);
568  #ifdef AXION_DEBUG_MODE
569  printf("%5.4f %1.6e %1.6e %1.6e\n", 0.0, t_intercept, kss, wpls);
570  #endif
571  }
572  // Calculate the necessary quantities -- T(r), kappa_s^2(r) and omega_pl^2(r) -- and store them internally.
573  for (int i = 0; i < pts; i++)
574  {
575  double sum = 0.0;
576  double ne = 0.0;
577  temperature.push_back((1.0E-3*K2eV)*data["temperature"][i]);
578  for (int j = 0; j < 29; j++)
579  {
580  double temp = data[j+6][i]*Z_vals[j]/A_vals[j];
581  ne += temp;
582  sum += temp*(1.0 + Z_vals[j]);
583  }
584  double kss = factor*sum*data["rho"][i]/temperature[i];
585  kappa_s_sq.push_back(kss);
586  double wpls = factor*ne*data["rho"][i]/(1.0E+6*m_electron);
587  w_pl_sq.push_back(wpls);
588  #ifdef AXION_DEBUG_MODE
589  printf("%5.4f %1.6e %1.6e %1.6e\n", data["radius"][i], temperature[i], kss, wpls);
590  #endif
591  }
592  // Set up the interpolating functions for temperature and screening scale.
593  accel[0] = gsl_interp_accel_alloc ();
594  linear_interp[0] = gsl_spline_alloc (gsl_interp_linear, pts);
595  const double* temp_vals = &temperature[0];
596  gsl_spline_init (linear_interp[0], radius, temp_vals, pts);
597  accel[1] = gsl_interp_accel_alloc ();
598  linear_interp[1] = gsl_spline_alloc (gsl_interp_linear, pts);
599  const double* kappa_squared_vals = &kappa_s_sq[0];
600  gsl_spline_init (linear_interp[1], radius, kappa_squared_vals, pts);
601  accel[2] = gsl_interp_accel_alloc ();
602  linear_interp[2] = gsl_spline_alloc (gsl_interp_linear, pts);
603  const double* omega_pl_squared_vals = &w_pl_sq[0];
604  gsl_spline_init (linear_interp[2], radius, omega_pl_squared_vals, pts);
605 
606  logger() << LogTags::info << "Initialisation of solar model from file '"+file+"' complete!" << std::endl;
607  logger() << LogTags::debug << "Entries in model file: " << pts << " for solar radius in [" << data["radius"][0] << ", " << data["radius"][pts-1] << "]." << EOM;
608  }
609 
610  // Move assignment operator
611  SolarModel& SolarModel::operator=(SolarModel &&model)
612  {
613  if (this != &model)
614  {
615  std::swap(data,model.data);
616  std::swap(accel,model.accel);
617  std::swap(linear_interp, model.linear_interp);
618  }
619  return *this;
620  }
621 
622  // Class destructor
623  SolarModel::~SolarModel()
624  {
625  for (auto interp : linear_interp)
626  gsl_spline_free (interp);
627  for (auto acc : accel)
628  gsl_interp_accel_free (acc);
629  }
630 
631  // Routine to return the temperature (in keV) of the zone around the distance r from the centre of the Sun.
632  double SolarModel::temperature_in_keV(double r) { return gsl_spline_eval(linear_interp[0], r, accel[0]); }
633 
634  // Routine to return the screening paramter kappa^2 in units of keV^2 (kappa^-1 = Debye-Hueckel radius).
635  double SolarModel::kappa_squared(double r)
636  {
637  // Interpolated value, directly from the Solar model.
638  return gsl_spline_eval(linear_interp[1], r, accel[1]);
639  }
640 
641  // Routine to return the plasma freqeuency squared (in keV^2) of the zone around the distance r from the centre of the Sun.
642  double SolarModel::omega_pl_squared(double r) { return gsl_spline_eval(linear_interp[2], r, accel[2]); }
643 
644  // Constant numbers for precision etc.
645  const double abs_prec = 1.0E-1, rel_prec = 1.0E-6;
646  const int method = 5;
647  // Auxillary structure for passing the model parameters to the gsl solver.
648  struct SolarModel_params1 {double erg; double rad; SolarModel* sol;};
649  struct SolarModel_params2 {double erg; double rs; SolarModel* sol;};
650  struct SolarModel_params3 {double rs; double ma0; SolarModel* sol; AxionInterpolator* eff_exp;};
651  struct SolarModel_params4 {double ma0; AxionInterpolator* eff_exp; AxionInterpolator* gaee_flux;};
652 
653  double rho_integrand (double rho, void * params)
654  {
655  // Retrieve parameters and other integration variables.
656  struct SolarModel_params1 * p1 = (struct SolarModel_params1 *)params;
657  double erg = (p1->erg);
658  double r = (p1->rad);
659  SolarModel* sol = (p1->sol);
660 
661  // Get kappa_s^2, omega_plasma^2 and the temperature.
662  double ks_sq = sol->kappa_squared(rho);
663  double w_pl_sq = sol->omega_pl_squared(rho);
664  double T_in_keV = sol->temperature_in_keV(rho);
665 
666  // Calculate the flux.
667  double x = 4.0*(erg*erg)/ks_sq;
668  double cylinder = rho*rho - r*r;
669  cylinder = rho/sqrt(cylinder);
670  double energy_factor = erg*sqrt(erg*erg - w_pl_sq)/gsl_expm1(erg/T_in_keV);
671  double rate = (ks_sq*T_in_keV)*((1.0 + 1.0/x)*gsl_log1p(x) - 1.0);
672 
673  return cylinder*energy_factor*rate;
674  }
675 
676  double rad_integrand(double rad, void * params)
677  {
678  // Retrieve and pass on parameters.
679  struct SolarModel_params2 * p2 = (struct SolarModel_params2 *)params;
680  SolarModel* sol = (p2->sol);
681  double rmax = std::min(1.0, sol->r_hi);
682  SolarModel_params1 p1 = {p2->erg, rad, sol};
683 
684  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1E6);
685  double result, error;
686 
687  gsl_function F;
688  F.function = &rho_integrand;
689  F.params = &p1;
690 
691  //gsl_set_error_handler_off();
692  gsl_integration_qag (&F, rad, rmax, 1e-2*abs_prec, 1e-2*rel_prec, 1E6, method, w, &result, &error);
693  //printf ("GSL status: %s\n", gsl_strerror (status));
694  //gsl_integration_qags(&F, rad, rmax, 1e-1*abs_prec, 1e-1*rel_prec, 1E6, w, &result, &error);
695  gsl_integration_workspace_free (w);
696 
697  result = rad*result;
698  return result;
699  }
700 
701  double erg_integrand(double erg, void * params)
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  }
731 
732  double alt_erg_integrand(double erg, void * params)
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  }
746 
747  // Provides a customised interpolation container for the CAST likelihoods.
748  class CAST_SolarModel_Interpolator
749  {
750  public:
751  CAST_SolarModel_Interpolator(std::string solar_model_gagg, std::string solar_model_gaee, std::string data_set);
752  CAST_SolarModel_Interpolator(CAST_SolarModel_Interpolator&&);
753  ~CAST_SolarModel_Interpolator();
754  // Delete copy constructor and assignment operator to avoid shallow copies
755  CAST_SolarModel_Interpolator(const CAST_SolarModel_Interpolator&) = delete;
756  CAST_SolarModel_Interpolator operator=(const CAST_SolarModel_Interpolator&) = delete;
757  std::vector<double> evaluate_gagg_contrib(double m_ax);
758  std::vector<double> evaluate_gaee_contrib(double m_ax);
759  private:
760  int n_bins;
761  ASCIItableReader gagg_data;
762  ASCIItableReader gaee_data;
763  ASCIItableReader eff_exp_data;
764  std::vector<gsl_interp_accel*> gagg_acc;
765  std::vector<gsl_interp_accel*> gaee_acc;
766  std::vector<gsl_spline*> gagg_linear_interp;
767  std::vector<gsl_spline*> gaee_linear_interp;
768  };
769 
770  // Class creators for CAST_SolarModel_Interpolator
771  // Needs path to pre-claculated data for the "default" option.
772  CAST_SolarModel_Interpolator::CAST_SolarModel_Interpolator(std::string solar_model_gagg, std::string solar_model_gaee, std::string data_set)
773  {
774  const std::string darkbitdata_path = GAMBIT_DIR "/DarkBit/data/";
775  bool user_gagg_file_missing = true, user_gaee_file_missing = true;
776  logger() << LogTags::info << "Using solar models '"+solar_model_gagg+"' (axion-photon int.) and '"+solar_model_gaee+"' (axion-electron int.) for experiment '"+data_set+"'." << EOM;
777 
778  // Check if a pre-computed a file for a given model exists.
779  user_gagg_file_missing = not(Utils::file_exists(darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gagg+"_gagg.dat"));
780  user_gaee_file_missing = not(Utils::file_exists(darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gaee+"_gaee.dat"));
781  if (not(user_gagg_file_missing)) { logger() << LogTags::info << "Found pre-calculated axion-photon counts file for experiment '"+data_set+"' and solar model '"+solar_model_gagg+"'. Skipping calculation step..." << EOM; }
782  if (not(user_gaee_file_missing)) { logger() << LogTags::info << "Found pre-calculated axion-electron counts file for experiment '"+data_set+"' and solar model '"+solar_model_gaee+"'. Skipping calculation step..." << EOM; }
783 
784  // If either file does not exists, compute it.
785  if (user_gagg_file_missing || user_gaee_file_missing)
786  {
787  // Define the list of logarithmic masses log10(m_ax/keV) for the interpolating tables.
788  const int n_mass_bins = 183;
789  const double log_masses [n_mass_bins] = {-3., -2.8, -2.6, -2.4, -2.2, -2.15, -2.1, -2.05, -2., -1.95, -1.9, -1.85, -1.8475, -1.84, -1.8325, -1.825, -1.8175, -1.81, -1.8025, -1.795, -1.7875, -1.78, -1.7725, -1.765, -1.7575, -1.75,
790  -1.7425, -1.735, -1.7275, -1.72, -1.7125, -1.705, -1.6975, -1.69, -1.6825, -1.675, -1.6675, -1.66, -1.6525, -1.645, -1.6375, -1.63, -1.6225, -1.615, -1.6075, -1.6, -1.5925, -1.585, -1.5775, -1.57,
791  -1.5625, -1.555, -1.5475, -1.54, -1.5325, -1.525, -1.5175, -1.51, -1.5025, -1.495, -1.4875, -1.48, -1.4725, -1.465, -1.4575, -1.45, -1.4425, -1.435, -1.4275, -1.42, -1.4125, -1.405, -1.3975,
792  -1.39, -1.3825, -1.375, -1.3675, -1.36, -1.3525, -1.345, -1.3375, -1.33, -1.3225, -1.315, -1.3075, -1.3, -1.2925, -1.285, -1.2775, -1.27, -1.2625, -1.255, -1.2475, -1.24, -1.2325, -1.225, -1.2175,
793  -1.21, -1.2025, -1.195, -1.1875, -1.18, -1.1725, -1.165, -1.1575, -1.15, -1.1425, -1.135, -1.1275, -1.12, -1.1125, -1.105, -1.0975, -1.09, -1.0825, -1.075, -1.0675, -1.06, -1.0525, -1.045,
794  -1.0375, -1.03, -1.0225, -1.015, -1.0075, -1., -0.9925, -0.985, -0.9775, -0.97, -0.9625, -0.955, -0.9475, -0.94, -0.9325, -0.925, -0.9175, -0.91, -0.9025, -0.895, -0.8875, -0.88, -0.8725, -0.865,
795  -0.8575, -0.85, -0.8425, -0.835, -0.8275, -0.82, -0.8125, -0.805, -0.7975, -0.79, -0.7825, -0.775, -0.7675, -0.76, -0.7525, -0.745, -0.7375, -0.73, -0.7225, -0.715, -0.7075, -0.7, -0.65, -0.6,
796  -0.55, -0.5, -0.4, -0.2, 0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2.};
797 
798  // Define quantities specific to CAST and the data set.
799  // prefactor_gagg = (keV/eV)^6 * (1 cm^2/eVcm^2) * (1 day/eVs) * (10^10 cm/eVcm) * (10^-19 eV^-1)^4 * ((9.26 m/eVm) * (9.0 T/(T/eV^2) ))^2 / (128 pi^3)
800  const double prefactor_gagg = 29302.30262;
801  // prefactor_gaee = 10^13 * (10^-19 eV^-1)^2 * ((9.26 m/eVm) * (9.0 T/(T/eV^2) ))^2 / 4
802  // 10^13 = normalisation factor of files
803  const double prefactor_gaee = 1.701818353e-4;
804  // Lowest energy bin (in keV), bin size (in keV), and max. integration radius.
805  double bin_lo = 2.0, bin_delta = 0.5, rs = 1.0;
806  // Number of bins
807  int n_bins = 10;
808  if (data_set=="CAST2007") { bin_lo = 0.8; bin_delta = 0.3; rs = 0.231738; n_bins = 20; };
809 
810  // Arrays to store the results in.
811  double gagg_counts [n_bins*n_mass_bins];
812  double gaee_counts [n_bins*n_mass_bins];
813 
814  // Load the solar model.
815  // Solar radius R_Sol and D_Sol (= 1 a.u.) in 10^10 cm.
816  double radius_sol = 6.9598, distance_sol = 1495.978707;
817  double temp = prefactor_gagg*gsl_pow_2(radius_sol/distance_sol)*radius_sol;
818 
819  SolarModel model_gagg;
820  if (user_gagg_file_missing)
821  {
822  if (Utils::file_exists(darkbitdata_path+"SolarModel_"+solar_model_gagg+".dat"))
823  {
824  model_gagg = std::move(SolarModel(darkbitdata_path+"SolarModel_"+solar_model_gagg+".dat"));
825  } else {
826  DarkBit_error().raise(LOCAL_INFO, "ERROR! No solar model file found for '"+solar_model_gagg+"'.\n"
827  " Check 'DarkBit/data' for files named 'SolarModel_*.dat' for available options *.");
828  }
829  }
830 
831  // Load and interpolate effective exposure and the data for the axion-electron spectrum (with its nasty peaks).
832  AxionInterpolator eff_exposure (darkbitdata_path+"CAST/"+data_set+"_EffectiveExposure.dat");
833  AxionInterpolator gaee_spectrum;
834  if (user_gaee_file_missing)
835  {
836  if (Utils::file_exists(darkbitdata_path+"CAST/"+"Axion_Spectrum_"+solar_model_gaee+"_gaee.dat"))
837  {
838  gaee_spectrum = std::move(AxionInterpolator(darkbitdata_path+"CAST/"+"Axion_Spectrum_"+solar_model_gaee+"_gaee.dat"));
839  } else {
840  DarkBit_error().raise(LOCAL_INFO, "ERROR! No spectrum file found for axion-electron interactions and model '"+solar_model_gaee+"'.\n"
841  " Check 'DarkBit/data' for files named 'Axion_Spectrum_*_gaee.dat' for available options *.");
842  }
843  }
844  double all_peaks [32] = {0.653029, 0.779074, 0.920547, 0.956836, 1.02042, 1.05343, 1.3497, 1.40807, 1.46949, 1.59487, 1.62314, 1.65075, 1.72461, 1.76286, 1.86037, 2.00007, 2.45281, 2.61233, 3.12669, 3.30616, 3.88237, 4.08163, 5.64394,
845  5.76064, 6.14217, 6.19863, 6.58874, 6.63942, 6.66482, 7.68441, 7.74104, 7.76785};
846 
847  // Prepare integration routine by defining the gsl functions etc.
848  gsl_function F;
849  F.function = &erg_integrand;
850  gsl_function G;
851  G.function = &alt_erg_integrand;
852 
853  double erg_lo = bin_lo, erg_hi = bin_lo;
854 
855  logger() << LogTags::info << "Calculating reference counts for dataset '"+data_set+"'..." << EOM;
856  #ifdef AXION_DEBUG_MODE
857  std::cout << "DEBUGGING INFO for solar model integration:\n"
858  "Using model '"+solar_model_gagg+"' for axion-photon interactions,"
859  "and model '"+solar_model_gaee+"' for axion-electron interactions.\n\n"
860  "coupling log10(m/eV) [erg_low/keV, erg_high/keV] log10(counts)" << std::endl;
861  #endif
862  for (int bin = 0; bin < n_bins; bin++)
863  {
864  erg_lo = erg_hi;
865  erg_hi += bin_delta;
866  gsl_integration_workspace * v = gsl_integration_workspace_alloc (1E6);
867  gsl_integration_workspace * w = gsl_integration_workspace_alloc (1E6);
868  // Only take into account the peaks relevant for the current energy bin.
869  std::vector<double> relevant_peaks;
870  relevant_peaks.push_back(erg_lo);
871  for (int i = 0; i < 32; i++)
872  {
873  double temp = all_peaks[i];
874  if ( (erg_lo < temp) && (temp < erg_hi) ) { relevant_peaks.push_back(temp); }
875  }
876  relevant_peaks.push_back(erg_hi);
877 
878  for (int i = 0; i < n_mass_bins; i++)
879  {
880  double gagg_result, gagg_error, gaee_result, gaee_error;
881  double m_ax = pow(10,log_masses[i]);
882  // Only perform integration if axion-photon counts file does not exist.
883  if (user_gagg_file_missing)
884  {
885  SolarModel_params3 p3 = {rs, m_ax, &model_gagg, &eff_exposure};
886  F.params = &p3;
887  gsl_integration_qag (&F, erg_lo, erg_hi, abs_prec, rel_prec, 1E6, method, v, &gagg_result, &gagg_error);
888 
889  #ifdef AXION_OMP_DEBUG_MODE
890  printf("gagg | % 6.4f [%3.2f, %3.2f] % 4.3e\n", log10(m_ax), erg_lo, erg_hi, log10(temp*gagg_result));
891  #endif
892 
893  gagg_counts[bin*n_mass_bins+i] = log10(temp*gagg_result);
894  }
895  // Only perform integration if axion-electron counts file does not exist.
896  if (user_gaee_file_missing)
897  {
898  SolarModel_params4 p4 = {m_ax, &eff_exposure, &gaee_spectrum};
899  G.params = &p4;
900  gsl_integration_qagp(&G, &relevant_peaks[0], relevant_peaks.size(), abs_prec, rel_prec, 1E6, w, &gaee_result, &gaee_error);
901 
902  #ifdef AXION_OMP_DEBUG_MODE
903  printf("gaee | % 6.4f [%3.2f, %3.2f] % 4.3e\n", log10(m_ax), erg_lo, erg_hi, log10(0.826*prefactor_gaee*gaee_result));
904  #endif
905 
906  // Include efficiency factor from not integrating over the full Solar disc in CAST2007 here:
907  if (data_set=="CAST2007") { gaee_result = 0.826*gaee_result; }
908  gaee_counts[bin*n_mass_bins+i] = log10(prefactor_gaee*gaee_result);
909  }
910  }
911  gsl_integration_workspace_free (v);
912  gsl_integration_workspace_free (w);
913  }
914 
915 
916  // Write the results to a file (if the file does not yet exist).
917  if (user_gagg_file_missing)
918  {
919  std::string header = "########################################################################\n"
920  "# Reference Counts for Solar Model "+solar_model_gagg+std::string(std::max(0,36-static_cast<int>(solar_model_gagg.length())),' ')+"#\n"
921  "# Column 1: log10(Axion mass in eV) #\n"
922  "# n>1: log10(Reference photon counts in energy bin n-1) #\n"
923  "########################################################################\n";
924  std::ofstream gagg_file (darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gagg+"_gagg.dat");
925  gagg_file << header;
926  gagg_file << std::fixed << std::scientific << std::setprecision(7);
927  for (int i = 0; i < n_mass_bins; i++)
928  {
929  gagg_file << log_masses[i];
930  for (int j = 0; j < n_bins; j++) { gagg_file << " " << gagg_counts[j*n_mass_bins+i]; }
931  if (i < n_mass_bins-1) { gagg_file << "\n"; }
932  }
933  gagg_file.close();
934  logger() << LogTags::info << "Output file '"+darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gagg+"_gagg.dat"+"' written for axion-photon interactions." << EOM;
935  }
936 
937  if (user_gaee_file_missing)
938  {
939  std::string header = "########################################################################\n"
940  "# Reference Counts for Solar Model "+solar_model_gaee+std::string(std::max(0,36-static_cast<int>(solar_model_gaee.length())),' ')+"#\n"
941  "# Column 1: log10(Axion mass in eV) #\n"
942  "# n>1: log10(Reference photon counts in energy bin n-1) #\n"
943  "########################################################################\n";
944  std::ofstream gaee_file (darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gaee+"_gaee.dat");
945  gaee_file << header;
946  gaee_file << std::fixed << std::scientific << std::setprecision(7);
947  for (int i = 0; i < n_mass_bins; i++)
948  {
949  gaee_file << log_masses[i];
950  for (int j = 0; j < n_bins; j++) { gaee_file << " " << gaee_counts[j*n_mass_bins+i]; }
951  if (i < n_mass_bins-1) { gaee_file << "\n"; }
952  }
953  gaee_file.close();
954  logger() << LogTags::info << "Output file '"+darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gaee+"_gagg.dat"+"' written for axion-electron interactions." << EOM;
955  }
956  }
957 
958  // Read in pre-integrated fluxes for the chosen models.
959  // 0-entry = mass values; remaining entries = counts in bins.
960  gagg_data = ASCIItableReader(darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gagg+"_gagg.dat");
961  gaee_data = ASCIItableReader(darkbitdata_path+"CAST/"+data_set+"_ReferenceCounts_"+solar_model_gaee+"_gaee.dat");
962  n_bins = gagg_data.getncol() - 1;
963 
964  // Point to the address of the first entry of masses.
965  const double* mass_gagg = &gagg_data[0][0];
966  const double* mass_gaee = &gaee_data[0][0];
967 
968  for (int bin = 0; bin < n_bins; bin++)
969  {
970  // Determine the number of interpolated mass values.
971  int gagg_pts = gagg_data[bin+1].size();
972  int gaee_pts = gaee_data[bin+1].size();
973  gagg_acc.push_back( gsl_interp_accel_alloc () );
974  gaee_acc.push_back( gsl_interp_accel_alloc () );
975  gagg_linear_interp.push_back( gsl_spline_alloc (gsl_interp_linear, gagg_pts) );
976  gaee_linear_interp.push_back( gsl_spline_alloc (gsl_interp_linear, gaee_pts) );
977  // Get flux values and initialise splines.
978  const double* flux_gagg = &gagg_data[bin+1][0];
979  const double* flux_gaee = &gaee_data[bin+1][0];
980  gsl_spline_init (gagg_linear_interp[bin], mass_gagg, flux_gagg, gagg_pts);
981  gsl_spline_init (gaee_linear_interp[bin], mass_gaee, flux_gaee, gaee_pts);
982  }
983  }
984 
985  // Move constructor
986  CAST_SolarModel_Interpolator::CAST_SolarModel_Interpolator(CAST_SolarModel_Interpolator &&interp) :
987  n_bins(std::move(interp.n_bins)),
988  gagg_data(std::move(interp.gagg_data)),
989  gaee_data(std::move(interp.gaee_data)),
990  eff_exp_data(std::move(interp.eff_exp_data)),
991  gagg_acc(std::move(interp.gagg_acc)),
992  gaee_acc(std::move(interp.gaee_acc)),
993  gagg_linear_interp(std::move(interp.gagg_linear_interp)),
994  gaee_linear_interp(std::move(interp.gaee_linear_interp))
995  {}
996 
997  // Class destructor
998  CAST_SolarModel_Interpolator::~CAST_SolarModel_Interpolator()
999  {
1000  for(auto gagg_li : gagg_linear_interp)
1001  gsl_spline_free (gagg_li);
1002  for(auto gagg_ac : gagg_acc)
1003  gsl_interp_accel_free (gagg_ac);
1004  for(auto gaee_li : gaee_linear_interp)
1005  gsl_spline_free (gaee_li);
1006  for(auto gaee_ac : gaee_acc)
1007  gsl_interp_accel_free (gaee_ac);
1008 
1009  }
1010 
1011  // Returns reference value counts for the photon-axion contribution.
1012  std::vector<double> CAST_SolarModel_Interpolator::evaluate_gagg_contrib(double m_ax)
1013  {
1014  std::vector<double> result;
1015  double lgm = log10(m_ax);
1016  // If m < 0.001 eV, just return the result for the result for the coherent limit.
1017  lgm = std::max(-3.0, lgm);
1018  // Only perform a calculation for valid masses.
1019  if (lgm < 2.0)
1020  {
1021  for (int i = 0; i < n_bins; i++) { result.push_back(gsl_spline_eval(gagg_linear_interp[i], lgm, gagg_acc[i])); }
1022  } else {
1023  for (int i = 0; i < n_bins; i++) { result.push_back(0.0); }
1024  }
1025 
1026  return result;
1027  }
1028 
1029  // Returns reference value counts for the photon-axion contribution.
1030  std::vector<double> CAST_SolarModel_Interpolator::evaluate_gaee_contrib(double m_ax)
1031  {
1032  std::vector<double> result;
1033  double lgm = log10(m_ax);
1034  // If m < 0.001 eV, just return the result for the result for the coherent limit.
1035  lgm = std::max(-3.0, lgm);
1036  // Only perform a calculation for valid masses.
1037  if (lgm < 2.0)
1038  {
1039  for (int i = 0; i < n_bins; i++) { result.push_back(gsl_spline_eval(gaee_linear_interp[i], lgm, gaee_acc[i])); }
1040  } else {
1041  for (int i = 0; i < n_bins; i++) { result.push_back(0.0); }
1042  }
1043 
1044  return result;
1045  }
1046 
1047  // Use simplified version of Gaussian likelihood from GAMBIT Utils.
1048  double gaussian_nuisance_lnL(double theo, double obs, double sigma) { return Stats::gaussian_loglikelihood(theo, obs, 0, sigma, false); }
1049 
1051  // //
1052  // Miscellaneous Theory Results //
1053  // //
1055 
1061  // Effective relatvistic degrees of freedom //
1064 
1065  // Function to provide the effective relativistic degrees of freedom (for the Standard Model).
1066  double gStar(double T)
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  }
1084 
1085  // Function to provide the effective relativistic entropic degrees of freedom (for the Standard Model).
1086  double gStar_S(double T)
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  }
1104 
1106  // QCD-axion mass relation //
1108 
1109  // Capability function to provide a simple Gaussian nuisance likelihood for
1110  // the zero-termperature mass of QCD axions.
1112  {
1114  double LambdaChi = *Param["LambdaChi"];
1115 
1116  // Results from NLO calculations (1511.02867).
1117  const double Lmu = 75.5;
1118  const double Lsigma = 0.5;
1119 
1120  result = gaussian_nuisance_lnL(Lmu, LambdaChi, Lsigma);
1121  }
1122 
1123  // Capability function to provide a simple Gaussian nuisance likelihood for
1124  // the model-independent contribution to the axion-photon coupling for QCD axions.
1126  {
1128  double CaggQCD = *Param["CaggQCD"];
1129 
1130  // Results from NLO calculations (1511.02867).
1131  const double CaggQCDmu = 1.92;
1132  const double CaggQCDsigma = 0.04;
1133 
1134  result = gaussian_nuisance_lnL(CaggQCDmu, CaggQCD, CaggQCDsigma);
1135  }
1136 
1137 
1138  // Auxillary function for QCD nuisance likelihood below.
1139  double log_chi (double T, double beta, double Tchi)
1140  {
1141  double result = 0.0;
1142  if (T > Tchi) { result = -beta*log10(T/Tchi); }
1143 
1144  return result;
1145  }
1146 
1147  // Capability function to provide a lieklihood for the temperature dependence of the QCD axion mass (doi:10.1038/nature20115).
1149  {
1151  double Tchi = *Param["Tchi"];
1152  double beta = *Param["beta"];
1153 
1154  // Results from lattice QCD (doi:10.1038/nature20115, Supplementary Material).
1155  // We normalised their findings by dividing out their value for chi(T=0) and removed its contribution to the error.
1156  const double temp_vals [20] = {100, 120, 140, 170, 200, 240, 290, 350, 420, 500, 600, 720, 860, 1000, 1200, 1500, 1800, 2100, 2500, 3000};
1157  const double log_chi_vals [20] = {0.00554625, 0.0255462, -0.0844538, -0.484454, -1.00445, -1.75445, -2.45445, -3.07445, -3.66445, -4.22445, -4.80445, -5.39445, -5.96445, -6.45445, -7.05445, -7.79445, -8.40445, -8.93445, -9.53445, -10.1545};
1158  const double log_chi_err_vals [20] = {0.014468, 0.0361846, 0.014468, 0.014468, 0.064104, 0.064104, 0.0510815, 0.0361846, 0.0510815, 0.064104, 0.0878027, 0.110042, 0.142159, 0.163124, 0.183873, 0.224965, 0.255557, 0.286023, 0.316401, 0.356804};
1159 
1160  double dummy = 0.0;
1161  for (int i = 0; i < 20; i++) { dummy = dummy + gaussian_nuisance_lnL(log_chi_vals[i], log_chi(temp_vals[i],beta,Tchi), log_chi_err_vals[i]); }
1162 
1163  result = dummy;
1164  }
1165 
1167  // //
1168  // Axion Experiments //
1169  // //
1171 
1177  // ALPS 1 experiment //
1180 
1181  // Generic functions to calculate the expected signal per frame(!) for any data run.
1182  // Input: laser power, gas coefficient nm1 = n-1; result in no. of photons.
1183  double ALPS1_signal_general(double power, double nm1, double m_ax, double gagg)
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  }
1208 
1209  // Specific capability to provide the expected signal from data run 1 (5 data frames).
1210  void calc_ALPS1_signal_vac(double &result)
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  }
1218 
1219  // Specific capability to provide the expected signal from data run 3 (8 data frames; 0.18 mbar).
1220  void calc_ALPS1_signal_gas(double &result)
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  }
1228 
1229  // General likelihood function for the ALPS 1 experiment (given expected signal s = obs - bkg).
1230  double ALPS1_lnL_general(double s, double mu, double sigma)
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  }
1235 
1236  // Capability to provide joint liklihood for all three data runs.
1237  void calc_lnL_ALPS1(double &result)
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  }
1259 
1261  // CAST experiment (vacuum runs only) //
1263 
1264  // Calculates the signal prediction for the CAST experiment (CCD detector 2004).
1265  void calc_CAST2007_signal_vac(std::vector<double> &result)
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  }
1291 
1292  // Calculates the signal prediction for the CAST experiment (all detectors in 1705.02290)
1293  void calc_CAST2017_signal_vac(std::vector<std::vector<double>> &result)
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  }
1339 
1340  // General binned Poisson likelihood for the CAST experiment.
1341  double CAST_lnL_general(std::vector<double> s, const std::vector<double> bkg_counts, const std::vector<int> sig_counts)
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  }
1354 
1355  // Capability to provide CAST likelihood for hep-ex/0702006 .
1356  void calc_lnL_CAST2007(double &result)
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  }
1378 
1379  // Capability to provide CAST likelihood for 1705.02290 .
1380  void calc_lnL_CAST2017(double &result)
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  }
1431 
1433  // Various haloscope experiments //
1435 
1436  // Capability to provide generic haloscope "signal" prediction.
1437  // All current haloscope likelihoods are approximated. We only need the predicted signal power up to a constant of proportionality.
1438  void calc_Haloscope_signal(double &result)
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  }
1452 
1456  // ADMX approximated likelihood function for data from publications from 1998 to 2009.
1457  void calc_lnL_Haloscope_ADMX1(double &result)
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  }
1494 
1495  // ADMX approximated likelihood function for data from 2018 paper (1804.05750).
1496  void calc_lnL_Haloscope_ADMX2(double &result)
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  }
1521 
1522  // University of Florida (UF) approximated likelihood function; Hagmann+ Phys. Rev. D 42, 1297(R) (1990).
1523  void calc_lnL_Haloscope_UF(double &result)
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  }
1545 
1546  // Rochester-Brookhaven-Fermi (RBF) approximated likelihood function; Phys. Rev. D 40, 3153 (1989).
1547  void calc_lnL_Haloscope_RBF(double &result)
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  }
1586 
1588  // //
1589  // Axion Cosmology //
1590  // //
1592 
1598  // Energy density in realignment axions today //
1601 
1602  /* Some auxillary functions for solving the necessary differential equations
1603  */
1604 
1605  // Provides function F1 for the change in variables time -> temperature (see 1810.07192).
1606  double SpecialFun1(double T)
1607  {
1608  // log10(T/GeV) required for interpolation.
1609  double lgT = log10(T) - 3.0;
1610  // Tabulated data: x = log10(T/GeV), y = F1(T); gR and gS from 0910.1066 .
1611  static AxionInterpolator F1 (GAMBIT_DIR "/DarkBit/data/Axion_DiffEqnFun1.dat", "linear");
1612  double res = -1.0;
1613  if ((lgT > 3.0) && (lgT < -5.0)) { res = F1.interpolate (lgT); }
1614  return res;
1615  }
1616 
1617  // Provides function F3 for the change in variables time -> temperature (see 1810.07192).
1618  double SpecialFun3(double T)
1619  {
1620  // log10(T/GeV) required for interpolation.
1621  double lgT = log10(T) - 3.0;
1622  // Tabulated data: x = log10(T/GeV), y = F3(T); gR and gS from 0910.1066 .
1623  static AxionInterpolator F3 (GAMBIT_DIR "/DarkBit/data/Axion_DiffEqnFun3.dat", "linear");
1624  double res = 0.0;
1625  if ((lgT > 3.0) && (lgT < -5.0)) { res = F3.interpolate (lgT); }
1626  return res;
1627  }
1628 
1629  // Auxillary function to calculate the Hubble parameter in a radiation-dominated universe.
1630  double hubble_rad_dom(double T)
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  }
1636 
1637  // General form of the temperature-dependent axion mass.
1638  double axion_mass_temp(double T, double beta, double Tchi)
1639  {
1640  double res = 1.0;
1641  if (T > Tchi) { res = pow(T/Tchi,-0.5*beta); }
1642  return res;
1643  }
1644 
1645  // Auxillary structure for passing the model parameters to the gsl solver.
1646  struct AxionEDT_params {double ma0; double beta; double Tchi; double thetai; double Tosc;};
1647 
1648  // Auxillary function with root Tosc, the temperature where the axion field oscillations start (defined by mA = 3H).
1649  // Note that this is only to set the temperature scale of the problem. The differential equation is solved numerically around
1650  // this point and the numerical factor in the definition is pure convention.
1651  double equation_Tosc(double T, void *params)
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  }
1664 
1665  // Capability function to solve equation_Tosc for Tosc.
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  }
1721 
1722  /* Differential equation solver to calculate the axion energy density today */
1723 
1724  // Initialise the quantities needed for the ODE solver from the gsl library.
1725  // Define the system of differential equations as a function of relative temperature.
1726  int scal_field_eq(double tau, const double y[], double f[], void *params)
1727  {
1728  struct AxionEDT_params * p = (struct AxionEDT_params *)params;
1729  double ma0 = (p->ma0);
1730  double beta = (p->beta);
1731  double Tchi= (p->Tchi);
1732  double Tosc = (p->Tosc);
1733  double thetai = (p->thetai);
1734  // f stores derivatives, y stores functions.
1735  f[0] = y[1];
1736  f[1] = -gsl_pow_2(SpecialFun1(-tau*Tosc) * ma0*axion_mass_temp(-tau*Tosc,beta,Tchi) / (hubble_rad_dom(-tau*Tosc) * (-tau))) * gsl_sf_sin(y[0]*thetai)/thetai;
1737 
1738  return GSL_SUCCESS;
1739  }
1740 
1741  // Define the Jacobian for the system of differential equations.
1742  int scal_field_eq_jac(double tau, const double y[], double *dfdy, double dfdt[], void *params)
1743  {
1744  //(void)(t); // avoid unused parameter warning.
1745  struct AxionEDT_params * p = (struct AxionEDT_params *)params;
1746  double ma0 = (p->ma0);
1747  double beta = (p->beta);
1748  double Tchi = (p->Tchi);
1749  double Tosc = (p->Tosc);
1750  double thetai = (p->thetai);
1751  gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
1752  gsl_matrix * m = &dfdy_mat.matrix;
1753  // (i, j) entries for matrix m; last entry = df[i]/df[j].
1754  gsl_matrix_set (m, 0, 0, 0);
1755  gsl_matrix_set (m, 0, 1, 1);
1756  gsl_matrix_set (m, 1, 0, -gsl_pow_2(SpecialFun1(-tau*Tosc) * ma0*axion_mass_temp(-tau*Tosc,beta,Tchi) / (hubble_rad_dom(-tau*Tosc) * (-tau))) * gsl_sf_cos(y[0]*thetai));
1757  gsl_matrix_set (m, 1, 1, -SpecialFun3(-tau*Tosc) / (-tau));
1758  dfdt[0] = 0.0;
1759  dfdt[1] = 0.0;
1760 
1761  return GSL_SUCCESS;
1762  }
1763 
1764  // Capability function to solve the differential equation for the energy density in axions today (in terms of the critical density).
1765  void RD_oh2_Axions(double &result)
1766  {
1767  using namespace Pipes::RD_oh2_Axions;
1768  double ma0 = *Param["ma0"];
1769  double beta = *Param["beta"];
1770  double Tchi = *Param["Tchi"];
1771  double thetai = *Param["thetai"];
1772  double fa = *Param["fa"];
1773  double Tosc = *Dep::AxionOscillationTemperature;
1774 
1775  if ( (thetai<-pi) || (thetai>3.0*pi) ) { DarkBit_error().raise(LOCAL_INFO, "ERROR! The parameter 'thetai' should be chosen from the interval [-pi,3pi]."); }
1776  // If thetai in (pi,3pi): map it back to its equivalent value in (-pi,pi]. This is to allow sampling around pi and easier averaging.
1777  if (thetai>pi) { thetai = thetai - 2.0*pi; }
1778 
1779  // Only do computations if thetai > 0.
1780  result = 0.0;
1781  if (fabs(thetai) > 0)
1782  {
1783  // TCMB in MeV.
1784  const double TCMB = *Dep::T_cmb*K2eV*1.0E-6;
1785  // Critical energy density today * h^2 (in eV^4).
1786  const double ede_crit_today = 3.0*2.69862E-11;
1787 
1788  struct AxionEDT_params p = {ma0, beta, Tchi, thetai, Tosc};
1789 
1790  // Function, Jacobian, number of dimensions + pointer to params.
1791  gsl_odeiv2_system sys = {scal_field_eq, scal_field_eq_jac, 2, &p};
1792  // Evolution from Temp = 1e5 x Tosc to Temp = 0.001 x Tosc.
1793  double tau2 = -0.001, tau1 = -1E5;
1794  // Initial conditions for (u and v = u') as functions of temperature:
1795  double y[2] = {1.0, 0.0};
1796  // Settings for the driver: pointer to the ODE system sys, the gsl method, initial step size,
1797  // absolute accuracy in y[0] and y[1], relative accuracy in y[0] and y[1].
1798  // Other possible choices: gsl_odeiv2_step_rk4 (classic), gsl_odeiv2_step_rk8pd, gsl_odeiv2_step_rkf45 (standard choices).
1799  gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_bsimp, -0.1*tau1, 1E-8, 1E-8);
1800 
1801  // Numerically solve ODE by continuing the integration well into the harmonic and adiabatic regime (stopping conditions
1802  // via check1 = |hat(theta)| and check2 = 3H/m need to be satisfied).
1803  double new_step;
1804  double check1 = 1.0, check2 = 1.0;
1805  int i = 0;
1806 
1807  #ifdef AXION_DEBUG_MODE
1808  std::cout << "DEBUGGING INFO for relic density calculation:\n"
1809  "#'temperature' theta dtheta/dtau" << std::endl;
1810  #endif
1811 
1812  do
1813  {
1814  i++;
1815  new_step = -pow(10.0, 1.0 + (log10(-tau2)-1.0)*i/1000.0);
1816  int status = gsl_odeiv2_driver_apply (d, &tau1, new_step, y);
1817  if (status != GSL_SUCCESS) { std::cout << "Error, return value = " << d << std::endl; }
1818  check1 = fabs(thetai)*sqrt(gsl_pow_2( fabs(y[0]) ) + gsl_pow_2( fabs((-new_step)*y[1]*hubble_rad_dom(-new_step*Tosc)/(ma0*axion_mass_temp(-new_step*Tosc,beta,Tchi))) ));
1819  check2 = 3.0*hubble_rad_dom(-new_step*Tosc)/(ma0*axion_mass_temp(-new_step*Tosc,beta,Tchi));
1820 
1821  #ifdef AXION_DEBUG_MODE
1822  std::cout << -new_step << " " << thetai*y[0] << " " << -tau2*thetai*y[1] << std::endl;
1823  #endif
1824 
1825  } while ( ((check1>1.0E-2) || (check2>1.0E-3)) && (i<1E3) );
1826 
1827  i++;
1828  if (i>=1E+3)
1829  {
1830  std::ostringstream buffer;
1831  buffer << "T_end: " << -new_step << " | theta_hat_val: " << check1 << ", theta_der: "<< -tau2*y[1]*thetai << ", 3H/m_osc: " << 3.0*hubble_rad_dom(Tosc)/(ma0*axion_mass_temp(-new_step*Tosc,beta,Tchi)) << ", 3H/m: " << check2 << " .\n";
1832  DarkBit_warning().raise(LOCAL_INFO, "WARNING! Maximum number of integration steps reached for energy density calculator!\n "+buffer.str());
1833  }
1834 
1835  // Calculate the axion energy density at the stopping point.
1836  double ede = 1E+18*gsl_pow_2(fa)*(0.5*gsl_pow_2(y[1]*thetai*hubble_rad_dom(-new_step*Tosc)*(-new_step)) + gsl_pow_2(ma0*axion_mass_temp(-new_step*Tosc,beta,Tchi))*(1.0 - gsl_sf_cos(y[0]*thetai)));
1837  // Use conservation of entropy to scale the axion energy density to its present value (relative to the critical energy density).
1838  double OmegaAh2 = ede*gsl_pow_3(TCMB/(-new_step*Tosc))*(gStar_S(TCMB)/gStar_S(-new_step*Tosc))*(1.0/axion_mass_temp(-new_step*Tosc,beta,Tchi))/ede_crit_today;
1839 
1840  gsl_odeiv2_driver_free (d);
1841 
1842  result = OmegaAh2;
1843  }
1844  }
1845 
1847  // //
1848  // Axion Astrophysics //
1849  // //
1851 
1857  // R-parameter //
1860 
1861  // Capability function to calculate the R-parameter (1512.08108).
1862  // Based and extending on Refs [11, 12, 13, 75] and 10.3204/DESY-PROC-2015-02/straniero_oscar in 1512.08108 .
1863  void calc_RParameter(double &result)
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  }
1885 
1886  // Capability function to calculate the likelihood for the R-parameter.
1887  void calc_lnL_RParameter(double &result)
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  }
1900 
1902  // White Dwarf cooling hints //
1904 
1905  // Capability function to compute the cooling likelihood of G117-B15A (1205.6180; observations from Kepler+ (2011)).
1906  void calc_lnL_WDVar_G117B15A(double &result)
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  }
1926 
1927  // Capability function to compute the cooling likelihood of R548 (1211.3389 using T = 11630 K; observations from Mukadam+ (2012)).
1928  void calc_lnL_WDVar_R548(double &result)
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  }
1947 
1948  // Capability function to compute the cooling likelihood of PG1351+489 (1605.07668 & 1406.6034; using observations from Redaelli+ (2011)).
1949  void calc_lnL_WDVar_PG1351489(double &result)
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  }
1970 
1971  // Capability function to compute the cooling likelihood of L192 (1605.06458 using l=1 & k=2; observations from Sullivan+Chote (2015)).
1972  void calc_lnL_WDVar_L192 (double &result)
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  }
1991 
1993  // SN 1987A limits (from axion-photon conversion in the B-field of the Milky Way or axion-photon decay) //
1995 
1996  // Capability function to calculate the likelihood for SN 1987A (based on data from 25 to 100 MeV photons, interpreted
1997  // by Chupp et al., Phys. Rev. Lett. 62, 505 (1989). Use 10 sec of data for conversion.
1998  void calc_lnL_SN1987A (double &result)
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  }
2010 
2012  // SN 1987A photon fluence (from axion-photon conversion in the B-field of the Milky Way) //
2014 
2015  // Capability function to calculate the photon fluence from SN 1987A as a result of axion-photon
2016  // conversion in the B-field of the Milky Way (based on arXiv:1410.3747).
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  }
2026 
2028  // Spectral distortions (H.E.S.S. telescope searches) //
2030 
2031  // Calculate the likelihood for H.E.S.S. data assuming conversion in the galaxy cluster magnetic field (GCMF, "Conservative" limits, 1311.3148).
2032  void calc_lnL_HESS_GCMF (double &result)
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  }
2049 
2050 
2051  }
2052 }
error & DarkBit_error()
START_MODEL thetai LambdaChi
Definition: Axions.hpp:49
greatScanData data
Definition: great.cpp:38
START_MODEL Tchi
Definition: Axions.hpp:35
void QCDAxion_TemperatureDependence_Nuisance_lnL(double &result)
Definition: Axions.cpp:1148
double gStar_S(double T)
Definition: Axions.cpp:1086
START_MODEL beta
Definition: Axions.hpp:35
const double m_electron
const double alpha_EM
double log_chi(double T, double beta, double Tchi)
Definition: Axions.cpp:1139
START_MODEL rs
Simple reader for ASCII tables.
void calc_lnL_CAST2007(double &result)
Definition: Axions.cpp:1356
double SpecialFun3(double T)
Definition: Axions.cpp:1618
#define LOCAL_INFO
Definition: local_info.hpp:34
int scal_field_eq(double tau, const double y[], double f[], void *params)
Definition: Axions.cpp:1726
void QCDAxion_ZeroTemperatureMass_Nuisance_lnL(double &result)
Definition: Axions.cpp:1111
void calc_Haloscope_signal(double &result)
Definition: Axions.cpp:1438
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
void RD_oh2_Axions(double &result)
Definition: Axions.cpp:1765
void calc_lnL_WDVar_L192(double &result)
Definition: Axions.cpp:1972
void calc_CAST2007_signal_vac(std::vector< double > &result)
Definition: Axions.cpp:1265
double erg_integrand(double erg, void *params)
Definition: Axions.cpp:701
void calc_lnL_Haloscope_ADMX2(double &result)
Definition: Axions.cpp:1496
void swap(Spectrum &first, Spectrum &second)
Swap resources of two Spectrum objects Note: Not a member function! This is an external function whic...
Definition: spectrum.cpp:57
void calc_RParameter(double &result)
Capabilities relating to astrophysical observations (R-parameter, H.E.S.S. telescope search...
Definition: Axions.cpp:1863
START_MODEL b
Definition: demo.hpp:235
int scal_field_eq_jac(double tau, const double y[], double *dfdy, double dfdt[], void *params)
Definition: Axions.cpp:1742
void calc_ALPS1_signal_vac(double &result)
Definition: Axions.cpp:1210
void calc_lnL_SN1987A(double &result)
Definition: Axions.cpp:1998
void QCDAxion_AxionPhotonConstant_Nuisance_lnL(double &result)
Definition: Axions.cpp:1125
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 atomic_mass_unit
void calc_lnL_Haloscope_ADMX1(double &result)
Definition: Axions.cpp:1457
General small utility functions.
void calc_ALPS1_signal_gas(double &result)
Definition: Axions.cpp:1220
void calc_CAST2017_signal_vac(std::vector< std::vector< double >> &result)
Definition: Axions.cpp:1293
const double hbar
const double m_planck_red
void calc_lnL_RParameter(double &result)
Definition: Axions.cpp:1887
GAMBIT error class.
Definition: exceptions.hpp:136
void calc_lnL_CAST2017(double &result)
Definition: Axions.cpp:1380
START_MODEL thetai thetai
Definition: Axions.hpp:49
const double sigma
Definition: SM_Z.hpp:43
Declarations of statistical utilities.
double hubble_rad_dom(double T)
Definition: Axions.cpp:1630
bool check2(const str &s1, const str &s2)
Sub-check for are_similar.
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
const double pi
EXPORT_SYMBOLS bool file_exists(const std::string &filename)
Check if a file exists.
double alt_erg_integrand(double erg, void *params)
Definition: Axions.cpp:732
const double mu
Definition: SM_Z.hpp:42
double G(const double x)
Definition: FlavBit.cpp:1589
const int method
Definition: Axions.cpp:646
double rho_integrand(double rho, void *params)
Definition: Axions.cpp:653
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
double CAST_lnL_general(std::vector< double > s, const std::vector< double > bkg_counts, const std::vector< int > sig_counts)
Definition: Axions.cpp:1341
Basic set of known mathematical and physical constants for GAMBIT.
Header file that includes all GAMBIT headers required for a module source file.
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
void calc_lnL_HESS_GCMF(double &result)
Definition: Axions.cpp:2032
double ALPS1_lnL_general(double s, double mu, double sigma)
Definition: Axions.cpp:1230
START_MODEL dNur_CMB r
void calc_lnL_WDVar_G117B15A(double &result)
Definition: Axions.cpp:1906
void calc_lnL_WDVar_R548(double &result)
Definition: Axions.cpp:1928
DS5_MSPCTM DS_INTDOF int
Rollcall header for module DarkBit.
const double gev2cm
warning & DarkBit_warning()
bool check1(const str &s1, const str &s2)
Sub-check for are_similar.
double pow(const double &a)
Outputs a^i.
double equation_Tosc(double T, void *params)
Definition: Axions.cpp:1651
double parabola(double x, const double params[])
H.E.S.S.-likelihood-related interpolation routines.
Definition: Axions.cpp:342
const double K2eV
void calc_lnL_WDVar_PG1351489(double &result)
Definition: Axions.cpp:1949
const double eV2g
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 rad_integrand(double rad, void *params)
Definition: Axions.cpp:676
void calc_PhotonFluence_SN1987A_Conversion(double &result)
Definition: Axions.cpp:2017
void calc_lnL_Haloscope_UF(double &result)
Definition: Axions.cpp:1523
double axion_mass_temp(double T, double beta, double Tchi)
Definition: Axions.cpp:1638
TODO: see if we can use this one:
Definition: Analysis.hpp:33
double E(const double x, const double y)
const double rel_prec
Definition: Axions.cpp:645
void calc_AxionOscillationTemperature(double &result)
Definition: Axions.cpp:1666
void calc_lnL_ALPS1(double &result)
Definition: Axions.cpp:1237
void calc_lnL_Haloscope_RBF(double &result)
Definition: Axions.cpp:1547
double gaussian_nuisance_lnL(double theo, double obs, double sigma)
Definition: Axions.cpp:1048
const double abs_prec
Definition: Axions.cpp:645
double SpecialFun1(double T)
Capabilities relating to axion cosmology. Currently only provides the energy density in axions today ...
Definition: Axions.cpp:1606
double intersect_parabola_line(double a, double b, double sign, const double pparams[])
Definition: Axions.cpp:345
Utility functions for DarkBit.