dmex.cpp

#include "Matrix.hpp"
#include "PhysicsConstants.hpp"
#include "Scattering.hpp"
#include "StandardHaloModel.hpp"
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

const std::string usage_text{R"(
Calculates DM-induced ionisation cross-section (electron scattering).

Usage:

Takes input directly from command line.
The first argument is name of the file containing Kion.
(File generated by ampsci in 'mat' format, assumed to be in atomic units).

The following input options are all optional, but must be given in order

  - Mediator mass, in MeV (may be 'inf')      [default: inf]
  - DM mass in GeV                            [default: 1.0]
  - Nuclear mass number (A)                   [default: 131]
  - cos(phi_orbit), for annual modulation     [default: 0.0]
  - alpha (Gaussian resolution parameter)     [default: 0.31]
  - beta (Gaussian resolution parameter)      [default: 0.0037]
  - gamma (detector efficiency parameter)     [default: 0.19]
  - delta (detector efficiency parameter)     [default: 3.3]
  - epsilon (detector efficiency parameter)   [default: 0.87]

e.g.,
  (1) ./dmex K_Xe_v_6_hp_orth_mat.txt > rates.txt
  (2) ./dmex K_Xe_v_6_hp_orth_mat.txt inf 0.5 > rates.txt
  (3) ./dmex K_Xe_v_6_hp_orth_mat.txt inf 1 131 0 0.31 0.0037 0.019 3.3 0.87

(1) Will calculate the rateswith all default inputs, using the 
    Kion form-factor stored in file: 'K_Xe_v_6_hp_orth_mat.txt'.
    The results are written to 'rates.txt', ready for plotting.

(2) Will calculate the rates for a heavy mediator, DM mass of 0.5 GeV.
    All other inputs will have their default values.
    The results are written to 'rates.txt', ready for plotting.

(3) Will calculate the rates, with all input options are set explicitly.
    Results written to screen.

)"};

//==============================================================================
//==============================================================================
//==============================================================================
int main(int argc, char *argv[]) {

  // Read input filename from command-line:
  const std::string in_filename = argc > 1 ? argv[1] : "";
  if (in_filename.empty()) {
    std::cout << "Error: No file given (pass as command-line arg)\n";
    std::cout << usage_text;
    return 1;
  }

  // Read in Kion from input file
  const auto [E_grid, q_grid, Kion] = Scattering::read_in_Kion(in_filename);
  if (E_grid.empty() || q_grid.empty()) {
    std::cout << in_filename << " - bad file, or file doesn't exist?\n";
    std::cout << usage_text;
    return 1;
  }

  // Check formats correct:
  assert(E_grid.size() == Kion.rows() && q_grid.size() == Kion.cols() &&
         "Input file must be in a valid format");

  // First input option: mediator mass (default: heavy mediator)
  const auto m_mediator_MeV = argc > 2 ? std::stod(argv[2]) : 1.0 / 0.0;
  std::cout << "# Mediator mass: " << m_mediator_MeV << " MeV\n";
  const auto m_mediator = m_mediator_MeV / Physics::Constants::m_e_MeV;

  // Define DM form-factor
  const auto Fx2 = [=](double q) {
    if (m_mediator == 1.0 / 0.0)
      return 1.0;
    const auto m2 = m_mediator * m_mediator;
    const auto al2 = Physics::Constants::alpha * Physics::Constants::alpha;
    const auto F = (m2 + al2) / (m2 + al2 * q * q);
    return F * F;
  };

  // Next input options: DM mass (default 1 GeV):
  const auto m_chi_GeV = argc > 3 ? std::stod(argv[3]) : 1.0;
  std::cout << "# DM mass: " << m_chi_GeV << " GeV\n";
  const auto m_chi = m_chi_GeV / Physics::Conversions::M_to_GeV;

  // Next input options: nuclear mass number, A (default 131 - for Xe):
  const double A_nuc = argc > 4 ? std::stod(argv[4]) : 131.0;

  // DM Velocity disctribution:

  // Construct SHM object: Note: this uses km/s units
  // for annual modulation: [-1,1]
  double cos_phi = argc > 5 ? std::stod(argv[5]) : 0.0;
  if (std::abs(cos_phi) > 1.0) {
    std::cout << "Warning: cos_phi=" << cos_phi
              << ", but should have |cosphi|<=1\n";
  }
  // nb: These are set to some defaults; they can be updated
  const double vobs = Astro::abs_v_obs(cos_phi, Astro::Constants::v_sun);
  const double v0 = Astro::Constants::v0;
  const double vesc = Astro::Constants::v_esc;
  std::cout << "# DM velocity distrobution:\n";
  std::cout << "# vobs = " << vobs << " km/s [cos(phi) = " << cos_phi
            << "], v0 = " << v0 << " km/s, vesc = " << vesc << " km/s\n";
  Astro::StandardHaloModel f_SHM(vobs, v0, vesc);

  // Define 'usable' f(v) function, in atomic units
  // Before: [f] = [1/(km/s)], since Integral(f)=1
  const auto Fv = [&](double v) {
    const auto v_kms = v * Physics::Conversions::V_to_kms;
    const auto fv = f_SHM.f(v_kms);
    return fv / (1.0 / Physics::Conversions::V_to_kms);
  };

  // Detector response:
  // Total observable event rate is calculated as:
  // dS(E) = eff(E) * Int[ dR(E') g(|E-E'|) dE']
  // where eff(E) is the detector efficiency
  // and g(E) is a Gaussian accounting for detector resolution
  // This is an over-simplified approach.

  // Default resolution parameters, from Xe1T paper: .....
  const double default_alpha = 0.310;
  const double default_beta = 0.0037;
  const auto alpha = argc > 6 ? std::stod(argv[6]) : default_alpha;
  const auto beta = argc > 7 ? std::stod(argv[7]) : default_beta;

  // Form detector resolution function (for XENON-like detector):
  std::cout << "# Resolution: alpha = " << alpha << ", beta = " << beta << "\n";
  auto resolution = [alpha, beta](double Eobs, double Eer) {
    const auto Eobs_keV = Eobs * Physics::Conversions::E_to_keV;
    // alpha and beta are in terms of keV - convert to atomic units
    const auto sigma = (alpha * std::sqrt(Eobs_keV) + beta * Eobs_keV) /
                       Physics::Conversions::E_to_keV;
    return Scattering::gaussian(sigma, Eobs - Eer);
  };

  // Default efficiency/acceptance parameters, from Xe1T paper: .....
  // Comes from a fit to Fig. X, from XENON collab. paper: ...
  const double default_gamma = 0.19;
  const double default_delta = 3.3;
  const double default_epsilon = 0.87;
  const auto gamma = argc > 8 ? std::stod(argv[8]) : default_gamma;
  const auto delta = argc > 9 ? std::stod(argv[9]) : default_delta;
  const auto epsilon = argc > 10 ? std::stod(argv[10]) : default_epsilon;

  // Efficiency as function of deposited energy (for XENON-like detector)
  const auto E_thresh = 0.5 / Physics::Conversions::E_to_keV;
  std::cout << "# Efficiency: gamma = " << gamma << ", delta = " << delta
            << ", epsilon = " << epsilon << "\n";
  auto efficiency = [=](double energy) {
    const auto e_keV = energy * Physics::Conversions::E_to_keV;
    // Full fit:
    // return (0.876 - 7.39e-04 * e_keV) /
    //        std::pow((1.0 + 0.104 * std::exp(-(e_keV - 1.98) / 0.360)), 2.03);
    // Simplified version (nearly as good):
    return energy < E_thresh
               ? 0.0
               : epsilon / (1.0 + gamma * std::exp(-delta * (e_keV - 2.0)));
  };

  // For simplicity, combine resolution and efficiency:
  auto detector_response = [&](double Eobs, double Eer) {
    return efficiency(Eobs) * resolution(Eobs, Eer);
  };

  //----------------------------------------------------------------------------
  // Do actual calculation:
  //----------------------------------------------------------------------------

  const int num_v_steps = 1000; // num-steps for velocity integration

  // calculate <ds.v>/dE, in units of sigma-bar_e (v and E in atomic units)
  const auto dsvde =
      Scattering::dsvdE(Kion, m_chi, Fx2, E_grid, q_grid, Fv, num_v_steps);

  // Mass of single nucleus, in kg
  const double mn_xe =
      A_nuc * (Physics::Constants::u_NMU * Physics::Constants::m_e_kg);

  // (rho / (mx*c^2)) - in units of cm^3
  const double rho_on_mxc2 =
      Astro::Constants::rhoDM_GeVcm3 / (m_chi * Physics::Conversions::M_to_GeV);

  using namespace VectorOverloads;

  // <ds.v>/dE in units: (sbe/cm^2)*cm^3/keV/day
  // Or, equivilantly, in units: cm^3/keV/day (with sigma-bar_e = 1 cm^2)
  const auto dsvde_cm3keVday = dsvde * Physics::Conversions::dsvdE_to_cm3keVday;

  // "Raw" rate, in units: (sbe/cm^2)/keV/kg/day
  const auto dRdE_keVkgday = dsvde_cm3keVday * (rho_on_mxc2 / mn_xe);

  // "Observable" rate: dSdE(E) = e(E) Int [dRdE(E') * g(|E-E'), dE']
  // in units: (sbe/cm^2)/keV/kg/day
  const auto dSdE_E =
      Scattering::convolvedRate(dRdE_keVkgday, detector_response, E_grid);

  // Output rate results to screen:
  std::cout << "\"E (keV)\" \"<ds.v>/dE ((sbe/cm^2)*cm^3/keV/day)\" \"dR/dE "
               "((sbe/cm^2)/keV/kg/day)\" \"dS/dE ((sbe/cm^2)/keV/kg/day)\"\n";
  for (std::size_t i = 0; i < E_grid.size(); ++i) {
    std::cout << E_grid.at(i) * Physics::Conversions::E_to_keV << " "
              << dsvde_cm3keVday.at(i) << " " << dRdE_keVkgday.at(i) << " "
              << dSdE_E.at(i) << "\n";
  }

  // integrate into bins: [0.5,1]
  const auto min_bin = 0.5 / Physics::Conversions::E_to_keV;
  const auto max_bin = 4.0 / Physics::Conversions::E_to_keV;
  const auto bin_width = 0.5 / Physics::Conversions::E_to_keV;
  double bin_start = min_bin;
  std::cout << "\n# Binned counts:\n";
  std::cout << "# Bin (keV)    | counts ((sbe/cm^2)/keV/kg/day)\n";
  while (bin_start < max_bin + bin_width) {
    const auto counts = Scattering::binnedCount(dSdE_E, bin_start,
                                                bin_start + bin_width, E_grid);
    printf("# [%.2f, %.2f] | %.3e\n",
           bin_start * Physics::Conversions::E_to_keV,
           (bin_start + bin_width) * Physics::Conversions::E_to_keV, counts);
    bin_start += bin_width;
  }
}