DarkBit.cpp

Go to the documentation of this file.

//   GAMBIT: Global and Modular BSM Inference Tool
//   *********************************************

#include "gambit/Elements/gambit_module_headers.hpp"
#include "gambit/DarkBit/DarkBit_rollcall.hpp"
#include "gambit/DarkBit/DarkBit_utils.hpp"

namespace Gambit
{
  namespace DarkBit
  {

    //
    //                 Simple DM property extractors
    //

    void mwimp_generic(double &result)
    {
      using namespace Pipes::mwimp_generic;
      result = Dep::TH_ProcessCatalog->getParticleProperty(*Dep::DarkMatter_ID).mass;
      if (result < 0.0) DarkBit_error().raise(LOCAL_INFO, "Negative WIMP mass detected.");
    }

    void sigmav_late_universe(double &result)
    {
      using namespace Pipes::sigmav_late_universe;
      std::string DMid = *Dep::DarkMatter_ID;
      std::string DMbarid = *Dep::DarkMatterConj_ID;
      TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMbarid);
      result = 0.0;
      // Add all the regular channels
      for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
          it != annProc.channelList.end(); ++it)
      {
        if ( it->nFinalStates == 2 )
        {
          // (sv)(v=0) for two-body final state
          double yield = it->genRate->bind("v")->eval(0.);
          if (yield >= 0.) result += yield;
        }
      }
      // Add invisible contributions
      double yield = annProc.genRateMisc->bind("v")->eval(0.);
      if (yield >= 0.) result += yield;
    }

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

      // set requested spectra to NULL since we don't need them
      double * SpNe=NULL,*SpNm=NULL,*SpNl=NULL;
      double * SpA=NULL,*SpE=NULL,*SpP=NULL;
      int err;

      // TODO: Add option to choose options for calculation from MicrOmegas (right now all effects are turned on)
      result = BEreq::calcSpectrum(byVal(3),byVal(SpA),byVal(SpE),byVal(SpP),byVal(SpNe),byVal(SpNm),byVal(SpNl) ,byVal(&err));

      if (err != 0 )
      {
         DarkBit_error().raise(LOCAL_INFO, "MicrOmegas spectrum calculation returned error code = " + std::to_string(err));
      }

    }

    void DM_process_from_ProcessCatalog(std::string &result)
    {
      using namespace Pipes::DM_process_from_ProcessCatalog;
      
      // Only need to check this once.
      static bool first = true;
      if (first)
      {
        // Can we find a process that is a decay?
        // (This logic used so that this capability does not depend on DarkMatterConj_ID)
        const TH_Process* p = (*Dep::TH_ProcessCatalog).find(*Dep::DarkMatter_ID);

        // If not, it's an annihilation.
        if (p == NULL) result = "annihilation";
        else result = "decay";
        first = false;
      }
      
    }


    //
    //        Extraction of local and global dark matter halo properties
    //


    double profile_gNFW(double rhos, double rs, double alpha, double beta, double gamma, double r)
    { return pow(2, (beta-gamma)/alpha)*rhos/pow(r/rs, gamma)/pow(1+pow(r/rs, alpha), (beta-gamma)/alpha); }

    double profile_Einasto(double rhos, double rs, double alpha, double r)
    { return rhos*exp((-2.0/alpha)*(pow(r/rs, alpha)-1)); }

    void GalacticHalo_gNFW(GalacticHaloProperties &result)
    {
      using namespace Pipes::GalacticHalo_gNFW;
      double rhos  = *Param["rhos"];
      double rs    = *Param["rs"];
      double r_sun = *Param["r_sun"];
      double alpha = *Param["alpha"];
      double beta  = *Param["beta"];
      double gamma = *Param["gamma"];
      result.DensityProfile = daFunk::func(profile_gNFW, rhos, rs, alpha, beta, gamma, daFunk::var("r"));
      result.r_sun = r_sun;
    }

    void GalacticHalo_Einasto(GalacticHaloProperties &result)
    {
      using namespace Pipes::GalacticHalo_Einasto;
      double rhos  = *Param["rhos"];
      double rs    = *Param["rs"];
      double r_sun = *Param["r_sun"];
      double alpha = *Param["alpha"];
      result.DensityProfile = daFunk::func(profile_Einasto, rhos, rs, alpha, daFunk::var("r"));
      result.r_sun = r_sun;
    }

    void ExtractLocalMaxwellianHalo(LocalMaxwellianHalo &result)
    {
      using namespace Pipes::ExtractLocalMaxwellianHalo;
      double rho0  = *Param["rho0"];
      double v0  = *Param["v0"];
      double vesc  = *Param["vesc"];
      double vrot  = *Param["vrot"];
      result.rho0 = rho0;
      result.v0 = v0;
      result.vesc = vesc;
      result.vrot = vrot;
    }



    //
    //                          DarkBit Unit Test
    //

    void UnitTest_DarkBit(int &result)
    {
      using namespace Pipes::UnitTest_DarkBit;
      using DarkBit_utils::gamma3bdy_limits;
      /* This function depends on all relevant DM observables (indirect and
       * direct) and dumps them into convenient files in YAML format, which
       * afterwards can be checked against the expectations.
       */

      double M_DM =
        Dep::TH_ProcessCatalog->getParticleProperty(*Dep::DarkMatter_ID).mass;
      double Gps = (*Dep::DD_couplings).gps;
      double Gpa = (*Dep::DD_couplings).gpa;
      double Gns = (*Dep::DD_couplings).gns;
      double Gna = (*Dep::DD_couplings).gna;
      double oh2 = *Dep::RD_oh2;

      std::string DMid = *Dep::DarkMatter_ID;
      std::string DMbarid = *Dep::DarkMatterConj_ID;
      TH_Process annProc = (*Dep::TH_ProcessCatalog).getProcess(DMid, DMbarid);
      daFunk::Funk spectrum = (*Dep::GA_Yield)->set("v", 0.);

      std::ostringstream filename;
      /*
      static unsigned int counter = 0;
      filename << runOptions->getValueOrDef<std::string>("UnitTest_DarkBit",
          "fileroot");
      filename << "_" << counter << ".yml";
      counter++;
      */
      filename << runOptions->getValueOrDef<std::string>("UnitTest.yaml", "filename");

      std::ofstream os;
      os.open(filename.str());
      if(os)
      {
        // Standard output.
        os << "# Direct detection couplings\n";
        os << "DDcouplings:\n";
        os << "  gps: " << Gps << "\n";
        os << "  gpa: " << Gpa << "\n";
        os << "  gns: " << Gns << "\n";
        os << "  gna: " << Gna << "\n";
        os << "\n";
        os << "# Particle masses [GeV] \n";
        os << "ParticleMasses:\n";
        os << "  Mchi: " << M_DM << "\n";
        os << "\n";
        os << "# Relic density Omega h^2\n";
        os << "RelicDensity:\n";
        os << "  oh2: " << oh2 << "\n";
        os << "\n";

        // Output gamma-ray spectrum (grid be set in YAML file).
        double x_min =
          runOptions->getValueOrDef<double>(0.1, "GA_AnnYield", "Emin");
        double x_max =
          runOptions->getValueOrDef<double>(10000, "GA_AnnYield", "Emax");
        int n = runOptions->getValueOrDef<double>(26, "GA_AnnYield", "nbins");
        // from 0.1 to 500 GeV
        std::vector<double> x = daFunk::logspace(log10(x_min), log10(x_max), n);
        std::vector<double> y = spectrum->bind("E")->vect(x);
        os << "# Annihilation spectrum dNdE [1/GeV]\n";
        os << "GammaRaySpectrum:\n";
        os << "  E: [";
        for (std::vector<double>::iterator it = x.begin(); it != x.end(); it++)
          os << *it << ", ";
        os  << "]\n";
        os << "  dNdE: [";
        for (std::vector<double>::iterator it = y.begin(); it != y.end(); it++)
          os << *it << ", ";
        os  << "]\n";
        os << std::endl;

        os << "# Annihilation rates\n";
        os << "AnnihilationRates:\n";
        for (std::vector<TH_Channel>::iterator it = annProc.channelList.begin();
            it != annProc.channelList.end(); ++it)
        {
          os << "  ";
          for (std::vector<std::string>::iterator
              jt = it->finalStateIDs.begin(); jt!=it->finalStateIDs.end(); jt++)
          {
            os << *jt << "";
          }
          if (it->finalStateIDs.size() == 2)
            os << ": " << it->genRate->bind("v")->eval(0);
          if (it->finalStateIDs.size() == 3)
          {
            double m1 = (*Dep::TH_ProcessCatalog).getParticleProperty(
                it->finalStateIDs[1]).mass;
            double m2 = (*Dep::TH_ProcessCatalog).getParticleProperty(
                it->finalStateIDs[2]).mass;
            double mass = M_DM;
            daFunk::Funk E1_low =  daFunk::func(gamma3bdy_limits<0>, daFunk::var("E"), mass, m1, m2);
            daFunk::Funk E1_high =  daFunk::func(gamma3bdy_limits<1>, daFunk::var("E"), mass, m1, m2);
            daFunk::Funk dsigmavde = it->genRate->gsl_integration("E1", E1_low, E1_high)->set_epsrel(1e-3)->set("v", 0);
            auto xgrid = daFunk::logspace(-2, 3, 1000);
            auto ygrid = daFunk::logspace(-2, 3, 1000);
            for ( size_t i = 0; i<xgrid.size(); i++ )
            {
              ygrid[i] = dsigmavde->bind("E")->eval(xgrid[i]);
            }
            auto interp = daFunk::interp("E", xgrid, ygrid);
            double svTOT = interp->gsl_integration("E", 10., 20.)->set_epsabs(1e-3)->bind()->eval();
            os << ": " << svTOT;
          }
          os << "\n";
        }
        os << std::endl;
      }
      else
      {
        logger() << "Warning: outputfile not open for writing." << EOM;
      }
      os.close();
      result = 0;
    }
  }
}