file src/ColliderBit_Higgs.cpp
[No description available] More…
Namespaces
| Name |
|---|
| Gambit TODO: see if we can use this one: |
| Gambit::ColliderBit |
Detailed Description
Author:
- Chris Rogan (crogan@cern.ch)
- Pat Scott (p.scott@imperial.ac.uk)
- James McKay (j.mckay14@imperial.ac.uk)
- Sanjay Bloor (sanjay.bloor12@imperial.ac.uk)
- Tomas Gonzalo (tomas.gonzalo@monash.edu)
Date:
- 2014 Aug
- 2015 May
- 2015 Jul
- 2016 Sep
- 2016 Sep
- 2020 Mar
- 2020 Mar
Functions of ColliderBit that deal exclusively with Higgs physics some functions were originally in CollderBit.cpp
Authors (add name and date if you modify):
Source code
// GAMBIT: Global and Modular BSM Inference Tool
// *********************************************
/// \file
///
/// Functions of ColliderBit that deal exclusively with Higgs physics
/// some functions were originally in CollderBit.cpp
///
/// *********************************************
///
/// Authors (add name and date if you modify):
///
/// \author Chris Rogan
/// (crogan@cern.ch)
/// \date 2014 Aug
/// \date 2015 May
///
/// \author Pat Scott
/// (p.scott@imperial.ac.uk)
/// \date 2015 Jul
/// \date 2016 Sep
///
/// \author James McKay
/// (j.mckay14@imperial.ac.uk)
/// \date 2016 Sep
///
/// \author Sanjay Bloor
/// (sanjay.bloor12@imperial.ac.uk)
/// \date 2020 Mar
///
/// \author Tomas Gonzalo
/// (tomas.gonzalo@monash.edu)
/// \date 2020 Mar
///
/// *********************************************
#include <cmath>
#include <string>
#include <iostream>
#include <fstream>
#include <memory>
#include <numeric>
#include <sstream>
#include <vector>
#include "gambit/Elements/gambit_module_headers.hpp"
#include "gambit/Utils/util_types.hpp"
#include "gambit/ColliderBit/ColliderBit_rollcall.hpp"
//#define COLLIDERBIT_DEBUG
namespace Gambit
{
namespace ColliderBit
{
/// Helper function to set HiggsBounds/Signals parameters cross-section ratios from a GAMBIT HiggsCouplingsTable
void set_CS(hb_ModelParameters &result, const HiggsCouplingsTable& couplings, int n_neutral_higgses)
{
for(int i = 0; i < n_neutral_higgses; i++)
{
result.CS_bg_hjb_ratio[i] = couplings.C_bb2[i];
result.CS_bb_hj_ratio[i] = couplings.C_bb2[i];
result.CS_lep_bbhj_ratio[i] = couplings.C_bb2[i];
result.CS_lep_tautauhj_ratio[i] = couplings.C_tautau2[i];
result.CS_lep_hjZ_ratio[i] = couplings.C_ZZ2[i];
result.CS_gg_hjZ_ratio[i] = 0.;
result.CS_dd_hjZ_ratio[i] = couplings.C_ZZ2[i];
result.CS_uu_hjZ_ratio[i] = couplings.C_ZZ2[i];
result.CS_ss_hjZ_ratio[i] = couplings.C_ZZ2[i];
result.CS_cc_hjZ_ratio[i] = couplings.C_ZZ2[i];
result.CS_bb_hjZ_ratio[i] = couplings.C_ZZ2[i];
result.CS_ud_hjWp_ratio[i] = couplings.C_WW2[i];
result.CS_cs_hjWp_ratio[i] = couplings.C_WW2[i];
result.CS_ud_hjWm_ratio[i] = couplings.C_WW2[i];
result.CS_cs_hjWm_ratio[i] = couplings.C_WW2[i];
result.CS_tev_vbf_ratio[i] = couplings.C_WW2[i];
result.CS_lhc7_vbf_ratio[i] = couplings.C_WW2[i];
result.CS_lhc8_vbf_ratio[i] = couplings.C_WW2[i];
result.CS_gg_hj_ratio[i] = couplings.C_gg2[i];
result.CS_tev_tthj_ratio[i] = couplings.C_tt2[i];
result.CS_lhc7_tthj_ratio[i] = couplings.C_tt2[i];
result.CS_lhc8_tthj_ratio[i] = couplings.C_tt2[i];
for(int j = 0; j < n_neutral_higgses; j++)
{
result.CS_lep_hjhi_ratio[i][j] = couplings.C_hiZ2[i][j];
}
}
// LEP H+ H- x-section ratio
result.CS_lep_HpjHmi_ratio[0] = 1.;
}
/// Helper function for populating a HiggsBounds/Signals ModelParameters object for SM-like Higgs.
void set_SMLikeHiggs_ModelParameters(const SubSpectrum& spec, const HiggsCouplingsTable& couplings, hb_ModelParameters &result)
{
// Retrieve the decays
const DecayTable::Entry& decays = couplings.get_neutral_decays(0);
// Set the CP of the lone higgs
result.CP[0] = couplings.CP[0];
// Set h mass and uncertainty
result.Mh[0] = spec.get(Par::Pole_Mass,25,0);
bool has_high_err = spec.has(Par::Pole_Mass_1srd_high, 25, 0);
bool has_low_err = spec.has(Par::Pole_Mass_1srd_low, 25, 0);
if (has_high_err and has_low_err)
{
double upper = spec.get(Par::Pole_Mass_1srd_high, 25, 0);
double lower = spec.get(Par::Pole_Mass_1srd_low, 25, 0);
result.deltaMh[0] = result.Mh[0] * std::max(upper,lower);
}
else
{
result.deltaMh[0] = 0.;
}
// Set the total h width
result.hGammaTot[0] = decays.width_in_GeV;
// Set the branching fractions
result.BR_hjss[0] = decays.BF("s", "sbar");
result.BR_hjcc[0] = decays.BF("c", "cbar");
result.BR_hjbb[0] = decays.BF("b", "bbar");
result.BR_hjmumu[0] = decays.BF("mu+", "mu-");
result.BR_hjtautau[0] = decays.BF("tau+", "tau-");
result.BR_hjWW[0] = decays.BF("W+", "W-");
result.BR_hjZZ[0] = decays.BF("Z0", "Z0");
result.BR_hjZga[0] = decays.BF("gamma", "Z0");
result.BR_hjgaga[0] = decays.BF("gamma", "gamma");
result.BR_hjgg[0] = decays.BF("g", "g");
// Add the invisibles
result.BR_hjinvisible[0] = 0.;
for (std::vector<std::pair<str,str>>::const_iterator it = couplings.invisibles.begin(); it != couplings.invisibles.end(); ++it)
{
result.BR_hjinvisible[0] += decays.BF(it->first, it->second);
}
// Retrieve cross-section ratios from the HiggsCouplingsTable
set_CS(result, couplings, 1);
// Zero all heavy neutral higgs masses, widths and effective couplings
for(int i = 1; i < 3; i++)
{
result.Mh[i] = 0.;
result.deltaMh[i] = 0.;
result.hGammaTot[i] = 0.;
result.CP[i] = 0.;
result.BR_hjss[i] = 0.;
result.BR_hjcc[i] = 0.;
result.BR_hjbb[i] = 0.;
result.BR_hjmumu[i] = 0.;
result.BR_hjtautau[i] = 0.;
result.BR_hjWW[i] = 0.;
result.BR_hjZZ[i] = 0.;
result.BR_hjZga[i] = 0.;
result.BR_hjgaga[i] = 0.;
result.BR_hjgg[i] = 0.;
result.BR_hjinvisible[i] = 0.;
result.CS_lep_hjZ_ratio[i] = 0.;
result.CS_lep_bbhj_ratio[i] = 0.;
result.CS_lep_tautauhj_ratio[i] = 0.;
result.CS_gg_hj_ratio[i] = 0.;
result.CS_bb_hj_ratio[i] = 0.;
result.CS_bg_hjb_ratio[i] = 0.;
result.CS_ud_hjWp_ratio[i] = 0.;
result.CS_cs_hjWp_ratio[i] = 0.;
result.CS_ud_hjWm_ratio[i] = 0.;
result.CS_cs_hjWm_ratio[i] = 0.;
result.CS_gg_hjZ_ratio[i] = 0.;
result.CS_dd_hjZ_ratio[i] = 0.;
result.CS_uu_hjZ_ratio[i] = 0.;
result.CS_ss_hjZ_ratio[i] = 0.;
result.CS_cc_hjZ_ratio[i] = 0.;
result.CS_bb_hjZ_ratio[i] = 0.;
result.CS_tev_vbf_ratio[i] = 0.;
result.CS_tev_tthj_ratio[i] = 0.;
result.CS_lhc7_vbf_ratio[i] = 0.;
result.CS_lhc7_tthj_ratio[i] = 0.;
result.CS_lhc8_vbf_ratio[i] = 0.;
result.CS_lhc8_tthj_ratio[i] = 0.;
for(int j = 0; j < 3; j++) result.BR_hjhihi[i][j] = 0.;
for(int j = 0; j < 3; j++) result.CS_lep_hjhi_ratio[i][j] = 0.;
}
// Zero all H+ masses, widths and effective couplings
result.MHplus[0] = 0.;
result.deltaMHplus[0] = 0.;
result.HpGammaTot[0] = 0.;
result.BR_tWpb = 0.;
result.BR_tHpjb[0] = 0.;
result.BR_Hpjcs[0] = 0.;
result.BR_Hpjcb[0] = 0.;
result.BR_Hptaunu[0] = 0.;
result.CS_lep_HpjHmi_ratio[0] = 0.;
}
/// SM-like (SM + possible invisibles) Higgs model parameters for HiggsBounds/Signals
void SMLikeHiggs_ModelParameters(hb_ModelParameters &result)
{
using namespace Pipes::SMLikeHiggs_ModelParameters;
dep_bucket<Spectrum>* spectrum_dependency = nullptr;
if (ModelInUse("ScalarSingletDM_Z2") or ModelInUse("ScalarSingletDM_Z2_running")) spectrum_dependency = &Dep::ScalarSingletDM_Z2_spectrum;
else if (ModelInUse("ScalarSingletDM_Z3") or ModelInUse("ScalarSingletDM_Z3_running")) spectrum_dependency = &Dep::ScalarSingletDM_Z3_spectrum;
else if (ModelInUse("StandardModel_Higgs") or ModelInUse("StandardModel_Higgs_running")) spectrum_dependency = &Dep::SM_spectrum;
else ColliderBit_error().raise(LOCAL_INFO, "No valid model for SMLikeHiggs_ModelParameters.");
const SubSpectrum& spec = (*spectrum_dependency)->get_HE();
set_SMLikeHiggs_ModelParameters(spec, *Dep::Higgs_Couplings, result);
}
/// MSSM-like (MSSM + NMSSM + ...) Higgs model parameters for HiggsBounds/Signals
void MSSMLikeHiggs_ModelParameters(hb_ModelParameters &result)
{
using namespace Pipes::MSSMLikeHiggs_ModelParameters;
// Set up neutral Higgses
int n_neutral_higgses = Dep::Higgs_Couplings->get_n_neutral_higgs();
// Set the CP of the Higgs states.
for (int i = 0; i < n_neutral_higgses; i++) result.CP[i] = Dep::Higgs_Couplings->CP[i];
// Retrieve higgs partial widths
const std::vector<const DecayTable::Entry*>& h0_widths = Dep::Higgs_Couplings->get_neutral_decays_array();
const DecayTable::Entry& H_plus_widths = Dep::Higgs_Couplings->get_charged_decays(0);
const DecayTable::Entry& t_widths = Dep::Higgs_Couplings->get_t_decays();
// Pick the correct spectrum and specify the Higgses
dep_bucket<Spectrum>* spectrum_dependency = nullptr;
std::vector<str> Higgses;
if (ModelInUse("MSSM63atMGUT") or ModelInUse("MSSM63atQ"))
{
spectrum_dependency = &Dep::MSSM_spectrum;
Higgses = initVector<str>("h0_1", "h0_2", "A0");
}
else ColliderBit_error().raise(LOCAL_INFO, "No valid model for MSSMLikeHiggs_ModelParameters.");
const SubSpectrum& spec = (*spectrum_dependency)->get_HE();
static const std::vector<str> sHneut(Higgses);
// Neutral higgs masses and errors
for(int i = 0; i < n_neutral_higgses; i++)
{
result.Mh[i] = spec.get(Par::Pole_Mass,sHneut[i]);
double upper = spec.get(Par::Pole_Mass_1srd_high,sHneut[i]);
double lower = spec.get(Par::Pole_Mass_1srd_low,sHneut[i]);
result.deltaMh[i] = result.Mh[i] * std::max(upper,lower);
}
// Loop over all neutral Higgses, setting their branching fractions and total widths.
for(int i = 0; i < n_neutral_higgses; i++)
{
result.hGammaTot[i] = h0_widths[i]->width_in_GeV;
result.BR_hjss[i] = h0_widths[i]->BF("s", "sbar");
result.BR_hjcc[i] = h0_widths[i]->BF("c", "cbar");
result.BR_hjbb[i] = h0_widths[i]->BF("b", "bbar");
result.BR_hjmumu[i] = h0_widths[i]->BF("mu+", "mu-");
result.BR_hjtautau[i] = h0_widths[i]->BF("tau+", "tau-");
result.BR_hjWW[i] = h0_widths[i]->has_channel("W+", "W-") ? h0_widths[i]->BF("W+", "W-") : 0.0;
result.BR_hjZZ[i] = h0_widths[i]->has_channel("Z0", "Z0") ? h0_widths[i]->BF("Z0", "Z0") : 0.0;
result.BR_hjZga[i] = h0_widths[i]->has_channel("gamma", "Z0") ? h0_widths[i]->BF("gamma", "Z0") : 0.0;
result.BR_hjgaga[i] = h0_widths[i]->BF("gamma", "gamma");
result.BR_hjgg[i] = h0_widths[i]->BF("g", "g");
// Do decays to invisibles
result.BR_hjinvisible[i] = 0.;
for (std::vector<std::pair<str,str>>::const_iterator it = Dep::Higgs_Couplings->invisibles.begin(); it != Dep::Higgs_Couplings->invisibles.end(); ++it)
{
result.BR_hjinvisible[i] += h0_widths[i]->BF(it->first, it->second);
}
// Do decays to other neutral higgses
for (int j = 0; j < n_neutral_higgses; j++)
{
if (2.*result.Mh[j] < result.Mh[i] and h0_widths[i]->has_channel(sHneut[j],sHneut[j]))
{
result.BR_hjhihi[i][j] = h0_widths[i]->BF(sHneut[j],sHneut[j]);
}
else
{
result.BR_hjhihi[i][j] = 0.;
}
}
}
// Charged higgs masses and errors
result.MHplus[0] = spec.get(Par::Pole_Mass,"H+");
double upper = spec.get(Par::Pole_Mass_1srd_high,"H+");
double lower = spec.get(Par::Pole_Mass_1srd_low,"H+");
result.deltaMHplus[0] = result.MHplus[0] * std::max(upper,lower);
// Set charged Higgs branching fractions and total width.
result.HpGammaTot[0] = H_plus_widths.width_in_GeV;
result.BR_Hpjcs[0] = H_plus_widths.BF("c", "sbar");
result.BR_Hpjcb[0] = H_plus_widths.BF("c", "bbar");
result.BR_Hptaunu[0] = H_plus_widths.BF("tau+", "nu_tau");
// Set top branching fractions
result.BR_tWpb = t_widths.BF("W+", "b");
result.BR_tHpjb[0] = t_widths.has_channel("H+", "b") ? t_widths.BF("H+", "b") : 0.0;
// Retrieve cross-section ratios from the HiggsCouplingsTable
set_CS(result, *Dep::Higgs_Couplings, n_neutral_higgses);
}
/// Get a LEP chisq from HiggsBounds
void calc_HB_LEP_LogLike(double &result)
{
using namespace Pipes::calc_HB_LEP_LogLike;
hb_ModelParameters ModelParam = *Dep::HB_ModelParameters;
Farray<double, 1,3, 1,3> CS_lep_hjhi_ratio;
Farray<double, 1,3, 1,3> BR_hjhihi;
// Transpose to get around Fortran matrix types
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
}
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
}
BEreq::HiggsBounds_neutral_input_part(&ModelParam.Mh[0], &ModelParam.hGammaTot[0], &ModelParam.CP[0],
&ModelParam.CS_lep_hjZ_ratio[0], &ModelParam.CS_lep_bbhj_ratio[0],
&ModelParam.CS_lep_tautauhj_ratio[0], &ModelParam.CS_lep_hjhi_ratio[0][0],
&ModelParam.CS_gg_hj_ratio[0], &ModelParam.CS_bb_hj_ratio[0],
&ModelParam.CS_bg_hjb_ratio[0], &ModelParam.CS_ud_hjWp_ratio[0],
&ModelParam.CS_cs_hjWp_ratio[0], &ModelParam.CS_ud_hjWm_ratio[0],
&ModelParam.CS_cs_hjWm_ratio[0], &ModelParam.CS_gg_hjZ_ratio[0],
&ModelParam.CS_dd_hjZ_ratio[0], &ModelParam.CS_uu_hjZ_ratio[0],
&ModelParam.CS_ss_hjZ_ratio[0], &ModelParam.CS_cc_hjZ_ratio[0],
&ModelParam.CS_bb_hjZ_ratio[0], &ModelParam.CS_tev_vbf_ratio[0],
&ModelParam.CS_tev_tthj_ratio[0], &ModelParam.CS_lhc7_vbf_ratio[0],
&ModelParam.CS_lhc7_tthj_ratio[0], &ModelParam.CS_lhc8_vbf_ratio[0],
&ModelParam.CS_lhc8_tthj_ratio[0], &ModelParam.BR_hjss[0],
&ModelParam.BR_hjcc[0], &ModelParam.BR_hjbb[0],
&ModelParam.BR_hjmumu[0], &ModelParam.BR_hjtautau[0],
&ModelParam.BR_hjWW[0], &ModelParam.BR_hjZZ[0],
&ModelParam.BR_hjZga[0], &ModelParam.BR_hjgaga[0],
&ModelParam.BR_hjgg[0], &ModelParam.BR_hjinvisible[0], &ModelParam.BR_hjhihi[0][0]);
// Transpose it back
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
}
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
}
BEreq::HiggsBounds_charged_input(&ModelParam.MHplus[0], &ModelParam.HpGammaTot[0], &ModelParam.CS_lep_HpjHmi_ratio[0],
&ModelParam.BR_tWpb, &ModelParam.BR_tHpjb[0], &ModelParam.BR_Hpjcs[0],
&ModelParam.BR_Hpjcb[0], &ModelParam.BR_Hptaunu[0]);
BEreq::HiggsBounds_set_mass_uncertainties(&ModelParam.deltaMh[0],&ModelParam.deltaMHplus[0]);
// run Higgs bounds 'classic'
double obsratio;
int HBresult, chan, ncombined;
BEreq::run_HiggsBounds_classic(HBresult,chan,obsratio,ncombined);
// extract the LEP chisq
double chisq_withouttheory,chisq_withtheory;
int chan2;
double theor_unc = 1.5; // theory uncertainty
BEreq::HB_calc_stats(theor_unc,chisq_withouttheory,chisq_withtheory,chan2);
// Catch HiggsBound's error value, chisq = -999
if( fabs(chisq_withouttheory - (-999.)) < 1e-6)
{
ColliderBit_warning().raise(LOCAL_INFO, "Got chisq=-999 from HB_calc_stats in HiggsBounds, indicating a cross-section outside tabulated range. Will use chisq=0.");
chisq_withouttheory = 0.0;
}
result = -0.5*chisq_withouttheory;
}
/// Get an LHC chisq from HiggsSignals
void calc_HS_LHC_LogLike(double &result)
{
using namespace Pipes::calc_HS_LHC_LogLike;
hb_ModelParameters ModelParam = *Dep::HB_ModelParameters;
Farray<double, 1,3, 1,3> CS_lep_hjhi_ratio;
Farray<double, 1,3, 1,3> BR_hjhihi;
// Transpose to get around Fortran matrix types
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
}
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
}
BEreq::HiggsBounds_neutral_input_part_HS(&ModelParam.Mh[0], &ModelParam.hGammaTot[0], &ModelParam.CP[0],
&ModelParam.CS_lep_hjZ_ratio[0], &ModelParam.CS_lep_bbhj_ratio[0],
&ModelParam.CS_lep_tautauhj_ratio[0], &ModelParam.CS_lep_hjhi_ratio[0][0],
&ModelParam.CS_gg_hj_ratio[0], &ModelParam.CS_bb_hj_ratio[0],
&ModelParam.CS_bg_hjb_ratio[0], &ModelParam.CS_ud_hjWp_ratio[0],
&ModelParam.CS_cs_hjWp_ratio[0], &ModelParam.CS_ud_hjWm_ratio[0],
&ModelParam.CS_cs_hjWm_ratio[0], &ModelParam.CS_gg_hjZ_ratio[0],
&ModelParam.CS_dd_hjZ_ratio[0], &ModelParam.CS_uu_hjZ_ratio[0],
&ModelParam.CS_ss_hjZ_ratio[0], &ModelParam.CS_cc_hjZ_ratio[0],
&ModelParam.CS_bb_hjZ_ratio[0], &ModelParam.CS_tev_vbf_ratio[0],
&ModelParam.CS_tev_tthj_ratio[0], &ModelParam.CS_lhc7_vbf_ratio[0],
&ModelParam.CS_lhc7_tthj_ratio[0], &ModelParam.CS_lhc8_vbf_ratio[0],
&ModelParam.CS_lhc8_tthj_ratio[0], &ModelParam.BR_hjss[0],
&ModelParam.BR_hjcc[0], &ModelParam.BR_hjbb[0],
&ModelParam.BR_hjmumu[0], &ModelParam.BR_hjtautau[0],
&ModelParam.BR_hjWW[0], &ModelParam.BR_hjZZ[0],
&ModelParam.BR_hjZga[0], &ModelParam.BR_hjgaga[0],
&ModelParam.BR_hjgg[0], &ModelParam.BR_hjinvisible[0], &ModelParam.BR_hjhihi[0][0]);
// Transpose it back
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
}
for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
{
ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
}
BEreq::HiggsBounds_charged_input_HS(&ModelParam.MHplus[0], &ModelParam.HpGammaTot[0], &ModelParam.CS_lep_HpjHmi_ratio[0],
&ModelParam.BR_tWpb, &ModelParam.BR_tHpjb[0], &ModelParam.BR_Hpjcs[0],
&ModelParam.BR_Hpjcb[0], &ModelParam.BR_Hptaunu[0]);
BEreq::HiggsSignals_neutral_input_MassUncertainty(&ModelParam.deltaMh[0]);
// add uncertainties to cross-sections and branching ratios
// double dCS[5] = {0.,0.,0.,0.,0.};
// double dBR[5] = {0.,0.,0.,0.,0.};
// BEreq::setup_rate_uncertainties(dCS,dBR);
// run HiggsSignals
int mode = 1; // 1- peak-centered chi2 method (recommended)
double csqmu, csqmh, csqtot, Pvalue;
int nobs;
BEreq::run_HiggsSignals(mode, csqmu, csqmh, csqtot, nobs, Pvalue);
result = -0.5*csqtot;
#ifdef COLLIDERBIT_DEBUG
std::ofstream f;
f.open ("HB_ModelParameters_contents.dat");
f<<"LHC log-likleihood";
for (int i = 0; i < 3; i++) f<<
" higgs index" <<
" "<<i<<":CP"<<
" "<<i<<":Mh"<<
" "<<i<<":hGammaTot"<<
" "<<i<<":CS_lep_hjZ_ratio"<<
" "<<i<<":CS_tev_vbf_ratio"<<
" "<<i<<":CS_lep_bbhj_ratio"<<
" "<<i<<":CS_lep_tautauhj_ratio"<<
" "<<i<<":CS_gg_hj_ratio"<<
" "<<i<<":CS_tev_tthj_ratio"<<
" "<<i<<":CS_lhc7_tthj_ratio"<<
" "<<i<<":CS_lhc8_tthj_ratio"<<
" "<<i<<":CS_lep_hjhi_ratio[0]"<<
" "<<i<<":CS_lep_hjhi_ratio[1]"<<
" "<<i<<":CS_lep_hjhi_ratio[2]"<<
" "<<i<<":BR_ss"<<
" "<<i<<":BR_cc"<<
" "<<i<<":BR_bb"<<
" "<<i<<":BR_mumu"<<
" "<<i<<":BR_tautau"<<
" "<<i<<":BR_WW"<<
" "<<i<<":BR_ZZ"<<
" "<<i<<":BR_Zga"<<
" "<<i<<":BR_gamgam"<<
" "<<i<<":BR_gg"<<
" "<<i<<":BR_invisible"<<
" "<<i<<":BR_hihi[0]"<<
" "<<i<<":BR_hihi[1]"<<
" "<<i<<":BR_hihi[2]";
f<<
" higgs index" <<
" "<<4<<"MHplus"<<
" "<<4<<":HpGammaTot"<<
" "<<4<<":CS_lep_HpjHmi_ratio"<<
" "<<4<<":BR_H+->cs"<<
" "<<4<<":BR_H+->cb"<<
" "<<4<<":BR_H+->taunu"<<
" "<<4<<":BR_t->W+b"<<
" "<<4<<":BR_t->H+b";
f << endl << std::setw(18) << result;
const int w = 24;
for (int i = 0; i < 3; i++)
{
f << std::setw(w) << i << std::setw(w) <<
ModelParam.CP[i] << std::setw(w) <<
ModelParam.Mh[i] << std::setw(w) <<
ModelParam.hGammaTot[i] << std::setw(w) <<
ModelParam.CS_lep_hjZ_ratio[i] << std::setw(w) <<
ModelParam.CS_tev_vbf_ratio[i] << std::setw(w) <<
ModelParam.CS_lep_bbhj_ratio[i] << std::setw(w) <<
ModelParam.CS_lep_tautauhj_ratio[i] << std::setw(w) <<
ModelParam.CS_gg_hj_ratio[i] << std::setw(w) <<
ModelParam.CS_tev_tthj_ratio[i] << std::setw(w) <<
ModelParam.CS_lhc7_tthj_ratio[i] << std::setw(w) <<
ModelParam.CS_lhc8_tthj_ratio[i];
for (int j = 0; j < 3; j++) f << std::setw(w) << ModelParam.CS_lep_hjhi_ratio[i][j];
f << std::setw(w) <<
ModelParam.BR_hjss[i] << std::setw(w) <<
ModelParam.BR_hjcc[i] << std::setw(w) <<
ModelParam.BR_hjbb[i] << std::setw(w) <<
ModelParam.BR_hjmumu[i] << std::setw(w) <<
ModelParam.BR_hjtautau[i] << std::setw(w) <<
ModelParam.BR_hjWW[i] << std::setw(w) <<
ModelParam.BR_hjZZ[i] << std::setw(w) <<
ModelParam.BR_hjZga[i] << std::setw(w) <<
ModelParam.BR_hjgaga[i] << std::setw(w) <<
ModelParam.BR_hjgg[i] << std::setw(w) <<
ModelParam.BR_hjinvisible[i];
for (int j = 0; j < 3; j++) f << std::setw(w) << ModelParam.BR_hjhihi[i][j];
}
f << std::setw(w) << 4 << std::setw(w) <<
ModelParam.MHplus[0] << std::setw(w) <<
ModelParam.HpGammaTot[0] << std::setw(w) <<
ModelParam.CS_lep_HpjHmi_ratio[0] << std::setw(w) <<
ModelParam.BR_Hpjcs[0] << std::setw(w) <<
ModelParam.BR_Hpjcb[0] << std::setw(w) <<
ModelParam.BR_Hptaunu[0] << std::setw(w) <<
ModelParam.BR_tWpb << std::setw(w) <<
ModelParam.BR_tHpjb[0];
f.close();
#endif
}
/// Higgs production cross-sections from FeynHiggs.
void FeynHiggs_HiggsProd(fh_HiggsProd_container &result)
{
using namespace Pipes::FeynHiggs_HiggsProd;
Farray<fh_real, 1,52> prodxs;
fh_HiggsProd_container HiggsProd;
int error;
fh_real sqrts;
// Tevatron
sqrts = 2.;
error = 1;
BEreq::FHHiggsProd(error, sqrts, prodxs);
if (error != 0)
{
std::ostringstream err;
err << "BEreq::FHHiggsProd raised error flag for Tevatron: " << error << ".";
invalid_point().raise(err.str());
}
for(int i = 0; i < 52; i++) HiggsProd.prodxs_Tev[i] = prodxs(i+1);
// LHC7
sqrts = 7.;
error = 1;
BEreq::FHHiggsProd(error, sqrts, prodxs);
if (error != 0)
{
std::ostringstream err;
err << "BEreq::FHHiggsProd raised error flag for LHC7: " << error << ".";
invalid_point().raise(err.str());
}
for(int i = 0; i < 52; i++) HiggsProd.prodxs_LHC7[i] = prodxs(i+1);
// LHC8
sqrts = 8.;
error = 1;
BEreq::FHHiggsProd(error, sqrts, prodxs);
if (error != 0)
{
std::ostringstream err;
err << "BEreq::FHHiggsProd raised error flag for LHC8: " << error << ".";
invalid_point().raise(err.str());
}
for(int i = 0; i < 52; i++) HiggsProd.prodxs_LHC8[i] = prodxs(i+1);
// The ttbar production cross-sections for the (BSM,SM) model can be found at (prodxs_X[h+27], prodxs_X[h+30]),
// where h is the higgs index (0 = h0_1, 1 = h0_2, 2 = A0) and X is one of Tev, LHC7 or LHC8.
result = HiggsProd;
}
}
}
Updated on 2023-06-26 at 21:36:56 +0000