file src/SunNeutrinos.cpp

[No description available] More…

Namespaces

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

Detailed Description

Author:

Date:

  • 2015 Apr 2018 Sep
  • 2016 Aug
  • 2018 Jan, Aug

Solar neutrino likelihoods.


Authors (add name and date if you modify):


Source code

//   GAMBIT: Global and Modular BSM Inference Tool
//   *********************************************
///  \file
///
///  Solar neutrino likelihoods.
///
///  *********************************************
///
///  Authors (add name and date if you modify):
///
///  \author Pat Scott
///          (pscott@imperial.ac.uk)
///  \date 2015 Apr
///        2018 Sep
///
///  \author Sebastian Wild
///          (sebastian.wild@ph.tum.de)
///  \date 2016 Aug
///
///  \author Ankit Beniwal
///          (ankit.beniwal@adelaide.edu.au)
///  \date 2018 Jan, Aug
///
///  *********************************************

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

//#define DARKBIT_DEBUG

namespace Gambit
{

  namespace DarkBit
  {

    //////////////////////////////////////////////////////////////////////////
    //
    //            Neutrino telescope likelihoods and observables
    //
    //////////////////////////////////////////////////////////////////////////

    /*! \brief Capture rate of regular dark matter in the Sun (no v-dependent
     *         or q-dependent cross-sections) (s^-1). DarkSUSY 5 version.
     */
    void capture_rate_Sun_const_xsec_DS5(double &result)
    {
      using namespace Pipes::capture_rate_Sun_const_xsec_DS5;

      if (BEreq::cap_Sun_v0q0_isoscalar.origin()=="DarkSUSY")
        if(!(*Dep::DarkSUSY5_PointInit_LocalHalo))
          DarkBit_error().raise(LOCAL_INFO,"DarkSUSY halo model not initialized!");

      // When calculating the solar capture rate, DarkSUSY assumes that the
      // proton and neutron scattering cross-sections are the same; we
      // assume that whichever backend has been hooked up here does so too.

      result = BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp, *Dep::sigma_SI_p, *Dep::sigma_SD_p);

    }

    /*! \brief Capture rate of regular dark matter in the Sun (no v-dependent
     *         or q-dependent cross-sections) (s^-1). DarkSUSY 6 version.
     */
    void capture_rate_Sun_const_xsec(double &result)
    {
      using namespace Pipes::capture_rate_Sun_const_xsec;

      if (BEreq::cap_Sun_v0q0_isoscalar.origin()=="DarkSUSY")
        if(!(*Dep::DarkSUSY_PointInit_LocalHalo))
          DarkBit_error().raise(LOCAL_INFO,"DarkSUSY halo model not initialized!");

      LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;
      double rho0 = LocalHaloParameters.rho0;
      double rho0_eff = (*Dep::RD_fraction)*rho0;

      // When calculating the solar capture rate, DarkSUSY assumes that the
      // proton and neutron scattering cross-sections are the same; we
      // assume that whichever backend has been hooked up here does so too.

      result = BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp, rho0_eff, *Dep::sigma_SI_p, *Dep::sigma_SD_p);

    }

    ///Alternative to the darkSusy fct, using captn_specific from capgen instead
    void capture_rate_Sun_const_xsec_capgen(double &result)
    {
      using namespace Pipes::capture_rate_Sun_const_xsec_capgen;
      double resultSD;
      double resultSI;
      double maxcap;

      maxcap = BEreq::cap_sun_saturation(*Dep::mwimp);
      BEreq::cap_Sun_v0q0_isoscalar(*Dep::mwimp,*Dep::sigma_SD_p,*Dep::sigma_SI_p,resultSD,resultSI);
      result = resultSI + resultSD;

      if (maxcap < result)
      {
        result = maxcap;
      }

    }

    ///Capture rate for v^n and q^n-dependent cross sections.
    ///Isoscalar (same proton/neutron coupling)
    ///SD only couples to Hydrogen.
    ///See DirectDetection.cpp to see how to define the cross sections sigma_SD_p, sigma_SI_pi
    void capture_rate_Sun_vnqn(double &result)
    {
      using namespace Pipes::capture_rate_Sun_vnqn;

      double resultSD;
      double resultSI;
      double capped;
      int qpow;
      int vpow;
      const int nelems = 29;
      double maxcap;

      maxcap = BEreq::cap_sun_saturation(*Dep::mwimp);

      resultSI = 0e0;
      resultSD = 0e0;
      typedef map_intpair_dbl::const_iterator it_type;

      //Spin-dependent:
      for(it_type iterator = Dep::sigma_SD_p->begin();
        iterator != Dep::sigma_SD_p->end();
        iterator++)
      {
         //don't capture anything if cross section is zero or all the DM is already capped
        if((iterator->second > 1e-90) && (resultSD < maxcap))
        {
          qpow = (iterator->first).first/2 ;
          vpow =  (iterator->first).second/2;

          //Capture
          BEreq::cap_Sun_vnqn_isoscalar(*Dep::mwimp,iterator->second,1,qpow,vpow,*Dep::spinwimpx2/2.,capped);
          resultSD = resultSD+capped;
        }
      }

      //Spin independent:
      for(it_type iterator = Dep::sigma_SI_p->begin();
        iterator != Dep::sigma_SI_p->end();
        iterator++)
      {
        if((iterator->second > 1e-90) && (resultSI+resultSD < maxcap))
        {
          qpow = (iterator->first).first/2 ;
          vpow =  (iterator->first).second/2;

          //Capture
          BEreq::cap_Sun_vnqn_isoscalar(*Dep::mwimp,iterator->second,nelems,qpow,vpow,*Dep::spinwimpx2/2.,capped);
          resultSI = resultSI+capped;
        }
      }
      result = resultSI+resultSD;

      logger() << "Capgen captured: SI: " << resultSI << " SD: " << resultSD << " total: " << result << "max = " << maxcap << "\n" << EOM;

      // If capture is above saturation, return saturation value.
      if (maxcap < result)
      {
        result = maxcap;
      }

    }

    //Capture rate for Non-Relataivistic Effective Operator (NREO)
    void capture_rate_Sun_NREO(double &result)
    {
      #ifdef DARKBIT_DEBUG
        cout << "Starting capture_rate_Sun_NREO ..." << endl;
      #endif
      using namespace Pipes::capture_rate_Sun_NREO;

      double capped;
      double maxcap;
      const int niso = 16;

      maxcap = BEreq::cap_sun_saturation(Dep::WIMP_properties->mass);

      /*
      use pipe to access parameters of model (0c1...1c15) here (3.2.3 of gambit paper)

      for loop through C++ array of [0c1,0c2,...] (initialized by INI in captn_gen.cpp)
      call populate Array with the value found in the C++ array and the position in the C++ array
      */

      int coupleNum;
      int maxCouplingIndex = -1;
      if (Dep::DD_nonrel_WCs->CPTbasis == 0) // if we are using the NREffectiveTheory basis
      {
        maxCouplingIndex = 15;
      }
      else if (Dep::DD_nonrel_WCs->CPTbasis == 1) // if we are using the NREFT_CPT basis
      {
        maxCouplingIndex = 12;
      }
      else
      {
        DarkBit_error().raise(LOCAL_INFO, "Got unexpected value for Dep::DD_nonrel_WCs->CPTbasis. Should be 0 or 1.");
      }
      for(int j=0; j<maxCouplingIndex; j++)
      {
        coupleNum = j + 1; // this is the coupling number, ranges 1 to 15 (but not 2)
        if (coupleNum != 2) // 2 is not an allowed coupling constant
        {
          BEreq::captn_populate_array(Dep::DD_nonrel_WCs->c0.at(coupleNum), coupleNum, 0);
          BEreq::captn_populate_array(Dep::DD_nonrel_WCs->c1.at(coupleNum), coupleNum, 1);
        }
      }


      /*
      Code to sum over all elements in solar model simultaneously.
      The third parameter tells captn_NREO how many isotopes to sum over,
      of which captn is currently set up to sum the 16 from arxiv:1501.03729.
      */
      BEreq::captn_NREO(Dep::WIMP_properties->mass,Dep::WIMP_properties->spinx2/2.,niso,capped);

      result = capped;

      logger() << "Capgen captured: total: " << result << "max = " << maxcap << "\n" << EOM;

      // If capture is above saturation, return saturation value.
      if (maxcap < result)
      {
        result = maxcap;
      }
    }

    /*! \brief Equilibration time for capture and annihilation of dark matter
     * in the Sun (s)
     */
    void equilibration_time_Sun(double &result)
    {
      using namespace Pipes::equilibration_time_Sun;

      double sigmav = 0;
      double T_Sun_core = 1.35e-6; // Sun's core temperature (GeV)

      std::string DMid = *Dep::DarkMatter_ID;
      std::string DMbarid = *Dep::DarkMatterConj_ID;

      // Make sure that we're not trying to work with decaying DM.
      const TH_Process* p = Dep::TH_ProcessCatalog->find(DMid, DMbarid);
      if (p == NULL) DarkBit_error().raise(LOCAL_INFO, "Sorry, decaying DM is not supported yet by the DarkBit neutrino routines.");
      TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMbarid);

      // 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 = sqrt(2T/mDM)) for two-body final state
          sigmav += it->genRate->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));
        }
      }

      // Add invisible contributions
      sigmav += annProc.genRateMisc->bind("v")->eval(sqrt(2.0*T_Sun_core/(*Dep::mwimp)));

      double ca = sigmav/6.6e28 * pow(*Dep::mwimp/20.0, 1.5);
      // Scale the annihilation rate down by a factor of two if the DM is not self-conjugate
      if (not (*Dep::TH_ProcessCatalog).getProcess(*Dep::DarkMatter_ID, *Dep::DarkMatterConj_ID).isSelfConj) ca *= 0.5;
      result = pow(*Dep::capture_rate_Sun * ca, -0.5);

      // std::cout << "v = " << sqrt(2.0*T_Sun_core/(*Dep::mwimp)) << " and sigmav inside equilibration_time_Sun = " << sigmav << std::endl;
      // std::cout << "capture_rate_Sun inside equilibration_time_Sun = " << *Dep::capture_rate_Sun << std::endl;
    }

    /// Annihilation rate of dark matter in the Sun (s^-1)
    void annihilation_rate_Sun(double &result)
    {
      using namespace Pipes::annihilation_rate_Sun;
      double tt_sun = 1.5e17 / *Dep::equilibration_time_Sun;
      result = *Dep::capture_rate_Sun * 0.5 * pow(tanh(tt_sun),2.0);
    }

    /// Neutrino yield function pointer and setup
    void nuyield_from_DS(nuyield_info &result)
    {

      using namespace Pipes::nuyield_from_DS;
      double annihilation_bf[29];
      double Higgs_decay_BFs_neutral[29][3];
      double Higgs_decay_BFs_charged[15];
      double Higgs_masses_neutral[3];
      double Higgs_mass_charged;

      // Set annihilation branching fractions
      std::string DMid = *Dep::DarkMatter_ID;
      std::string DMbarid = *Dep::DarkMatterConj_ID;
      TH_Process annProc = Dep::TH_ProcessCatalog->getProcess(DMid, DMbarid);
      std::vector< std::vector<str> > neutral_channels = BEreq::get_DS_neutral_h_decay_channels();
      // the missing channel
      const std::vector<str> adhoc_chan = initVector<str>("W-", "H+");

      for (int i=0; i<29; i++)
      {
        const TH_Channel* channel = annProc.find(neutral_channels[i]);

        if (channel == NULL or i == 26) // Channel 26 has not been implemented in DarkSUSY.
        {
          annihilation_bf[i] = 0.;
        }
        else
        {
          annihilation_bf[i] = channel->genRate->bind("v")->eval(0.);
          if (i == 10) // Add W- H+ for this channel
          {
            channel = annProc.find(adhoc_chan);
            if (channel == NULL) DarkBit_error().raise(LOCAL_INFO,
                "W+H- exists in process catalog but not W-H+."
                " That's some suspiciously severe CP violation yo.");
            annihilation_bf[i] += channel->genRate->bind("v")->eval(0.);
          }
          annihilation_bf[i] /= *Dep::sigmav;

          // Check that having this channel turned on makes sense at all.
          #ifdef DARKBIT_DEBUG
            double mtot = 0;
            cout << "Particles and masses in DM annihilation final state: " << endl;
            for (auto p = neutral_channels[i].begin(); p != neutral_channels[i].end(); ++p)
            {
              double m = Dep::TH_ProcessCatalog->getParticleProperty(*p).mass;
              cout << "  " << *p << " " << m << endl;
              mtot += m;
            }
            cout << "Sqrt(s) vs total mass of final states: " << 2 * *Dep::mwimp << " vs. " << mtot << endl;
            cout << "Branching fraction in v=0 limit: " << annihilation_bf[i] << endl << endl;
            if (mtot > 2 * *Dep::mwimp and annihilation_bf[i] > 0.0)
             DarkBit_error().raise(LOCAL_INFO, "Channel is open in process catalog but should not be kinematically allowed.");
          #endif
        }

      }

      // Set Higgs masses
      if (Dep::TH_ProcessCatalog->hasParticleProperty("h0_1"))
        Higgs_masses_neutral[1] = Dep::TH_ProcessCatalog->getParticleProperty("h0_1").mass;
      else
        DarkBit_error().raise(LOCAL_INFO, "No SM-like Higgs in ProcessCatalog!");
      Higgs_masses_neutral[0] =
        Dep::TH_ProcessCatalog->hasParticleProperty("h0_2") ?
        Dep::TH_ProcessCatalog->getParticleProperty("h0_2").mass : 0.;
      Higgs_masses_neutral[2] =
        Dep::TH_ProcessCatalog->hasParticleProperty("A0") ?
        Dep::TH_ProcessCatalog->getParticleProperty("A0").mass : 0.;
      Higgs_mass_charged =
        Dep::TH_ProcessCatalog->hasParticleProperty("H+") ?
        Dep::TH_ProcessCatalog->getParticleProperty("H+").mass : 0.;

      // Find out which Higgs exist and have decay data in the process
      // catalog.
      const TH_Process* h0_decays[3];
      h0_decays[0] = Dep::TH_ProcessCatalog->find("h0_2");
      h0_decays[1] = Dep::TH_ProcessCatalog->find("h0_1");
      h0_decays[2] = Dep::TH_ProcessCatalog->find("A0");
      const TH_Process* Hplus_decays = Dep::TH_ProcessCatalog->find("H+");
      const TH_Process* Hminus_decays = Dep::TH_ProcessCatalog->find("H-");
      if (Hplus_decays != NULL and Hminus_decays == NULL) DarkBit_error().raise(
          LOCAL_INFO, "H+ decays exist in process catalog but not H-.");
      if (Hplus_decays == NULL and Hminus_decays != NULL) DarkBit_error().raise(
          LOCAL_INFO, "H- decays exist in process catalog but not H+.");

      // Set the neutral Higgs decay branching fractions
      for (int i=0; i<3; i++)       // Loop over the three neutral Higgs
      {

        // If this Higgs exists, set its decay properties.
        if (h0_decays[i] != NULL)
        {

          // Get the total decay width, for normalising partial widths to BFs.
          // TODO: Replace when BFs become directly available.
          double totalwidth = 0.0;
          for (std::vector<TH_Channel>::const_iterator
              it = h0_decays[i]->channelList.begin();
              it != h0_decays[i]->channelList.end(); ++it)
          {
            // decay width in GeV for two-body final state
            if ( it->nFinalStates == 2 ) totalwidth +=
              it->genRate->bind()->eval();
          }

          std::vector<str> neutral_channel;
          // Loop over the decay channels for neutral scalars
          for (int j=0; j<29; j++)
          {
            neutral_channel.clear();
            neutral_channel.push_back(DarkBit_utils::str_flav_to_mass((neutral_channels[j])[0]));
            neutral_channel.push_back(DarkBit_utils::str_flav_to_mass((neutral_channels[j])[1]));
            const TH_Channel* channel = h0_decays[i]->find(neutral_channel);
            // If this Higgs can decay into this channel, set the BF.
            if (channel != NULL)
            {
              Higgs_decay_BFs_neutral[j][i] = channel->genRate->bind()->eval();
              if (j == 10)          // Add W- H+ for this channel
              {
                channel = h0_decays[i]->find(adhoc_chan);
                if (channel == NULL) DarkBit_error().raise(LOCAL_INFO,
                    "W+H- exists in process catalog but not W-H+."
                    " That's some suspiciously severe CP violation yo.");
                Higgs_decay_BFs_neutral[j][i]
                  += channel->genRate->bind()->eval();
              }
              // This channel has not been implemented in DarkSUSY.
              if (j == 26) Higgs_decay_BFs_neutral[j][i] = 0.;
              Higgs_decay_BFs_neutral[j][i] /= totalwidth;
            }
            else
            {
              Higgs_decay_BFs_neutral[j][i] = 0.;
            }
          }

        }

        else
        {
          // Loop over the decay channels for neutral scalars, setting all BFs
          // for this Higgs to zero.
          for (int j=0; j<29; j++) Higgs_decay_BFs_neutral[j][i] = 0.;
        }

      }

      // If they exist, set the charged Higgs decay branching fractions
      // (DarkSUSY assumes that H+/H- decays are CP-invariant)
      if (Hplus_decays != NULL)
      {

        // Define the charged Higgs decay channels
        std::vector< std::vector<str> > charged_channels = BEreq::get_DS_charged_h_decay_channels();

        // Get the total decay width, for normalising partial widths to BFs.
        // TODO: Replace when BFs become directly available.
        double totalwidth = 0.0;
        for (std::vector<TH_Channel>::const_iterator
            it = Hplus_decays->channelList.begin();
            it != Hplus_decays->channelList.end(); ++it)
        {
          // decay width in GeV for two-body final state
          if (it->nFinalStates == 2) totalwidth += it->genRate->bind()->eval();
        }

        std::vector<str> charged_channel;
        // Loop over the decay channels for charged scalars
        for (int j=0; j<15; j++)
        {
          charged_channel.clear();
          charged_channel.push_back(DarkBit_utils::str_flav_to_mass((charged_channels[j])[0]));
          charged_channel.push_back(DarkBit_utils::str_flav_to_mass((charged_channels[j])[1]));
          const TH_Channel* channel = Hplus_decays->find(charged_channel);
          // If this Higgs can decay into this channel, set the BF.
          if (channel != NULL)
          {
            Higgs_decay_BFs_charged[j] = channel->genRate->bind()->eval();
            Higgs_decay_BFs_charged[j] /= totalwidth;
          }
          else
          {
            Higgs_decay_BFs_charged[j] = 0.;
          }
        }

      }

      else
      {
        // Loop over the decay channels for charged scalars, setting all BFs
        // for this Higgs to zero.
        for (int j=0; j<15; j++) Higgs_decay_BFs_charged[j] = 0.;
      }

      // Debug output
      #ifdef DARKBIT_DEBUG
        for (int j=0; j<29; j++)
        {
          cout<< "annihilation bfs: " << j << " " << annihilation_bf[j] << endl;
        }
        cout<< endl;
        for (int i=0; i<3; i++)
        {
          for (int j=0; j<29; j++)
          {
            cout<< "higgs neutral bfs: " << i << " " << j << " "
              << Higgs_decay_BFs_neutral[j][i] << endl;
          }
        }
        cout<< endl;
        for (int j=0; j<15; j++)
        {
          cout<< "higgs charged bfs: " << j << " "
            << Higgs_decay_BFs_charged[j] << endl;
        }
        cout<< endl;
        for (int j=0; j<3; j++)
        {
          cout<< "higgs masses neutral: " << j << " "
            << Higgs_masses_neutral[j] << endl;
        }
        cout<< endl;
        cout<< "higgs charged mass: " << Higgs_mass_charged << endl;
        cout<< endl;
        cout<< "*Dep::mwimp: " << *Dep::mwimp << endl;
        cout<< endl;
        cout<< "*Dep::sigmav: " << *Dep::sigmav << endl;
        cout<< endl;
      #endif

      // Set up DarkSUSY to do neutrino yields for this particular WIMP
      BEreq::DS_nuyield_setup(annihilation_bf, Higgs_decay_BFs_neutral,
          Higgs_decay_BFs_charged, Higgs_masses_neutral,
          Higgs_mass_charged, *Dep::mwimp);

      // Hand back the pointer to the DarkSUSY neutrino yield function
      result.pointer = BEreq::nuyield.pointer();

      //TODO: change below to >= when version numbers are available as ints
      // Treat the yield function as threadsafe only if the loaded version of DarkSUSY supports it.
      result.threadsafe = (BEreq::nuyield.version() == "5.1.3");

      // Avoid OpenMP with gfortran 6.x and later, as the implementation seems to be buggy (causes stack overflows).
      #ifdef __GNUC__
        #define GCC_VERSION (__GNUC__ * 10000 \
                     + __GNUC_MINOR__ * 100 \
                     + __GNUC_PATCHLEVEL__)
        #if GCC_VERSION > 60000
          result.threadsafe = false;
        #endif
        #undef GCC_VERSION
      #endif

    }

    /// \brief Likelihood calculators for different IceCube event samples
    /// These functions all include the likelihood of the background-only model for the respective sample.
    /// We define the final log-likelihood as delta = sum over analyses of (lnL_model - lnL_BG), conservatively
    /// forbidding delta > 0 in order to always just use the neutrino likelihood as a limit.  This ignores small
    /// low-E excesses caused by impending breakdown of approximations used in IceCube response data and the nulike
    /// likelihood at very low E. This implies conditioning on all but one parameter (e.g. the cross-section),
    /// such that including any particular IC analysis adds just *one* additional degree of freedom to the fit.
    /// @{

    /// \brief 22-string IceCube sample: predicted signal and background
    /// counts, observed counts and likelihoods.
    void IC22_full(nudata &result)
    {
      using namespace Pipes::IC22_full;
      static bool first = true;
      const double bgloglike = -808.4581;
      double sigpred, bgpred, lnLike, pval;
      int totobs;
      char experiment[300] = "IC-22";
      void* context = NULL;
      double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
      /// Option nulike_speed<int>: Speed setting for nulike backend (default 3)
      int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
      BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
          byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
          theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
      piped_invalid_point.check();
      piped_errors.check(DarkBit_error());
      piped_warnings.check(DarkBit_warning());
      result.signal = sigpred;
      result.bg = bgpred;
      result.nobs = totobs;
      result.loglike = lnLike;
      result.pvalue = pval;
      if (first)
      {
        result.bgloglike = bgloglike;
        first = false;
      }
    }

    /// \brief 79-string IceCube WH sample: predicted signal and background
    /// counts, observed counts and likelihoods.
    void IC79WH_full(nudata &result)
    {
      static bool first = true;
      using namespace Pipes::IC79WH_full;
      const double bgloglike = -11874.8689;
      double sigpred, bgpred, lnLike, pval;
      int totobs;
      char experiment[300] = "IC-79 WH";
      void* context = NULL;
      double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
      /// Option nulike_speed<int>: Speed setting for nulike backend (default 3)
      int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
      BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
          byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
          theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
      piped_invalid_point.check();
      piped_errors.check(DarkBit_error());
      piped_warnings.check(DarkBit_warning());
      result.signal = sigpred;
      result.bg = bgpred;
      result.nobs = totobs;
      result.loglike = lnLike;
      result.pvalue = pval;
      if (first)
      {
        result.bgloglike = bgloglike;
        first = false;
      }
    }

    /// \brief 79-string IceCube WL sample: predicted signal and background
    /// counts, observed counts and likelihoods.
    void IC79WL_full(nudata &result)
    {
      static bool first = true;
      using namespace Pipes::IC79WL_full;
      const double bgloglike = -1813.4503;
      double sigpred, bgpred, lnLike, pval;
      int totobs;
      char experiment[300] = "IC-79 WL";
      void* context = NULL;
      double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
      /// Option nulike_speed<int>: Speed setting for nulike backend (default 3)
      int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
      BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
          byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
          theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
      piped_invalid_point.check();
      piped_errors.check(DarkBit_error());
      piped_warnings.check(DarkBit_warning());
      result.signal = sigpred;
      result.bg = bgpred;
      result.nobs = totobs;
      result.loglike = lnLike;
      result.pvalue = pval;
      if (first)
      {
        result.bgloglike = bgloglike;
        first = false;
      }
    }

    /// \brief 79-string IceCube SL sample: predicted signal and background
    /// counts, observed counts and likelihoods.
    void IC79SL_full(nudata &result)
    {
      static bool first = true;
      using namespace Pipes::IC79SL_full;
      const double bgloglike = -5015.6474;
      double sigpred, bgpred, lnLike, pval;
      int totobs;
      char experiment[300] = "IC-79 SL";
      void* context = NULL;
      double theoryError = (*Dep::mwimp > 100.0 ? 0.05*sqrt(*Dep::mwimp*0.01) : 0.05);
      /// Option nulike_speed<int>: Speed setting for nulike backend (default 3)
      int speed = runOptions->getValueOrDef<int>(3,"nulike_speed");
      BEreq::nubounds(experiment[0], *Dep::mwimp, *Dep::annihilation_rate_Sun,
          byVal(Dep::nuyield_ptr->pointer), sigpred, bgpred, totobs, lnLike, pval, 4,
          theoryError, speed, false, 0.0, 0.0, context, Dep::nuyield_ptr->threadsafe);
      piped_invalid_point.check();
      piped_errors.check(DarkBit_error());
      piped_warnings.check(DarkBit_warning());
      result.signal = sigpred;
      result.bg = bgpred;
      result.nobs = totobs;
      result.loglike = lnLike;
      result.pvalue = pval;
      if (first)
      {
        result.bgloglike = bgloglike;
        first = false;
      }
    }
    /// @}

    /// 22-string extractor module functions
    /// @{
    void IC22_signal (double &result)   {
      result = Pipes::IC22_signal ::Dep::IC22_data->signal;      }
    void IC22_bg     (double &result)   {
      result = Pipes::IC22_bg     ::Dep::IC22_data->bg;          }
    void IC22_nobs   (int    &result)   {
      result = Pipes::IC22_nobs   ::Dep::IC22_data->nobs;        }
    void IC22_loglike(double &result)   {
      result = Pipes::IC22_loglike::Dep::IC22_data->loglike;     }
    void IC22_bgloglike(double &result) {
      result = Pipes::IC22_bgloglike::Dep::IC22_data->bgloglike; }
    void IC22_pvalue (double &result)   {
      result = Pipes::IC22_pvalue ::Dep::IC22_data->pvalue;      }
    /// @}

    /// 79-string WH extractor module functions
    /// @{
    void IC79WH_signal (double &result)   {
      result = Pipes::IC79WH_signal ::Dep::IC79WH_data->signal;      }
    void IC79WH_bg     (double &result)   {
      result = Pipes::IC79WH_bg     ::Dep::IC79WH_data->bg;          }
    void IC79WH_nobs   (int    &result)   {
      result = Pipes::IC79WH_nobs   ::Dep::IC79WH_data->nobs;        }
    void IC79WH_loglike(double &result)   {
      result = Pipes::IC79WH_loglike::Dep::IC79WH_data->loglike;     }
    void IC79WH_bgloglike(double &result) {
      result = Pipes::IC79WH_bgloglike::Dep::IC79WH_data->bgloglike; }
    void IC79WH_pvalue (double &result)   {
      result = Pipes::IC79WH_pvalue ::Dep::IC79WH_data->pvalue;      }
    /// @}

    /// 79-string WL extractor module functions
    /// @{
    void IC79WL_signal (double &result)   {
      result = Pipes::IC79WL_signal ::Dep::IC79WL_data->signal;      }
    void IC79WL_bg     (double &result)   {
      result = Pipes::IC79WL_bg     ::Dep::IC79WL_data->bg;          }
    void IC79WL_nobs   (int    &result)   {
      result = Pipes::IC79WL_nobs   ::Dep::IC79WL_data->nobs;        }
    void IC79WL_loglike(double &result)   {
      result = Pipes::IC79WL_loglike::Dep::IC79WL_data->loglike;     }
    void IC79WL_bgloglike(double &result) {
      result = Pipes::IC79WL_bgloglike::Dep::IC79WL_data->bgloglike; }
    void IC79WL_pvalue (double &result)   {
      result = Pipes::IC79WL_pvalue ::Dep::IC79WL_data->pvalue;      }
    /// @}

    /// 79-string SL extractor module functions
    /// @{
    void IC79SL_signal (double &result)   {
      result = Pipes::IC79SL_signal ::Dep::IC79SL_data->signal;      }
    void IC79SL_bg     (double &result)   {
      result = Pipes::IC79SL_bg     ::Dep::IC79SL_data->bg;          }
    void IC79SL_nobs   (int    &result)   {
      result = Pipes::IC79SL_nobs   ::Dep::IC79SL_data->nobs;        }
    void IC79SL_loglike(double &result)   {
      result = Pipes::IC79SL_loglike::Dep::IC79SL_data->loglike;     }
    void IC79SL_bgloglike(double &result) {
      result = Pipes::IC79SL_bgloglike::Dep::IC79SL_data->bgloglike; }
    void IC79SL_pvalue (double &result)   {
      result = Pipes::IC79SL_pvalue ::Dep::IC79SL_data->pvalue;      }
    /// @}

    /// Composite IceCube 79-string likelihood function.
    void IC79_loglike(double &result)
    {
      using namespace Pipes::IC79_loglike;
      result = *Dep::IC79SL_loglike - *Dep::IC79SL_bgloglike +
               *Dep::IC79WL_loglike - *Dep::IC79WL_bgloglike +
               *Dep::IC79WH_loglike - *Dep::IC79WH_bgloglike;
      if (result > 0.0) result = 0.0;
    }

    /// Complete composite IceCube likelihood function.
    void IC_loglike(double &result)
    {
      using namespace Pipes::IC_loglike;
      result = *Dep::IC22_loglike   - *Dep::IC22_bgloglike   +
               *Dep::IC79SL_loglike - *Dep::IC79SL_bgloglike +
               *Dep::IC79WL_loglike - *Dep::IC79WL_bgloglike +
               *Dep::IC79WH_loglike - *Dep::IC79WH_bgloglike;
      if (result > 0.0) result = 0.0;
#ifdef DARKBIT_DEBUG
      cout << "IC likelihood: " << result << endl;
      cout << "IC79SL contribution: " << *Dep::IC79SL_loglike - *Dep::IC79SL_bgloglike << endl;
      cout << "IC79WL contribution: " << *Dep::IC79WL_loglike - *Dep::IC79WL_bgloglike << endl;
      cout << "IC79WH contribution: " << *Dep::IC79WH_loglike - *Dep::IC79WH_bgloglike << endl;
      cout << "IC22   contribution: " << *Dep::IC22_loglike   - *Dep::IC22_bgloglike   << endl;
#endif
    }

    /// Function to set Local Halo Parameters in DarkSUSY (DS5 only)
    void DarkSUSY5_PointInit_LocalHalo_func(bool &result)
    {
      using namespace Pipes::DarkSUSY5_PointInit_LocalHalo_func;

      LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;

      double rho0 = LocalHaloParameters.rho0;
      double rho0_eff = (*Dep::RD_fraction)*rho0;
      double vrot = LocalHaloParameters.vrot;
      double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
      double vesc = LocalHaloParameters.vesc;
      /// Option v_earth<double>: Keplerian velocity of the Earth around the Sun in km/s (default 29.78)
      double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");

      BEreq::dshmcom->rho0 = rho0;
      BEreq::dshmcom->v_sun = vrot;
      BEreq::dshmcom->v_earth = v_earth;
      BEreq::dshmcom->rhox = rho0_eff;

      BEreq::dshmframevelcom->v_obs = vrot;

      BEreq::dshmisodf->vd_3d = vd_3d;
      BEreq::dshmisodf->vgalesc = vesc;

      BEreq::dshmnoclue->vobs = vrot;

      logger() << LogTags::debug
               << "Updating DarkSUSY halo parameters:" << std::endl
               << "    rho0 [GeV/cm^3] = " << rho0 << std::endl
               << "    rho0_eff [GeV/cm^3] = " << rho0_eff << std::endl
               << "    v_sun [km/s]  = " << vrot<< std::endl
               << "    v_earth [km/s]  = " << v_earth << std::endl
               << "    v_obs [km/s]  = " << vrot << std::endl
               << "    vd_3d [km/s]  = " << vd_3d << std::endl
               << "    v_esc [km/s]  = " << vesc << EOM;

      result = true;

      return;
    }

    /// Function to set Local Halo Parameters in DarkSUSY (DS 6)
    void DarkSUSY_PointInit_LocalHalo_func(bool &result)
    {
      using namespace Pipes::DarkSUSY_PointInit_LocalHalo_func;

      LocalMaxwellianHalo LocalHaloParameters = *Dep::LocalHalo;

      double rho0 = LocalHaloParameters.rho0;
      double vrot = LocalHaloParameters.vrot;
      double vd_3d = sqrt(3./2.)*LocalHaloParameters.v0;
      double vesc = LocalHaloParameters.vesc;
      /// Option v_earth<double>: Keplerian velocity of the Earth around the Sun in km/s (default 29.78)
      double v_earth = runOptions->getValueOrDef<double>(29.78, "v_earth");

      BEreq::dshmcom->rho0 = rho0;
      BEreq::dshmcom->v_sun = vrot;
      BEreq::dshmcom->v_earth = v_earth;

      BEreq::dshmframevelcom->v_obs = vrot;

      BEreq::dshmisodf->vd_3d = vd_3d;
      BEreq::dshmisodf->vgalesc = vesc;

      BEreq::dshmnoclue->vobs = vrot;

      logger() << LogTags::debug
               << "Updating DarkSUSY halo parameters:" << std::endl
               << "    rho0 [GeV/cm^3] = " << rho0 << std::endl
               << "    v_sun [km/s]  = " << vrot<< std::endl
               << "    v_earth [km/s]  = " << v_earth << std::endl
               << "    v_obs [km/s]  = " << vrot << std::endl
               << "    vd_3d [km/s]  = " << vd_3d << std::endl
               << "    v_esc [km/s]  = " << vesc << EOM;

      result = true;

      return;
    }

  }
}

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