file src/Xray.cpp

[No description available] More…

Namespaces

Name
Gambit
TODO: see if we can use this one:
Gambit::DarkBit

Detailed Description

Author:

Date:

  • 2019 Sep
  • 2021 April, May
  • 2023 Aug

GAMBIT: Global and Modular BSM Inference Tool


Xray likelihoods for DarkBit. From Cirelli et al. https://arxiv.org/abs/2303.08854


Authors (add name and date if you modify):


Source code

///  GAMBIT: Global and Modular BSM Inference Tool
///  *********************************************
///  \file
///
///  Xray likelihoods for DarkBit.
///  From Cirelli et al. https://arxiv.org/abs/2303.08854
///
///  *********************************************
///
///  Authors (add name and date if you modify):
///
///  \author Patrick Stöcker
///          (stoecker@physik.rwth-aachen.de)
///  \date 2019 Sep
///
///  \author Iñigo Saez Casares
///          (inigo.saez_casares@ens-paris-saclay.fr)
///  \date 2021 April, May
///
///  \author Chris Cappiello
///          (cvc1@queensu.ca)
///  \date 2023 Aug
///
///  *********************************************

// TODO: Temporarily disabled until project is ready
/*
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_min.h>
*/

#include <gsl/gsl_odeiv2.h>

#include "gambit/Elements/gambit_module_headers.hpp"
#include "gambit/DarkBit/DarkBit_rollcall.hpp"
#include "gambit/DarkBit/DarkBit_utils.hpp"
#include "gambit/DarkBit/DarkBit_types.hpp"
#include "gambit/Utils/interp_collection.hpp"
#include "gambit/Utils/numerical_constants.hpp"
#include "gambit/Utils/util_functions.hpp"
#include "gambit/Utils/ascii_table_reader.hpp"
#include "gambit/Utils/statistics.hpp"


namespace Gambit
{
  namespace DarkBit
  {

    /// Xray loglikelihood from Cirelli et al (MISSINGREF)
    void Xray_loglikes_Cirelli(double &result)
    {
      using Interpolator2D = Utils::interp2d_gsl_collection;
      using namespace Pipes::Xray_loglikes_Cirelli;

      std::string DM_ID = Dep::WIMP_properties->name;
      std::string DMbar_ID = Dep::WIMP_properties->conjugate;
      double DM_mass = Dep::WIMP_properties->mass;
      LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
      double rho0_resc = LocalHaloParameters.rho0/0.3;
      double suppression = *Dep::ID_suppression;
      TH_Process process = Dep::TH_ProcessCatalog->getProcess(DM_ID, DMbar_ID);

      double logchisqre;
      double logchisqrmu;
      double logchisqrpi;
      double logchisqr;

      static Interpolator2D xraygride(
        "xraygride",
        GAMBIT_DIR "/DarkBit/data/Xray/xraylikelihoodse.dat",
        { "me","sigmave","loglikee"}
      );
      static Interpolator2D xraygridmu(
        "xraygridmu",
        GAMBIT_DIR "/DarkBit/data/Xray/xraylikelihoodsmu.dat",
        { "mmu","sigmavmu","loglikemu"}
      );
      static Interpolator2D xraygridpi(
        "xraygridpi",
        GAMBIT_DIR "/DarkBit/data/Xray/xraylikelihoodspi.dat",
        { "mpi","sigmavpi","loglikepi"}
      );

      double m_mine = pow(10,xraygride.x_min-3);
      double m_maxe = pow(10,xraygride.x_max-3);
      double sigmav_mine = pow(10,xraygride.y_min);
      double sigmav_maxe = pow(10,xraygride.y_max);
      double m_minmu = pow(10,xraygridmu.x_min-3);
      double m_maxmu = pow(10,xraygridmu.x_max-3);
      double sigmav_minmu = pow(10,xraygridmu.y_min);
      double sigmav_maxmu = pow(10,xraygridmu.y_max);
      double m_minpi = pow(10,xraygridpi.x_min-3);
      double m_maxpi = pow(10,xraygridpi.x_max-3);
      double sigmav_minpi = pow(10,xraygridpi.y_min);
      double sigmav_maxpi = pow(10,xraygridpi.y_max);

      double sve=0;
      double svmu=0;
      double svpi=0;
      std::vector<std::string> fs;
      std::string finalStates;
      double rate=0;

      for (std::vector<TH_Channel>::iterator it = process.channelList.begin(); it!=process.channelList.end();it++)
      {
        fs = it->finalStateIDs;
        finalStates = fs[0] + " " + fs[1];
        rate = it->genRate->bind("v")->eval(0.);
        if(fs[0]=="e+_1")
        {
          sve += rate*suppression * rho0_resc * rho0_resc;
        }
        else if(fs[0]=="e+_2")
        {
          svmu += rate*suppression * rho0_resc * rho0_resc;
        }
        else if(fs[0]=="pi+")
        {
          svpi += rate*suppression * rho0_resc * rho0_resc;
        }
      }

      if (DM_mass > m_maxe || DM_mass < m_mine)
      {
        logchisqre=-3;
      }
      else if(sve > sigmav_maxe)
      {
        logchisqre = xraygride.eval(log10(DM_mass)+3, log10(sigmav_maxe));
      }
      else if(sve < sigmav_mine)
      {
        logchisqre=-3;
      }
      else
      {
        logchisqre = xraygride.eval(log10(DM_mass)+3, log10(sve));
      }
//
      if (DM_mass > m_maxmu || DM_mass < m_minmu)
      {
        logchisqrmu=-3;
      }
      else if(svmu > sigmav_maxmu)
      {
        logchisqrmu = xraygridmu.eval(log10(DM_mass)+3, log10(sigmav_maxmu));
      }
      else if(svmu < sigmav_minmu)
      {
        logchisqrmu=-3;
      }
      else
      {
        logchisqrmu = xraygridmu.eval(log10(DM_mass)+3, log10(svmu));
      }
//
      if (DM_mass > m_maxpi || DM_mass < m_minpi)
      {
        logchisqrpi=-3;
      }
      else if(svpi > sigmav_maxpi)
      {
        logchisqrpi = xraygridpi.eval(log10(DM_mass)+3, log10(sigmav_maxpi));
      }
      else if(svpi < sigmav_minpi)
      {
        logchisqrpi=-3;
      }
      else
      {
        logchisqrpi = xraygridpi.eval(log10(DM_mass)+3, log10(svpi));
      }
      logchisqr = std::max(logchisqrpi,std::max(logchisqre,logchisqrmu));
      result = -pow(10,logchisqr)/2;
    }

///  *********************************************

// TODO: Temporarily disabled until project is ready
/*

    /////////////////////////////////////////////////////////////////
    //      Auxillary functions and classes for interpolation      //
    /////////////////////////////////////////////////////////////////

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

    // XrayInterpolator class: Provides a general 1-D interpolation container based on the gsl library.
    // Can be declared static for efficiency & easy one-time initialisation of interpolating functions.
    // (This is the twin sibling of AxionInterpolator in DarkBit/src/Axions.cpp)
    class XrayInterpolator
    {
      public:
        // Overloaded class creators for the XrayInterpolator class using the init function below.
        XrayInterpolator();
        XrayInterpolator(const std::vector<double> x, const std::vector<double> y, std::string type);
        XrayInterpolator(const std::vector<double> x, const std::vector<double> y);
        XrayInterpolator(std::string file, std::string type);
        XrayInterpolator(std::string file);
        XrayInterpolator& operator=(XrayInterpolator&&);
        // Destructor.
        ~XrayInterpolator();
        // Delete copy constructor and assignment operator to avoid shallow copies.
        XrayInterpolator(const XrayInterpolator&) = delete;
        XrayInterpolator operator=(const XrayInterpolator&) = delete;
        // Routine to access interpolated values.
        double interpolate(double x);
        // Routine to access upper and lower boundaries of available data.
        double lower();
        double upper();
      private:
        // Initialiser for the XrayInterpolator class.
        void init(std::string file, std::string type);
        void init(const std::vector<double> x, const std::vector<double> y, std::string type);
        // The gsl objects for the interpolating functions.
        gsl_interp_accel *acc;
        gsl_spline *spline;
        // Upper and lower boundaries available for the interpolating function.
        double lo;
        double up;
    };

    // Default constructor.
    XrayInterpolator::XrayInterpolator() {};

    // Initialiser for the XrayInterpolator class.
    void XrayInterpolator::init(const std::vector<double> x, const std::vector<double> y, std::string type)
    {
      int pts = x.size();
      // Get first and last value of the "x" component.
      lo = x.front();
      up = x.back();
      acc = gsl_interp_accel_alloc ();
      if (type == "cspline")
      {
        spline = gsl_spline_alloc (gsl_interp_cspline, pts);
      }
      else if (type == "linear")
      {
        spline = gsl_spline_alloc (gsl_interp_linear, pts);
      }
      else
      {
        DarkBit_error().raise(LOCAL_INFO, "ERROR! Interpolation type '"+type+"' not known to class XrayInterpolator.\n       Available types: 'linear' and 'cspline'.");
      };

      gsl_spline_init (spline, &x[0], &y[0], pts);
    };

    // Overloaded class creators for the XrayInterpolator class using the init function above.
    XrayInterpolator::XrayInterpolator(const std::vector<double> x, const std::vector<double> y, std::string type) { init(x, y, type); };
    XrayInterpolator::XrayInterpolator(const std::vector<double> x, const std::vector<double> y) { init(x, y, "linear"); };

    // Initialiser for the XrayInterpolator class.
    void XrayInterpolator::init(std::string file, std::string type)
    {
      // Check if file exists.
      if (not(Utils::file_exists(file)))
      {
        DarkBit_error().raise(LOCAL_INFO, "ERROR! File '"+file+"' not found!");
      } else {
        logger() << LogTags::debug << "Reading data from file '"+file+"' and interpolating it with '"+type+"' method." << EOM;
      };
      // Read numerical values from data file.
      ASCIItableReader tab (file);
      tab.setcolnames("x", "y");

      // for (int idx=1; idx < tab["x"].size(); idx++)
      // {
      //   std::cout << "x[" << idx << "] = " << tab["x"][idx] << "; dx = " << tab["x"][idx] -tab["x"][idx-1] << std::endl;
      //   if (tab["x"][idx] -tab["x"][idx-1] <= 0.0) std::cout << "OH NO" << std::endl;
      // }

      init(tab["x"],tab["y"],type);
    };

    // Overloaded class creators for the XrayInterpolator class using the init function above.
    XrayInterpolator::XrayInterpolator(std::string file, std::string type) { init(file, type); };
    XrayInterpolator::XrayInterpolator(std::string file) { init(file, "linear"); };

    // Move assignment operator
    XrayInterpolator& XrayInterpolator::operator=(XrayInterpolator&& interp)
    {
      if(this != &interp)
      {
        std::swap(acc,interp.acc);
        std::swap(spline,interp.spline);
        std::swap(lo,interp.lo);
        std::swap(up,interp.up);
      }
      return *this;
    }

    // Destructor
    XrayInterpolator::~XrayInterpolator()
    {
      gsl_spline_free (spline);
      gsl_interp_accel_free (acc);
    }

    // Routine to access interpolated values.
    double XrayInterpolator::interpolate(double x) { return gsl_spline_eval(spline, x, acc); };

    // Routines to return upper and lower boundaries of interpolating function
    double XrayInterpolator::lower() { return lo; };
    double XrayInterpolator::upper() { return up; };

    ////////////////////////////////////////////////////
    //               Xray likelihoods                 //
    // -- based on likehoods for sterile neutrinos -- //
    ////////////////////////////////////////////////////

    void compute_lnL_Xray_WISPy(double& result)
    {
      using namespace Pipes::compute_lnL_Xray_WISPy;

      static XrayInterpolator WISPy_bound;
      static bool first = true;
      static std::pair<double,double> xlim;

      if (first)
      {
        WISPy_bound = std::move(XrayInterpolator(GAMBIT_DIR "/DarkBit/data/WISPy_bound.dat","linear"));
        xlim.first = WISPy_bound.lower();
        xlim.second = WISPy_bound.upper();
        first = false;
      }

      double t_universe = *Dep::age_universe; // Age of the Universe in seconds

      double logm = log10(*Param["mass"]) + 9; // In "DecayingDM_mixture", the mass is given in GeV. Need to convert it into eV
      double tau = *Param["lifetime"];   // lifetime is already in untis of s. No tranformation needed.
      double frac = *Param["fraction"];
      double BR_ph = *Param["BR_ph"];

      if (logm <= xlim.first || logm >= xlim.second)
      {
        // Bound can only be applied if log10(mass) is within the range of the table.
        result = 0.0;
      }
      else
      {
        double tau_bound = pow(10.,WISPy_bound.interpolate(logm));
        bool excluded = ((1./frac)*exp(t_universe/tau)*BR_ph*tau < tau_bound);
        result = (excluded ? -9.0 : 0.0);
      }
    }

    void compute_lnL_Xray_Integral_SPI_sterile_nu(double& result)
    {
      using namespace Pipes::compute_lnL_Xray_Integral_SPI_sterile_nu;

      static XrayInterpolator sin2_2t_bound;
      static bool first = true;
      static std::pair<double,double> xlim;

      if (first)
      {
        sin2_2t_bound = std::move(XrayInterpolator(GAMBIT_DIR "/DarkBit/data/Integral_sterile_nu_bound.dat","linear"));
        xlim.first = sin2_2t_bound.lower();
        xlim.second = sin2_2t_bound.upper();
        first = false;
      }

      double t_universe = *Dep::age_universe; // Age of the Universe in seconds

      double mass = *Param["mass"] * 1e6; // In "DecayingDM_mixture", the mass is given in GeV. Need to convert it into keV
      double tau = *Param["lifetime"];   // lifetime is already in untis of s. No tranformation needed.
      double frac = *Param["fraction"];
      double BR_ph = *Param["BR_ph"];

      if (mass <= xlim.first || mass >= xlim.second)
      {
        // Bound can only be applied if the mass is within the range of the table.
        result = 0.0;
      }
      else
      {
        double sin2_2t = sin2_2t_bound.interpolate(mass);
        double tau_bound = 2./1.36038e-32 * pow(1e10*sin2_2t,-1.)*pow(mass,-5.);
        bool excluded = ((1./frac)*exp(t_universe/tau)*BR_ph*tau < tau_bound);
        result = (excluded ? -9.0 : 0.0);
      }
    }

    void compute_lnL_Xray_M31_sterile_nu(double& result)
    {
      using namespace Pipes::compute_lnL_Xray_M31_sterile_nu;

      static XrayInterpolator sin2_2t_bound;
      static bool first = true;
      static std::pair<double,double> xlim;

      if (first)
      {
        sin2_2t_bound = std::move(XrayInterpolator(GAMBIT_DIR "/DarkBit/data/M31_sterile_nu_bound.dat","linear"));
        xlim.first = sin2_2t_bound.lower();
        xlim.second = sin2_2t_bound.upper();
        first = false;
      }

      double t_universe = *Dep::age_universe; // Age of the Universe in seconds

      double mass = *Param["mass"] * 1e6; // In "DecayingDM_mixture", the mass is given in GeV. Need to convert it into keV
      double tau = *Param["lifetime"];   // lifetime is already in untis of s. No tranformation needed.
      double frac = *Param["fraction"];
      double BR_ph = *Param["BR_ph"];

      if (mass <= xlim.first || mass >= xlim.second)
      {
        // Bound can only be applied if the mass is within the range of the table.
        result = 0.0;
      }
      else
      {
        double sin2_2t = sin2_2t_bound.interpolate(mass);
        double tau_bound = 2./1.36038e-32 * pow(1e10*sin2_2t,-1.)*pow(mass,-5.);
        bool excluded = ((1./frac)*exp(t_universe/tau)*BR_ph*tau < tau_bound);
        result = (excluded ? -9.0 : 0.0);
      }
    }

    void compute_lnL_Xray_NuSTAR_sterile_nu(double& result)
    {
      using namespace Pipes::compute_lnL_Xray_NuSTAR_sterile_nu;

      static XrayInterpolator sin2_2t_bound;
      static bool first = true;
      static std::pair<double,double> xlim;

      if (first)
      {
        sin2_2t_bound = std::move(XrayInterpolator(GAMBIT_DIR "/DarkBit/data/NuSTAR_sterile_nu_bound.dat","linear"));
        xlim.first = sin2_2t_bound.lower();
        xlim.second = sin2_2t_bound.upper();
        first = false;
      }

      double t_universe = *Dep::age_universe; // Age of the Universe in seconds

      double mass = *Param["mass"] * 1e6; // In "DecayingDM_mixture", the mass is given in GeV. Need to convert it into keV
      double tau = *Param["lifetime"];   // lifetime is already in untis of s. No tranformation needed.
      double frac = *Param["fraction"];
      double BR_ph = *Param["BR_ph"];

      if (mass <= xlim.first || mass >= xlim.second)
      {
        // Bound can only be applied if the mass is within the range of the table.
        result = 0.0;
      }
      else
      {
        double sin2_2t = sin2_2t_bound.interpolate(mass);
        double tau_bound = 2./1.36038e-32 * pow(1e10*sin2_2t,-1.)*pow(mass,-5.);
        bool excluded = ((1./frac)*exp(t_universe/tau)*BR_ph*tau < tau_bound);
        result = (excluded ? -9.0 : 0.0);
      }
    }

    //------------- Numerical constants and other useful things -------------//

    // masses
    const double Mp = Gambit::m_planck; // Planck mass [GeV]

    // mathematical constants
    const double pi=Gambit::pi;

    // physical constants
    const double hbar_GeV = Gambit::hbar; // reduced Planck constant [GeV.s]
    const double cs = Gambit::s2cm; // speed of light [cm/s]
    const double Mpc_2_km = 3.0857e19; // Mpc to km

    // Minimum finite result returnable from log(double x);
    const double logmin = log(std::numeric_limits<double>::min());

    ////////////////////////////////////////////////////////////////////
    //         Support class to handle X-ray experiments              //
    ////////////////////////////////////////////////////////////////////

    //------------- Class declaration -------------//

    class Xray
    {
      public:

        Xray(std::string experiment, double J_factor);
        double solidAngle(std::vector<double> lRange, std::vector<double> bRange);
        void set_deltaOmega();
        double getDeltaOmega() const;
        double getJ() const;
        double getEmin() const;
        double getEmax() const;
        double getDeltaE() const;
        int getFluxOrigin() const;

        double flux(double const& E);
        double sigma(double const& E);
        double fluxIntegrated(double const& E);
        double sigmaIntegrated(double const& E);
        double deltaE(double const& E);
        ~Xray();

      protected:

        double m_J; // astrophysical factor for the predicted photon flux from decaying DM
        double m_Emin; // minimum energy of the observations
        double m_Emax; // maximum energy of the observations
        double m_deltaOmega; // total solid angle of observation
        double m_deltaE; // energy resolution in percentage of the energy scale
        int m_fluxOrigin; // origin of observed flux: galactic (1), extra-galactic (2) or both (3)
        std::vector<std::vector<double> > m_lRange; // observation region in galactic coordinates (degrees)
        std::vector<std::vector<double> > m_bRange;
        std::string m_experiment;
        std::map<std::string, int> m_experimentMap;
    };

    // constructor
    Xray::Xray(std::string experiment, double J_factor) : m_J(J_factor), m_experiment(experiment), m_experimentMap({{"INTEGRAL", 1}, {"HEAO", 2}})
    {
      switch(m_experimentMap[m_experiment])
      {
        case 1 :
          m_Emin = 20e3;
          m_Emax = 2e6;
          m_lRange = { {-30., 30.} };
          m_bRange = { {-15., 15.} };
          m_deltaE = 8e3;
          m_deltaOmega = 0.542068;
          m_fluxOrigin = 1;
          break;

        case 2 :
          m_Emin = 3e3;
          m_Emax = 60e3;
          m_lRange = { {58., 109.}, {238., 289.} };
          m_bRange = { {-90., -20.}, {20., 90.} };
          m_deltaE = 0.3;
          m_deltaOmega = 1.17135;
          m_fluxOrigin = 3;
          break;

        default :
          throw std::runtime_error("Wrong experiment name in Xray object");
          break;
      }
      //set_deltaOmega();
    }

    //------------- Function returning the energy dispersion of the instrument -------------//

    double Xray::deltaE (double const& E)
    {
      switch(m_experimentMap[m_experiment])
      {
        case 1 :
          return m_deltaE;
          break;

        case 2 :
          return m_deltaE*E;
          break;

        default :
          return 1.;
          break;
      }
    }

    // ----------- Functions to compute the solid angle of observation -------------//

    // auxiliary function for gsl integration
    double deltaOmega (double x[], size_t dim, void *p)
    {
      (void)(p);
      (void)(dim);
      return cos(x[1]);
    }

    // computes the solide angle for a given galactic coordinates range (in degrees)
    double Xray::solidAngle(std::vector<double> lRange, std::vector<double> bRange)
    {
      const size_t dim = 2, calls = 1e8;
      const double xl[dim] = {lRange[0]*pi/180., bRange[0]*pi/180.}, xu[dim] = {lRange[1]*pi/180., bRange[1]*pi/180.};
      double result, abserr;

      gsl_monte_plain_state *s = gsl_monte_plain_alloc(dim);
      gsl_monte_plain_init(s);
      gsl_rng *r = gsl_rng_alloc(gsl_rng_taus2);

      gsl_monte_function F;

      F.f = &deltaOmega;
      F.dim = dim;
      F.params = 0;

      gsl_monte_plain_integrate(&F, xl, xu, dim, calls, r, s, &result, &abserr);

      gsl_monte_plain_free(s);

      return result;
    }

    // sets the total solid angle of observation for the experiment
    void Xray::set_deltaOmega()
    {
      double result(0);
      for (size_t i=0; i<m_lRange.size(); ++i)
      {
        result += solidAngle(m_lRange[i], m_bRange[i]);
      }
      m_deltaOmega = result;
    }

    //------------- Functions to compute the photon flux and its standard deviation -------------//

    // differential photon flux [photons/keV/cm²/s]
    double Xray::flux(double const& E)
    {
      switch(m_experimentMap[m_experiment])
      {
        case 1 :
          // return 4.8e-8*pow(E/100e3,-1.55) + 6.6e-8*exp(-(E-50e3)/7.5e3);
          return 1.6e-7*exp(-(E-50e3)/7.7e3) + 0.92e-7*pow(E/100e3, -1.79) + 0.34e-7*pow(E/100e3, -0.95)*exp(-(E-100e3)/3411e3) + 67.3e-7/E;
          break;

        case 2 :
          return 7.877*pow(10., 0.87)*pow(E, -0.29)*exp(-(E/41.13e3))/E*m_deltaOmega;
          break;

        default :
          return 0.;
          break;
      }
    }

    // auxiliary function for gsl integration
    double flux_gsl(double x, void *p)
    {
      Xray *experiment = static_cast<Xray*>(p);
      return experiment->flux(x);
    }

    const double int_factor(1.1); // integration range = int_factor*energy dispersion instrument

    double Xray::fluxIntegrated(double const& E)
    {
      size_t n = 1e4;

      gsl_integration_workspace *w = gsl_integration_workspace_alloc(n);

      double epsabs = 0.;
      double epsrel = 1e-2;
      size_t limit = 1e3;
      double result, abserr;
      int key = 6;
      double delta = int_factor*deltaE(E);

      gsl_function F;
      F.function = &flux_gsl;
      F.params = this;

      gsl_integration_qag(&F, E-delta/2., E+delta/2., epsabs, epsrel, limit, key, w, &result, &abserr);

      gsl_integration_workspace_free(w);

      return result;
    }

    // standard deviation of the differential photon flux [photons/keV/cm²/s]
    double Xray::sigma(double const& E)
    {

      switch(m_experimentMap[m_experiment])
      {
        case 1 :
          // return sqrt(pow(pow(E/100e3,-1.55),2)*pow(0.6e-5,2) + pow(4.8e-8*1.55*pow(E/100e3,-2.55),2)*pow(0.25,2) + exp(-2*(E-50e3)/7.5e3)*pow(0.5e-8, 2.) + pow(6.6e-8, 2.)*pow((E-50e3)/pow(7.5e3, 2.), 2.)*exp(-2*(E-50e3)/7.5e3)*pow(1e3, 2.));
          return sqrt( pow(14.6e-4/E, 2) + pow(1.6e-7*(E-50e3)/7.7e3, 2) * pow(0.7e3, 2) * exp(-2.*(E-50e3)/7.7e3) + pow(0.4e-7, 2) * exp(-2.*(E-50e3)/7.7e3) + pow(0.34e-7*pow(E/100e3, -0.95), 2) * pow((E-100e3)/3411e3, 2) * exp(-(E-100e3)/3411e3) * pow(2371e3, 2));
          break;

        case 2 :
          return (flux(E)*E/m_deltaOmega)*0.02*m_deltaOmega/E;
          break;

        default :
          return 0.;
          break;
      }
    }

    // auxiliary function for gsl integration
    double sigma_gsl (double x, void *p)
    {
      Xray *experiment = static_cast<Xray*>(p);
      return pow(experiment->sigma(x), 2.);
    }


    double Xray::sigmaIntegrated(double const& E)
    {
      size_t n = 1e4;

      gsl_integration_workspace *w = gsl_integration_workspace_alloc(n);

      double epsabs = 0.;
      double epsrel = 1e-2;
      size_t limit = 1e3;
      double result, abserr;
      int key = 6;
      double delta = int_factor*deltaE(E);

      gsl_function F;
      F.function = &sigma_gsl;
      F.params = this;

      gsl_integration_qag(&F, E-delta/2., E+delta/2., epsabs, epsrel, limit, key, w, &result, &abserr);

      gsl_integration_workspace_free(w);

      return result;
    }

    //------------- Elevator functions -------------//

    double Xray::getDeltaOmega() const { return m_deltaOmega; }

    double Xray::getJ() const { return m_J; }

    double Xray::getEmin() const { return m_Emin; }

    double Xray::getEmax() const { return m_Emax; }

    double Xray::getDeltaE() const { return m_deltaE; }

    int Xray::getFluxOrigin() const { return m_fluxOrigin; }

    // destructor
    Xray::~Xray() { }

    //------------- Functions to compute the age of the Universe at a given redshift -------------//

    // useful structure
    struct cosmology_params {double OmegaM; double OmegaR; double OmegaLambda; double H0;};

    // auxiliary function for gsl integration
    double age_f (double x, void *p)
    {
      cosmology_params *params = static_cast<cosmology_params*>(p);
      double OmegaM = params->OmegaM;
      double OmegaLambda = params->OmegaLambda;
      double OmegaR = params->OmegaR;
      double OmegaK = 0.;

      return 1./sqrt( OmegaM*pow(1+x, 5.) + OmegaLambda*pow(1+x, 2.) + OmegaK*pow(1+x, 4.) + OmegaR*pow(1+x, 6.) );
    }

    // computes the age of the Universe at a given redshift ([0] age, [1] abserr)
    std::vector<double> ageUniverse (double redshift, double OmegaM, double OmegaR, double OmegaLambda, double H0)
    {
      size_t n = 1e4;

      gsl_integration_workspace *w =  gsl_integration_workspace_alloc(n);

      double epsabs = 0.;
      double epsrel = 1e-3;
      size_t limit = 1e3;
      double result, abserr;

      cosmology_params params = {OmegaM, OmegaR, OmegaLambda, H0};

      gsl_function F;
      F.function = &age_f;
      F.params = &params;

      gsl_integration_qagiu(&F, redshift, epsabs, epsrel, limit, w, &result, &abserr);

      gsl_integration_workspace_free(w);

      return {result/H0, abserr/H0};
    }

    //------------- Functions to compute X-ray likelihoods -------------//


    // useful structure
    struct XrayLikelihood_params {double mass; double tau; double gamma_ph; double fraction; Xray experiment; double OmegaDM; daFunk::Funk H_z; daFunk::Funk t_z; double ageUniverse;};

    // extra-galactic contribution to the differential photon flux [photons/eV/cm²/s]
    double dPhiEG(double const& E, XrayLikelihood_params *params)
    {
      double mass = params->mass;
      double x = mass/2./E;
      double z = x - 1.;

      if (z<0) { return 0.; }

      else
      {
      Xray experiment = params->experiment;

      double tau = params->tau, gamma_ph = params->gamma_ph, fraction = params->fraction;
      double OmegaDM = params->OmegaDM;

      daFunk::Funk H_z = params->H_z;
      boost::shared_ptr<daFunk::FunkBound> H_z_bound = H_z->bind("z");

      double H0 = H_z_bound->eval(0)/Mpc_2_km; // H0 in 1/s
      double rhoC = 3*pow(H0, 2)*pow(Mp, 2)/(8*pi)/hbar_GeV/pow(cs, 3)*1e9; // critical density un ev/cm^3
      double density = fraction*OmegaDM*rhoC;
      double H = H_z_bound->eval(x)/Mpc_2_km; // H(x) in 1/s

      daFunk::Funk t_z = params->t_z;
      boost::shared_ptr<daFunk::FunkBound> t_z_bound = t_z->bind("z");

      double t = t_z_bound->eval(z);

      return experiment.getDeltaOmega()*2.*1./(4*pi)*(gamma_ph*density*cs*exp(-t/tau))/(mass*E)/H;
      }
    }

    const double s(1./3.); // standard deviation of the gaussian for the galactic emission line = s*energy dispersion instrument

    // galactic (Milky Way) contribution to the differential photon flux [photons/eV/cm²/s]
    double dPhiG(double const& E, XrayLikelihood_params *params)
    {
      double mass = params->mass, tau = params->tau, gamma_ph = params->gamma_ph, fraction = params->fraction;
      Xray experiment = params->experiment;

      double J = experiment.getJ();

      double t0 = params->ageUniverse;

      double sigma = s*experiment.deltaE(E); // standard deviation of the gaussian modelling the enery dispersion of the instrument

      return 2.*(gamma_ph*J*fraction*exp(-t0/tau))/(4.*pi*mass)/sqrt(2*pi*sigma*sigma)*exp(-pow(E-mass/2.,2)/(2*sigma*sigma));
    }

    // total predicted differential photon flux for a given X-ray experiment [photons/eV/cm²/s]
    double XrayPrediction(double const& E, XrayLikelihood_params *params)
    {
      Xray experiment = params->experiment;
      switch(experiment.getFluxOrigin())
      {
        case 1 :
          return dPhiG(E, params);

        case 2 :
          return dPhiEG(E, params);

        case 3 :
          return dPhiG(E, params) + dPhiEG(E, params);

        default :
          throw std::runtime_error("Wrong value for m_fluxOrigin in Xray class, allowed values are 1 (galactic flux), 2 (extra-galactic flux) and 3 (both)");
          break;
      }
    }

    // auxiliary function for gsl integration
    double XrayPredictionIntegrated_gsl(double x, void *p)
    {
      XrayLikelihood_params *params = static_cast<XrayLikelihood_params*>(p);

      return XrayPrediction(x, params);
    }

    // total predicted photon flux integrated over an interval deltaE, centered around E [photons/cm²/s]
    double XrayPredictionIntegrated(double const& E, XrayLikelihood_params *params)
    {
      size_t n = 1e4;

      gsl_integration_workspace *w =  gsl_integration_workspace_alloc(n);

      double epsabs = 0.;
      double epsrel = 1e-2;
      size_t limit = 1e3;
      double result, abserr;
      int key = 6;
      Xray experiment = params->experiment;
      double delta = int_factor*experiment.deltaE(E);

      gsl_function F;
      F.function = &XrayPredictionIntegrated_gsl;
      F.params = params;

      gsl_integration_qag(&F, E-delta/2., E+delta/2., epsabs, epsrel, limit, key, w, &result, &abserr);

      gsl_integration_workspace_free(w);

      return result;
    }

    // auxiliary function for gsl minimization returning the log-likelihood for a given X-ray experiment
    double XrayLogLikelihood(double E, void *p)
    {
      XrayLikelihood_params *params = static_cast<XrayLikelihood_params*>(p);
      Xray experiment = params->experiment;

      double data = experiment.fluxIntegrated(E);
      double sigma = experiment.sigmaIntegrated(E);
      double prediction = XrayPredictionIntegrated(E, params);

      return (prediction>=data) ? -pow(data-prediction,2.)/(2.*sigma*sigma) : 0.;
    }

    // computes the energy E which minimizes the log-likelihood
    double minimizeLogLikelihood(XrayLikelihood_params *params)
    {
      int status;
      int iter = 0, max_iter = 100;
      const gsl_min_fminimizer_type *T;
      gsl_min_fminimizer *s;
      Xray experiment = params->experiment;
      double mass = params->mass;
      double Emin = experiment.getEmin(), Emax = experiment.getEmax();
      double a = Emin+experiment.deltaE(Emin), b = fmin(mass/2., Emax-experiment.deltaE(Emax));
      double m = (a+b)/2.;
      gsl_function F;

      F.function = &XrayLogLikelihood;
      F.params = params;
      T = gsl_min_fminimizer_brent;
      s = gsl_min_fminimizer_alloc(T);
      gsl_min_fminimizer_set (s, &F, m, a, b);

      do
      {
        iter++;
        status = gsl_min_fminimizer_iterate (s);

        m = gsl_min_fminimizer_x_minimum (s);
        a = gsl_min_fminimizer_x_lower (s);
        b = gsl_min_fminimizer_x_upper (s);

        status = gsl_min_test_interval (a, b, 0.1, 0);
      } while (status == GSL_CONTINUE && iter < max_iter);

      gsl_min_fminimizer_free (s);

      return m;
    }

    // Linear interpolation in lin-log space.
    double interpolate(double x, const std::vector<double> & xlist,
            const std::vector<double> & ylist, bool zerobound)
    {
        double x0, x1, y0, y1;
        int i = 1;
        if (zerobound)
        {
            if (x<xlist.front()) return 0;
            if (x>xlist.back()) return 0;
        }
        else
        {
            if (x<xlist.front()) return ylist.front();
            if (x>xlist.back()) return ylist.back();
        }
        // Find min i such that xlist[i]>=x.
        for (; xlist[i] < x; i++) {};
        x0 = xlist[i-1];
        x1 = xlist[i];
        y0 = ylist[i-1];
        y1 = ylist[i];
        // lin-vs-log interpolation for lnL vs flux
        return y0 + (y1-y0) * log(x/x0) / log(x1/x0);
    }

    void get_J_factor_INTEGRAL_CO (double &result)
    {
      using namespace Pipes::get_J_factor_INTEGRAL_CO;

      GalacticHaloProperties halo = *Dep::GalacticHalo;

      daFunk::Funk profile = halo.DensityProfile;

      std::vector<double> rho;
      auto r = daFunk::logspace(-3, 2, 100);
      double r_sun = halo.r_sun;

      for ( size_t i = 0; i<r.size(); i++ )
      {
        rho.push_back(profile->bind("r")->eval(r[i]));
      }

      std::vector<double> phi_pre;
      std::vector<double> intensity;

      BEreq::los_integral(byVal(r), byVal(rho), byVal(r_sun), phi_pre, intensity);

      auto emission = std::pair< std::vector<double>, std::vector<double> > (phi_pre, intensity);

      ASCIItableReader ROI = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/INTEGRAL/ROI_CO.txt");
      ROI.setcolnames({"phi", "weight"});
      std::vector<double> phi = ROI["phi"], weight = ROI["weight"];

      double J = 0;
      for ( size_t i = 0; i < phi.size(); i++ )
      {
        J += interpolate(phi[i], emission.first, emission.second, true)*weight[i]*3.0856775814913684e21;// J in Gev/cm^2
      }

      result = J;
    }

    // gsl error handler
    void handler (const char * reason, const char * file, int line, int gsl_errno)
    {
      if (gsl_errno == 4)
      {
        throw gsl_errno;
      }
      else { std::cerr << "gsl: " << file << ":" << line << ": ERROR: " << reason << std::endl; abort(); }
    }

    // capability function to compute the X-ray log-likelihood from the INTEGRAL experiment
    void calc_lnL_INTEGRAL_CO(double &result)
    {
      using namespace Pipes::calc_lnL_INTEGRAL_CO;

      double tau = *Param["lifetime"];

      double gamma_ph = 1/tau * *Param["BR_ph"];

      double mass = *Param["mass"]*1e9; // mass in eV

      double fraction = *Param["fraction"];

      double t0 = *Dep::age_universe;

      double J_factor = *Dep::J_factor_INTEGRAL_CO*1e9; //J in eV/cm^2

      static Xray experiment = Xray("INTEGRAL", J_factor);

      XrayLikelihood_params params = {mass, tau, gamma_ph, fraction, experiment, 0., daFunk::zero("z"), daFunk::zero("z"), t0};

      double Emin = experiment.getEmin(), Emax = experiment.getEmax(), E, lik1, lik2;

      // no constraints available above the electron threshold, we need to take into account the decay into charged particles
      if (mass >= 1e6) { result = 0; }

      else if (mass > 2.*Emin)
      {
        // modifies the gsl error handler and stores the default one
        gsl_error_handler_t *old_handler = gsl_set_error_handler (&handler);
        try
        {
          E = minimizeLogLikelihood(&params);
          result = XrayLogLikelihood(E, &params);
        }

        catch (int gsl_errno)
        {
          lik1 = XrayLogLikelihood(Emin+experiment.deltaE(Emin), &params);
          lik2 = XrayLogLikelihood(fmin(mass/2., Emax-experiment.deltaE(Emax)), &params);
          result = fmin(lik1, lik2);
        }
        // restores the default gsl error handler
        gsl_set_error_handler (old_handler);
      }

      else { result = 0; }
    }

    // function returning the decay photon flux in [photons/cm²/s] (assuming DM decays into a monochromatic line)
    // only used for the INTEGRAL_ang_b/l likelihoods
    double DecayFluxG (double gamma_ph, double fraction, double mass, double tau, double t0, double J_factor)
    {
      return 2.*(gamma_ph*fraction*exp(-t0/tau))/(4.*pi*mass)*J_factor;
    }

    void get_J_factor_INTEGRAL_ang_b (std::vector<double> &result)
    {
      using namespace Pipes::get_J_factor_INTEGRAL_ang_b;

      GalacticHaloProperties halo = *Dep::GalacticHalo;

      daFunk::Funk profile = halo.DensityProfile;

      std::vector<double> rho;
      auto r = daFunk::logspace(-3, 2, 100);
      double r_sun = halo.r_sun;

      for ( size_t i = 0; i<r.size(); i++ )
      {
        rho.push_back(profile->bind("r")->eval(r[i]));
      }

      std::vector<double> phi_pre;
      std::vector<double> intensity;

      BEreq::los_integral(byVal(r), byVal(rho), byVal(r_sun), phi_pre, intensity);

      auto emission = std::pair< std::vector<double>, std::vector<double> > (phi_pre, intensity);

      ASCIItableReader ROI_1 = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/INTEGRAL/ROI_ang_b_1.txt");
      ASCIItableReader ROI_2 = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/INTEGRAL/ROI_ang_b_2.txt");

      ROI_1.setcolnames({"phi", "weight"});
      ROI_2.setcolnames({"phi", "weight"});

      std::vector<double> phi_1 = ROI_1["phi"], weight_1 = ROI_1["weight"];
      std::vector<double> phi_2 = ROI_2["phi"], weight_2 = ROI_2["weight"];

      double J_1 = 0, J_2 = 0;

      for ( size_t i = 0; i < phi_1.size(); i++ )
      {
        J_1 += interpolate(phi_1[i], emission.first, emission.second, true)*weight_1[i]*3.0856775814913684e21; // J in Gev/cm^2
      }

      for ( size_t i = 0; i < phi_2.size(); i++ )
      {
        J_2 += interpolate(phi_2[i], emission.first, emission.second, true)*weight_2[i]*3.0856775814913684e21; // J in Gev/cm^2
      }

      result = {J_1, J_2};
    }

    void calc_lnL_INTEGRAL_ang_b (double &result)
    {
      using namespace Pipes::calc_lnL_INTEGRAL_ang_b;

      double tau = *Param["lifetime"];

      double gamma_ph = 1/tau * *Param["BR_ph"];

      double mass = *Param["mass"]; // mass in GeV

      double mass_keV = mass*1e6; // mass in keV

      double fraction = *Param["fraction"];

      double t0 = *Dep::age_universe; // Age of the Universe in seconds

      std::vector<double> J_factor = *Dep::J_factor_INTEGRAL_ang_b;

      double FluxG = DecayFluxG(gamma_ph, fraction, mass, tau, t0, J_factor[0]);

      static ASCIItableReader INTEGRAL = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/INTEGRAL/INTEGRAL_b.dat");

      INTEGRAL.setcolnames({"Emin", "Emax", "Flux", "Sigma"});

      static std::vector<double> Emin = INTEGRAL["Emin"], Emax = INTEGRAL["Emax"], Flux = INTEGRAL["Flux"], Sigma = INTEGRAL["Sigma"];

      std::vector<double> Omega = {1.6119, 4.1858};

      // no constraints available above the electron threshold, we need to take into account the decay into charged particles
      if (mass_keV >= 1e3) { result = 0; }

      else if (mass_keV < 2.**std::min_element(Emin.begin(), Emin.end())) { result = 0; }

      else
      {
        double loglik = 0.;

        double PredictedFlux, ObservedFlux, Error;

        for (size_t i = 0; i < Emin.size()-1; ++i)
        {
          PredictedFlux = ( (mass_keV >= 2*Emin[i]) && (mass_keV < 2*Emax[i]) ) ? FluxG/Omega[0] : 0;
          ObservedFlux = Flux[i];
          Error = Sigma[i];
          loglik += (PredictedFlux < ObservedFlux) ? 0 : -pow(ObservedFlux - PredictedFlux, 2)/(2.*pow(Error, 2));
        }

        result = loglik;
      }
    }

    void get_J_factor_INTEGRAL_ang_l (std::vector<double> &result)
    {
      using namespace Pipes::get_J_factor_INTEGRAL_ang_l;

      GalacticHaloProperties halo = *Dep::GalacticHalo;

      daFunk::Funk profile = halo.DensityProfile;

      std::vector<double> rho;
      auto r = daFunk::logspace(-3, 2, 100);
      double r_sun = halo.r_sun;

      for ( size_t i = 0; i<r.size(); i++ )
      {
        rho.push_back(profile->bind("r")->eval(r[i]));
      }

      std::vector<double> phi_pre;
      std::vector<double> intensity;

      BEreq::los_integral(byVal(r), byVal(rho), byVal(r_sun), phi_pre, intensity);

      auto emission = std::pair< std::vector<double>, std::vector<double> > (phi_pre, intensity);

      ASCIItableReader ROI_1 = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/INTEGRAL/ROI_ang_l_1.txt");
      ASCIItableReader ROI_2 = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/INTEGRAL/ROI_ang_l_2.txt");

      ROI_1.setcolnames({"phi", "weight"});
      ROI_2.setcolnames({"phi", "weight"});

      std::vector<double> phi_1 = ROI_1["phi"], weight_1 = ROI_1["weight"];
      std::vector<double> phi_2 = ROI_2["phi"], weight_2 = ROI_2["weight"];

      double J_1 = 0, J_2 = 0;

      for ( size_t i = 0; i < phi_1.size(); i++ )
      {
        J_1 += interpolate(phi_1[i], emission.first, emission.second, true)*weight_1[i]*3.0856775814913684e21; // J in Gev/cm^2
      }

      for ( size_t i = 0; i < phi_2.size(); i++ )
      {
        J_2 += interpolate(phi_2[i], emission.first, emission.second, true)*weight_2[i]*3.0856775814913684e21; // J in Gev/cm^2
      }

      result = {J_1, J_2};
    }

    void calc_lnL_INTEGRAL_ang_l (double &result)
    {
      using namespace Pipes::calc_lnL_INTEGRAL_ang_l;

      double tau = *Param["lifetime"];

      double gamma_ph = 1/tau * *Param["BR_ph"];

      double mass = *Param["mass"]; // mass in GeV

      double mass_keV = mass*1e6; // mass in keV

      double fraction = *Param["fraction"];

      double t0 = *Dep::age_universe; // Age of the Universe in seconds

      std::vector<double> J_factor = *Dep::J_factor_INTEGRAL_ang_l;

      double FluxG = DecayFluxG(gamma_ph, fraction, mass, tau, t0, J_factor[0]);

      static ASCIItableReader INTEGRAL = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/INTEGRAL/INTEGRAL_l.dat");

      INTEGRAL.setcolnames({"Emin", "Emax", "Flux", "Sigma"});

      static std::vector<double> Emin = INTEGRAL["Emin"], Emax = INTEGRAL["Emax"], Flux = INTEGRAL["Flux"], Sigma = INTEGRAL["Sigma"];

      std::vector<double> Omega = {1.4224, 1.7919};

      // no constraints available above the electron threshold, we need to take into account the decay into charged particles
      if (mass_keV >= 1e3) { result = 0; }

      else if (mass_keV < 2.**std::min_element(Emin.begin(), Emin.end())) { result = 0; }

      else
      {
        double loglik = 0.;

        double PredictedFlux, ObservedFlux, Error;

        for (size_t i = 0; i < Emin.size()-1; ++i)
        {
          PredictedFlux = ( (mass_keV >= 2*Emin[i]) && (mass_keV < 2*Emax[i]) ) ? FluxG/Omega[0] : 0;
          ObservedFlux = Flux[i];
          Error = Sigma[i];
          loglik += (PredictedFlux < ObservedFlux) ? 0 : -pow(ObservedFlux - PredictedFlux, 2)/(2.*pow(Error, 2));
        }

        result = loglik;
      }
    }

    // for some reason this is not giving the correct value of J, need to fix it!
    void get_J_factor_HEAO (double &result)
    {
      using namespace Pipes::get_J_factor_HEAO;

      GalacticHaloProperties halo = *Dep::GalacticHalo;

      daFunk::Funk profile = halo.DensityProfile;

      std::vector<double> rho;
      auto r = daFunk::logspace(-3, 2, 100);
      double r_sun = halo.r_sun;

      for ( size_t i = 0; i<r.size(); i++ )
      {
        rho.push_back(profile->bind("r")->eval(r[i]));
        // rho.push_back(pow(profile->bind("r")->eval(r[i]), 2));
      }

      std::vector<double> phi_pre;
      std::vector<double> intensity;

      BEreq::los_integral(byVal(r), byVal(rho), byVal(r_sun), phi_pre, intensity);

      auto emission = std::pair< std::vector<double>, std::vector<double> > (phi_pre, intensity);

      ASCIItableReader ROI = ASCIItableReader(GAMBIT_DIR "/DarkBit/data/HEAO/ROI.txt");
      ROI.setcolnames({"phi", "weight"});
      std::vector<double> phi = ROI["phi"], weight = ROI["weight"];

      double J = 0;
      for ( size_t i = 0; i < phi.size(); i++ )
      {
        J += interpolate(phi[i], emission.first, emission.second, true)*weight[i]*3.0856775814913684e21;// J in Gev/cm^2
      }

      result = J;

      // std::cout << "J = " << result/r_sun/rho_sun << std::endl;
    }

    void calc_lnL_HEAO(double &result)
    {
      using namespace Pipes::calc_lnL_HEAO;

      double OmegaDM = *Dep::Omega0_cdm;

      daFunk::Funk H_z = *Dep::H_at_z;
      daFunk::Funk t_z = *Dep::time_at_z;

      double tau = *Param["lifetime"];

      double gamma_ph = 1/tau * *Param["BR_ph"];

      double mass = *Param["mass"]*1e9; // mass in eV

      double fraction = *Param["fraction"];

      double t0 = *Dep::age_universe;

      static Xray experiment = Xray("HEAO", 9.894*1e9*3.0856775814913684e21); // J in ev / cm^2

      XrayLikelihood_params params = {mass, tau, gamma_ph, fraction, experiment, OmegaDM, H_z, t_z, t0};

      const double Emin = experiment.getEmin(), Emax = experiment.getEmax();
      double E, lik1, lik2;

      // no constraints available above the electron threshold, we need to take into account the decay into charged particles
      if (mass >= 1e6) { result = 0; }

      else if (mass > 2.*Emin)
      {
        // modifies the gsl error handler and stores the default one
        gsl_error_handler_t *old_handler = gsl_set_error_handler (&handler);
        try
        {
          E = minimizeLogLikelihood(&params);
          result = XrayLogLikelihood(E, &params);
        }

        catch (int gsl_errno)
        {
          lik1 = XrayLogLikelihood(Emin+experiment.deltaE(Emin), &params);
          lik2 = XrayLogLikelihood(fmin(mass/2., Emax-experiment.deltaE(Emax)), &params);
          result = fmin(lik1, lik2);
        }
        // restores the default gsl error handler
        gsl_set_error_handler (old_handler);
      }

      else { result = 0; }
    }
*/
  }
}

Updated on 2024-07-18 at 13:53:34 +0000