file src/SpecBit_DMEFT.cpp

[No description available] More…

Namespaces

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

Detailed Description

Author: The GAMBIT Collaboration

Date: 12:32PM on October 15, 2019

Implementation of SpecBit routines for DMEFT.

Authors (add name and date if you modify):

*** Automatically created by GUM ***


Source code

//   GAMBIT: Global and Modular BSM Inference Tool
//   *********************************************
///  \file
///
///  Implementation of SpecBit routines for 
///  DMEFT.
///
///  Authors (add name and date if you modify):    
///       *** Automatically created by GUM ***     
///                                                
///  \author The GAMBIT Collaboration             
///  \date 12:32PM on October 15, 2019
///                                                
///  ********************************************* 

#include "gambit/Elements/gambit_module_headers.hpp"
#include "gambit/Elements/spectrum.hpp"
#include "gambit/Utils/stream_overloads.hpp"
#include "gambit/Utils/util_macros.hpp"
#include "gambit/SpecBit/SpecBit_rollcall.hpp"
#include "gambit/SpecBit/SpecBit_helpers.hpp"
#include "gambit/SpecBit/QedQcdWrapper.hpp"
#include "gambit/Models/SimpleSpectra/DMEFTSimpleSpec.hpp"
#include "gambit/Models/SimpleSpectra/SMHiggsSimpleSpec.hpp"

namespace Gambit
{
  
  namespace SpecBit
  {
    using namespace LogTags;
    
    /// Get a simple wrapper for Spectrum object.
    void get_DMEFT_spectrum(Spectrum& result)
    {
      namespace myPipe = Pipes::get_DMEFT_spectrum;
      const SMInputs& sminputs = *myPipe::Dep::SMINPUTS;
      
      // Initialise model object 
      Models::DMEFTModel DMEFTmodel;
      
      // BSM parameters
      DMEFTmodel.DMEFT_Lambda = *myPipe::Param.at("Lambda");
      DMEFTmodel.DMEFT_C51 = *myPipe::Param.at("C51");
      DMEFTmodel.DMEFT_C52 = *myPipe::Param.at("C52");
      DMEFTmodel.DMEFT_C61 = *myPipe::Param.at("C61");
      DMEFTmodel.DMEFT_C62 = *myPipe::Param.at("C62");
      DMEFTmodel.DMEFT_C63 = *myPipe::Param.at("C63");
      DMEFTmodel.DMEFT_C64 = *myPipe::Param.at("C64");
      DMEFTmodel.DMEFT_C71 = *myPipe::Param.at("C71");
      DMEFTmodel.DMEFT_C72 = *myPipe::Param.at("C72");
      DMEFTmodel.DMEFT_C73 = *myPipe::Param.at("C73");
      DMEFTmodel.DMEFT_C74 = *myPipe::Param.at("C74");
      DMEFTmodel.DMEFT_C75 = *myPipe::Param.at("C75");
      DMEFTmodel.DMEFT_C76 = *myPipe::Param.at("C76");
      DMEFTmodel.DMEFT_C77 = *myPipe::Param.at("C77");
      DMEFTmodel.DMEFT_C78 = *myPipe::Param.at("C78");
      DMEFTmodel.DMEFT_C79 = *myPipe::Param.at("C79");
      DMEFTmodel.DMEFT_C710 = *myPipe::Param.at("C710");
      // Pole mass inputs (mh is a nuiisance parameter?)
      DMEFTmodel.DMEFT_chi_Pole_Mass = *myPipe::Param.at("mchi");
      DMEFTmodel.DMEFT_h0_1_Pole_Mass = *myPipe::Param.at("mH");
      // running top mass input, must be standard model msbar mt(mt)
      DMEFTmodel.mtrun = *myPipe::Param.at("mtrunIN");
      // Invalidate point if the EFT is violated for DM annihilation i.e. 2*m_DM > Lambda
      // Default: true
      if (myPipe::runOptions->getValueOrDef<bool>(true,"impose_EFT_validity"))
      {
        if (DMEFTmodel.DMEFT_Lambda < (2*DMEFTmodel.DMEFT_chi_Pole_Mass))
        {
          std::ostringstream msg;
          msg << "Parameter point [mchi, Lambda] = [" << DMEFTmodel.DMEFT_chi_Pole_Mass << " GeV, "
              << DMEFTmodel.DMEFT_Lambda << " GeV] does not satisfy the EFT.";
          invalid_point().raise(msg.str());
        }
      }
      else {}
      
      // quantities needed to fill container spectrum
      double alpha_em = 1.0 / sminputs.alphainv;
      double C = alpha_em * pi / (sminputs.GF * pow(2,0.5));
      double sinW2 = 0.5 - pow( 0.25 - C/pow(sminputs.mZ,2) , 0.5);
      double cosW2 = 0.5 + pow( 0.25 - C/pow(sminputs.mZ,2) , 0.5);
      double e = pow( 4*pi*( alpha_em ),0.5);
      double vev = 1. / sqrt(sqrt(2.)*sminputs.GF);
      
      // Gauge couplings
      DMEFTmodel.vev = vev;
      DMEFTmodel.g1 = e / sqrt(sinW2);
      DMEFTmodel.g2 = e / sqrt(cosW2);
      DMEFTmodel.g3 = pow( 4*pi*( sminputs.alphaS ),0.5);
      
      // Yukawas
      double sqrt2v = pow(2.0,0.5)/vev;
      DMEFTmodel.Yu[0][0] = sqrt2v * sminputs.mU;
      DMEFTmodel.Yu[1][1] = sqrt2v * sminputs.mCmC;
      // top quark is treated at one-loop with different running and pole mass
      DMEFTmodel.Yu[2][2] = sqrt2v * DMEFTmodel.mtrun;
      DMEFTmodel.Ye[0][0] = sqrt2v * sminputs.mE;
      DMEFTmodel.Ye[1][1] = sqrt2v * sminputs.mMu;
      DMEFTmodel.Ye[2][2] = sqrt2v * sminputs.mTau;
      DMEFTmodel.Yd[0][0] = sqrt2v * sminputs.mD;
      DMEFTmodel.Yd[1][1] = sqrt2v * sminputs.mS;
      DMEFTmodel.Yd[2][2] = sqrt2v * sminputs.mBmB;

      // Create a SubSpectrum object wrapper
      Models::DMEFTSimpleSpec spec(DMEFTmodel);
      
      // Retrieve any mass cuts
      static const Spectrum::mc_info mass_cut = myPipe::runOptions->getValueOrDef<Spectrum::mc_info>(Spectrum::mc_info(), "mass_cut");
      static const Spectrum::mr_info mass_ratio_cut = myPipe::runOptions->getValueOrDef<Spectrum::mr_info>(Spectrum::mr_info(), "mass_ratio_cut");

      // We have decided to calculate mT pole from the input running mt
      // using only one-loop QCD corrections
      // See footnote 9 https://link.springer.com/article/10.1007/JHEP11(2019)150

      // Approximate alpha_S(mtop) as alpha_S(mZ) because even a
      // one-loop correction here will lead to O(alpha_S^2) correction on mT
      // Via correspondence with the authors we checked this
      // approximation matches what is done in
      // https://link.springer.com/article/10.1007/JHEP11(2019)150
      // which we also validated numerically
      
      const double alpha_S_mtop = sminputs.alphaS;
      // Now extract pole mass from running mass
      const double mtop_MSBAR_mtop = DMEFTmodel.mtrun; // must be SM MSbar mt(mt)
      const double mtop_pole = mtop_MSBAR_mtop * (1. + 4. / 3.
                          * alpha_S_mtop / M_PI);

      // Make a local sminputs so we can change pole mT to the one we calculated
      SMInputs localsminputs = sminputs;
      localsminputs.mT = mtop_pole;
      // We don't supply a LE subspectrum here; an SMSimpleSpec will therefore be automatically created from 'localsminputs'
      result = Spectrum(spec,localsminputs,&myPipe::Param,mass_cut,mass_ratio_cut);

    }

    // Declaration: print spectrum out
    void fill_map_from_DMEFT_spectrum(std::map<std::string,double>&, const Spectrum&);
    
    void get_DMEFT_spectrum_as_map(std::map<std::string,double>& specmap)
    {
      namespace myPipe = Pipes::get_DMEFT_spectrum_as_map;
      const Spectrum& spec(*myPipe::Dep::DMEFT_spectrum);
      fill_map_from_DMEFT_spectrum(specmap, spec);
    }
    
    void fill_map_from_DMEFT_spectrum(std::map<std::string, double>& specmap, const Spectrum& spec)
    {
      /// Use SpectrumContents routines to automate
      static const SpectrumContents::DMEFT contents;
      static const std::vector<SpectrumParameter> required_parameters = contents.all_parameters();
      
      for(std::vector<SpectrumParameter>::const_iterator it = required_parameters.begin(); it != required_parameters.end(); ++it)
      {
        const Par::Tags        tag   = it->tag();
        const std::string      name  = it->name();
        const std::vector<int> shape = it->shape();
        
        // Scalar case
        if(shape.size()==1 and shape[0]==1)
        {
          std::ostringstream label;
          label << name <<" "<< Par::toString.at(tag);
          specmap[label.str()] = spec.get_HE().get(tag,name);
        }
        // Vector case
        else if(shape.size()==1 and shape[0]>1)
        {
          for(int i = 1; i<=shape[0]; ++i)
          {
            std::ostringstream label;
            label << name <<"_"<<i<<" "<< Par::toString.at(tag);
            specmap[label.str()] = spec.get_HE().get(tag,name,i);
          }
        }
        // Matrix case
        else if(shape.size()==2)
        {
          for(int i = 1; i<=shape[0]; ++i)
          {
            for(int j = 1; j<=shape[0]; ++j)
            {
              std::ostringstream label;
              label << name <<"_("<<i<<","<<j<<") "<<Par::toString.at(tag);
              specmap[label.str()] = spec.get_HE().get(tag,name,i,j);
            }
          }
        }
        // Deal with all other cases
        else
        {
          // ERROR
          std::ostringstream errmsg;
          errmsg << "Invalid parameter received while converting DMEFT_spectrum to map of strings!";
          errmsg << "Problematic parameter was: "<< tag <<", " << name << ", shape="<< shape;
          utils_error().forced_throw(LOCAL_INFO,errmsg.str());
        }

        // Include the pole mass of the top (which is a derived parameter and not in SpectrumContents::DMEFT)
        std::ostringstream label;
        label << "t" <<" "<< Par::toString.at(Par::Pole_Mass);
        specmap[label.str()] = spec.get(Par::Pole_Mass,"t");
      }
    }
    
  }
  
}

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