gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
ColliderBit_Higgs.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
35 
36 #include <cmath>
37 #include <string>
38 #include <iostream>
39 #include <fstream>
40 #include <memory>
41 #include <numeric>
42 #include <sstream>
43 #include <vector>
44 
48 
49 //#define COLLIDERBIT_DEBUG
50 
51 
52 namespace Gambit
53 {
54 
55  namespace ColliderBit
56  {
57 
59  void set_CS(hb_ModelParameters &result, const HiggsCouplingsTable& couplings, int n_neutral_higgses)
60  {
61  for(int i = 0; i < n_neutral_higgses; i++)
62  {
63  result.CS_bg_hjb_ratio[i] = couplings.C_bb2[i];
64  result.CS_bb_hj_ratio[i] = couplings.C_bb2[i];
65  result.CS_lep_bbhj_ratio[i] = couplings.C_bb2[i];
66 
67  result.CS_lep_tautauhj_ratio[i] = couplings.C_tautau2[i];
68 
69  result.CS_lep_hjZ_ratio[i] = couplings.C_ZZ2[i];
70  result.CS_gg_hjZ_ratio[i] = 0.;
71  result.CS_dd_hjZ_ratio[i] = couplings.C_ZZ2[i];
72  result.CS_uu_hjZ_ratio[i] = couplings.C_ZZ2[i];
73  result.CS_ss_hjZ_ratio[i] = couplings.C_ZZ2[i];
74  result.CS_cc_hjZ_ratio[i] = couplings.C_ZZ2[i];
75  result.CS_bb_hjZ_ratio[i] = couplings.C_ZZ2[i];
76 
77  result.CS_ud_hjWp_ratio[i] = couplings.C_WW2[i];
78  result.CS_cs_hjWp_ratio[i] = couplings.C_WW2[i];
79  result.CS_ud_hjWm_ratio[i] = couplings.C_WW2[i];
80  result.CS_cs_hjWm_ratio[i] = couplings.C_WW2[i];
81 
82  result.CS_tev_vbf_ratio[i] = couplings.C_WW2[i];
83  result.CS_lhc7_vbf_ratio[i] = couplings.C_WW2[i];
84  result.CS_lhc8_vbf_ratio[i] = couplings.C_WW2[i];
85 
86  result.CS_gg_hj_ratio[i] = couplings.C_gg2[i];
87 
88  result.CS_tev_tthj_ratio[i] = couplings.C_tt2[i];
89  result.CS_lhc7_tthj_ratio[i] = couplings.C_tt2[i];
90  result.CS_lhc8_tthj_ratio[i] = couplings.C_tt2[i];
91 
92  for(int j = 0; j < n_neutral_higgses; j++)
93  {
94  result.CS_lep_hjhi_ratio[i][j] = couplings.C_hiZ2[i][j];
95  }
96  }
97  // LEP H+ H- x-section ratio
98  result.CS_lep_HpjHmi_ratio[0] = 1.;
99  }
100 
102  void set_SMLikeHiggs_ModelParameters(const SubSpectrum& spec, const HiggsCouplingsTable& couplings, hb_ModelParameters &result)
103  {
104  // Retrieve the decays
105  const DecayTable::Entry& decays = couplings.get_neutral_decays(0);
106 
107  // Set the CP of the lone higgs
108  result.CP[0] = couplings.CP[0];
109 
110  // Set h mass and uncertainty
111  result.Mh[0] = spec.get(Par::Pole_Mass,25,0);
112  bool has_high_err = spec.has(Par::Pole_Mass_1srd_high, 25, 0);
113  bool has_low_err = spec.has(Par::Pole_Mass_1srd_low, 25, 0);
114  if (has_high_err and has_low_err)
115  {
116  double upper = spec.get(Par::Pole_Mass_1srd_high, 25, 0);
117  double lower = spec.get(Par::Pole_Mass_1srd_low, 25, 0);
118  result.deltaMh[0] = result.Mh[0] * std::max(upper,lower);
119  }
120  else
121  {
122  result.deltaMh[0] = 0.;
123  }
124 
125  // Set the total h width
126  result.hGammaTot[0] = decays.width_in_GeV;
127 
128  // Set the branching fractions
129  result.BR_hjss[0] = decays.BF("s", "sbar");
130  result.BR_hjcc[0] = decays.BF("c", "cbar");
131  result.BR_hjbb[0] = decays.BF("b", "bbar");
132  result.BR_hjmumu[0] = decays.BF("mu+", "mu-");
133  result.BR_hjtautau[0] = decays.BF("tau+", "tau-");
134  result.BR_hjWW[0] = decays.BF("W+", "W-");
135  result.BR_hjZZ[0] = decays.BF("Z0", "Z0");
136  result.BR_hjZga[0] = decays.BF("gamma", "Z0");
137  result.BR_hjgaga[0] = decays.BF("gamma", "gamma");
138  result.BR_hjgg[0] = decays.BF("g", "g");
139 
140  // Add the invisibles
141  result.BR_hjinvisible[0] = 0.;
142  for (std::vector<std::pair<str,str>>::const_iterator it = couplings.invisibles.begin(); it != couplings.invisibles.end(); ++it)
143  {
144  result.BR_hjinvisible[0] += decays.BF(it->first, it->second);
145  }
146 
147  // Retrieve cross-section ratios from the HiggsCouplingsTable
148  set_CS(result, couplings, 1);
149 
150  // Zero all heavy neutral higgs masses, widths and effective couplings
151  for(int i = 1; i < 3; i++)
152  {
153  result.Mh[i] = 0.;
154  result.deltaMh[i] = 0.;
155  result.hGammaTot[i] = 0.;
156  result.CP[i] = 0.;
157  result.BR_hjss[i] = 0.;
158  result.BR_hjcc[i] = 0.;
159  result.BR_hjbb[i] = 0.;
160  result.BR_hjmumu[i] = 0.;
161  result.BR_hjtautau[i] = 0.;
162  result.BR_hjWW[i] = 0.;
163  result.BR_hjZZ[i] = 0.;
164  result.BR_hjZga[i] = 0.;
165  result.BR_hjgaga[i] = 0.;
166  result.BR_hjgg[i] = 0.;
167  result.BR_hjinvisible[i] = 0.;
168  result.CS_lep_hjZ_ratio[i] = 0.;
169  result.CS_lep_bbhj_ratio[i] = 0.;
170  result.CS_lep_tautauhj_ratio[i] = 0.;
171  result.CS_gg_hj_ratio[i] = 0.;
172  result.CS_bb_hj_ratio[i] = 0.;
173  result.CS_bg_hjb_ratio[i] = 0.;
174  result.CS_ud_hjWp_ratio[i] = 0.;
175  result.CS_cs_hjWp_ratio[i] = 0.;
176  result.CS_ud_hjWm_ratio[i] = 0.;
177  result.CS_cs_hjWm_ratio[i] = 0.;
178  result.CS_gg_hjZ_ratio[i] = 0.;
179  result.CS_dd_hjZ_ratio[i] = 0.;
180  result.CS_uu_hjZ_ratio[i] = 0.;
181  result.CS_ss_hjZ_ratio[i] = 0.;
182  result.CS_cc_hjZ_ratio[i] = 0.;
183  result.CS_bb_hjZ_ratio[i] = 0.;
184  result.CS_tev_vbf_ratio[i] = 0.;
185  result.CS_tev_tthj_ratio[i] = 0.;
186  result.CS_lhc7_vbf_ratio[i] = 0.;
187  result.CS_lhc7_tthj_ratio[i] = 0.;
188  result.CS_lhc8_vbf_ratio[i] = 0.;
189  result.CS_lhc8_tthj_ratio[i] = 0.;
190  for(int j = 0; j < 3; j++) result.BR_hjhihi[i][j] = 0.;
191  for(int j = 0; j < 3; j++) result.CS_lep_hjhi_ratio[i][j] = 0.;
192  }
193 
194  // Zero all H+ masses, widths and effective couplings
195  result.MHplus[0] = 0.;
196  result.deltaMHplus[0] = 0.;
197  result.HpGammaTot[0] = 0.;
198  result.BR_tWpb = 0.;
199  result.BR_tHpjb[0] = 0.;
200  result.BR_Hpjcs[0] = 0.;
201  result.BR_Hpjcb[0] = 0.;
202  result.BR_Hptaunu[0] = 0.;
203  result.CS_lep_HpjHmi_ratio[0] = 0.;
204  }
205 
207  void SMLikeHiggs_ModelParameters(hb_ModelParameters &result)
208  {
209  using namespace Pipes::SMLikeHiggs_ModelParameters;
210  dep_bucket<Spectrum>* spectrum_dependency = nullptr;
211  if (ModelInUse("ScalarSingletDM_Z2") or ModelInUse("ScalarSingletDM_Z2_running")) spectrum_dependency = &Dep::ScalarSingletDM_Z2_spectrum;
212  else if (ModelInUse("ScalarSingletDM_Z3") or ModelInUse("ScalarSingletDM_Z3_running")) spectrum_dependency = &Dep::ScalarSingletDM_Z3_spectrum;
213  else if (ModelInUse("StandardModel_Higgs") or ModelInUse("StandardModel_Higgs_running")) spectrum_dependency = &Dep::SM_spectrum;
214  else ColliderBit_error().raise(LOCAL_INFO, "No valid model for SMLikeHiggs_ModelParameters.");
215  const SubSpectrum& spec = (*spectrum_dependency)->get_HE();
216  set_SMLikeHiggs_ModelParameters(spec, *Dep::Higgs_Couplings, result);
217  }
218 
220  void MSSMLikeHiggs_ModelParameters(hb_ModelParameters &result)
221  {
222  using namespace Pipes::MSSMLikeHiggs_ModelParameters;
223 
224  // Set up neutral Higgses
225  int n_neutral_higgses = Dep::Higgs_Couplings->get_n_neutral_higgs();
226 
227  // Set the CP of the Higgs states.
228  for (int i = 0; i < n_neutral_higgses; i++) result.CP[i] = Dep::Higgs_Couplings->CP[i];
229 
230  // Retrieve higgs partial widths
231  const std::vector<const DecayTable::Entry*>& h0_widths = Dep::Higgs_Couplings->get_neutral_decays_array();
232  const DecayTable::Entry& H_plus_widths = Dep::Higgs_Couplings->get_charged_decays(0);
233  const DecayTable::Entry& t_widths = Dep::Higgs_Couplings->get_t_decays();
234 
235  // Pick the correct spectrum and specify the Higgses
236  dep_bucket<Spectrum>* spectrum_dependency;
237  std::vector<str> Higgses;
238 
239  if (ModelInUse("MSSM63atMGUT") or ModelInUse("MSSM63atQ"))
240  {
241  spectrum_dependency = &Dep::MSSM_spectrum;
242  Higgses = initVector<str>("h0_1", "h0_2", "A0");
243  }
244  else ColliderBit_error().raise(LOCAL_INFO, "No valid model for MSSMLikeHiggs_ModelParameters.");
245 
246 
247  const SubSpectrum& spec = (*spectrum_dependency)->get_HE();
248  static const std::vector<str> sHneut(Higgses);
249 
250  // Neutral higgs masses and errors
251  for(int i = 0; i < n_neutral_higgses; i++)
252  {
253  result.Mh[i] = spec.get(Par::Pole_Mass,sHneut[i]);
254  double upper = spec.get(Par::Pole_Mass_1srd_high,sHneut[i]);
255  double lower = spec.get(Par::Pole_Mass_1srd_low,sHneut[i]);
256  result.deltaMh[i] = result.Mh[i] * std::max(upper,lower);
257  }
258 
259  // Loop over all neutral Higgses, setting their branching fractions and total widths.
260  for(int i = 0; i < n_neutral_higgses; i++)
261  {
262  result.hGammaTot[i] = h0_widths[i]->width_in_GeV;
263  result.BR_hjss[i] = h0_widths[i]->BF("s", "sbar");
264  result.BR_hjcc[i] = h0_widths[i]->BF("c", "cbar");
265  result.BR_hjbb[i] = h0_widths[i]->BF("b", "bbar");
266  result.BR_hjmumu[i] = h0_widths[i]->BF("mu+", "mu-");
267  result.BR_hjtautau[i] = h0_widths[i]->BF("tau+", "tau-");
268  result.BR_hjWW[i] = h0_widths[i]->has_channel("W+", "W-") ? h0_widths[i]->BF("W+", "W-") : 0.0;
269  result.BR_hjZZ[i] = h0_widths[i]->has_channel("Z0", "Z0") ? h0_widths[i]->BF("Z0", "Z0") : 0.0;
270  result.BR_hjZga[i] = h0_widths[i]->has_channel("gamma", "Z0") ? h0_widths[i]->BF("gamma", "Z0") : 0.0;
271  result.BR_hjgaga[i] = h0_widths[i]->BF("gamma", "gamma");
272  result.BR_hjgg[i] = h0_widths[i]->BF("g", "g");
273  // Do decays to invisibles
274  result.BR_hjinvisible[i] = 0.;
275  for (std::vector<std::pair<str,str>>::const_iterator it = Dep::Higgs_Couplings->invisibles.begin(); it != Dep::Higgs_Couplings->invisibles.end(); ++it)
276  {
277  result.BR_hjinvisible[i] += h0_widths[i]->BF(it->first, it->second);
278  }
279  // Do decays to other neutral higgses
280  for (int j = 0; j < n_neutral_higgses; j++)
281  {
282  if (2.*result.Mh[j] < result.Mh[i] and h0_widths[i]->has_channel(sHneut[j],sHneut[j]))
283  {
284  result.BR_hjhihi[i][j] = h0_widths[i]->BF(sHneut[j],sHneut[j]);
285  }
286  else
287  {
288  result.BR_hjhihi[i][j] = 0.;
289  }
290  }
291  }
292 
293  // Charged higgs masses and errors
294  result.MHplus[0] = spec.get(Par::Pole_Mass,"H+");
295  double upper = spec.get(Par::Pole_Mass_1srd_high,"H+");
296  double lower = spec.get(Par::Pole_Mass_1srd_low,"H+");
297  result.deltaMHplus[0] = result.MHplus[0] * std::max(upper,lower);
298 
299  // Set charged Higgs branching fractions and total width.
300  result.HpGammaTot[0] = H_plus_widths.width_in_GeV;
301  result.BR_Hpjcs[0] = H_plus_widths.BF("c", "sbar");
302  result.BR_Hpjcb[0] = H_plus_widths.BF("c", "bbar");
303  result.BR_Hptaunu[0] = H_plus_widths.BF("tau+", "nu_tau");
304 
305  // Set top branching fractions
306  result.BR_tWpb = t_widths.BF("W+", "b");
307  result.BR_tHpjb[0] = t_widths.has_channel("H+", "b") ? t_widths.BF("H+", "b") : 0.0;
308 
309  // Retrieve cross-section ratios from the HiggsCouplingsTable
310  set_CS(result, *Dep::Higgs_Couplings, n_neutral_higgses);
311  }
312 
314  void calc_HB_LEP_LogLike(double &result)
315  {
316  using namespace Pipes::calc_HB_LEP_LogLike;
317 
318  hb_ModelParameters ModelParam = *Dep::HB_ModelParameters;
319 
320  Farray<double, 1,3, 1,3> CS_lep_hjhi_ratio;
321  Farray<double, 1,3, 1,3> BR_hjhihi;
322  // Transpose to get around Fortran matrix types
323  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
324  {
325  CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
326  BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
327  }
328  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
329  {
330  ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
331  ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
332  }
333 
334 
335  BEreq::HiggsBounds_neutral_input_part(&ModelParam.Mh[0], &ModelParam.hGammaTot[0], &ModelParam.CP[0],
336  &ModelParam.CS_lep_hjZ_ratio[0], &ModelParam.CS_lep_bbhj_ratio[0],
337  &ModelParam.CS_lep_tautauhj_ratio[0], &ModelParam.CS_lep_hjhi_ratio[0][0],
338  &ModelParam.CS_gg_hj_ratio[0], &ModelParam.CS_bb_hj_ratio[0],
339  &ModelParam.CS_bg_hjb_ratio[0], &ModelParam.CS_ud_hjWp_ratio[0],
340  &ModelParam.CS_cs_hjWp_ratio[0], &ModelParam.CS_ud_hjWm_ratio[0],
341  &ModelParam.CS_cs_hjWm_ratio[0], &ModelParam.CS_gg_hjZ_ratio[0],
342  &ModelParam.CS_dd_hjZ_ratio[0], &ModelParam.CS_uu_hjZ_ratio[0],
343  &ModelParam.CS_ss_hjZ_ratio[0], &ModelParam.CS_cc_hjZ_ratio[0],
344  &ModelParam.CS_bb_hjZ_ratio[0], &ModelParam.CS_tev_vbf_ratio[0],
345  &ModelParam.CS_tev_tthj_ratio[0], &ModelParam.CS_lhc7_vbf_ratio[0],
346  &ModelParam.CS_lhc7_tthj_ratio[0], &ModelParam.CS_lhc8_vbf_ratio[0],
347  &ModelParam.CS_lhc8_tthj_ratio[0], &ModelParam.BR_hjss[0],
348  &ModelParam.BR_hjcc[0], &ModelParam.BR_hjbb[0],
349  &ModelParam.BR_hjmumu[0], &ModelParam.BR_hjtautau[0],
350  &ModelParam.BR_hjWW[0], &ModelParam.BR_hjZZ[0],
351  &ModelParam.BR_hjZga[0], &ModelParam.BR_hjgaga[0],
352  &ModelParam.BR_hjgg[0], &ModelParam.BR_hjinvisible[0], &ModelParam.BR_hjhihi[0][0]);
353 
354  // Transpose it back
355  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
356  {
357  CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
358  BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
359  }
360  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
361  {
362  ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
363  ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
364  }
365 
366 
367  BEreq::HiggsBounds_charged_input(&ModelParam.MHplus[0], &ModelParam.HpGammaTot[0], &ModelParam.CS_lep_HpjHmi_ratio[0],
368  &ModelParam.BR_tWpb, &ModelParam.BR_tHpjb[0], &ModelParam.BR_Hpjcs[0],
369  &ModelParam.BR_Hpjcb[0], &ModelParam.BR_Hptaunu[0]);
370 
371  BEreq::HiggsBounds_set_mass_uncertainties(&ModelParam.deltaMh[0],&ModelParam.deltaMHplus[0]);
372 
373  // run Higgs bounds 'classic'
374  double obsratio;
375  int HBresult, chan, ncombined;
376  BEreq::run_HiggsBounds_classic(HBresult,chan,obsratio,ncombined);
377 
378  // extract the LEP chisq
379  double chisq_withouttheory,chisq_withtheory;
380  int chan2;
381  double theor_unc = 1.5; // theory uncertainty
382  BEreq::HB_calc_stats(theor_unc,chisq_withouttheory,chisq_withtheory,chan2);
383 
384  // Catch HiggsBound's error value, chisq = -999
385  if( fabs(chisq_withouttheory - (-999.)) < 1e-6)
386  {
387  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.");
388  chisq_withouttheory = 0.0;
389  }
390 
391  result = -0.5*chisq_withouttheory;
392  }
393 
395  void calc_HS_LHC_LogLike(double &result)
396  {
397  using namespace Pipes::calc_HS_LHC_LogLike;
398 
399  hb_ModelParameters ModelParam = *Dep::HB_ModelParameters;
400 
401  Farray<double, 1,3, 1,3> CS_lep_hjhi_ratio;
402  Farray<double, 1,3, 1,3> BR_hjhihi;
403  // Transpose to get around Fortran matrix types
404  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
405  {
406  CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
407  BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
408  }
409  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
410  {
411  ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
412  ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
413  }
414 
415 
416  BEreq::HiggsBounds_neutral_input_part_HS(&ModelParam.Mh[0], &ModelParam.hGammaTot[0], &ModelParam.CP[0],
417  &ModelParam.CS_lep_hjZ_ratio[0], &ModelParam.CS_lep_bbhj_ratio[0],
418  &ModelParam.CS_lep_tautauhj_ratio[0], &ModelParam.CS_lep_hjhi_ratio[0][0],
419  &ModelParam.CS_gg_hj_ratio[0], &ModelParam.CS_bb_hj_ratio[0],
420  &ModelParam.CS_bg_hjb_ratio[0], &ModelParam.CS_ud_hjWp_ratio[0],
421  &ModelParam.CS_cs_hjWp_ratio[0], &ModelParam.CS_ud_hjWm_ratio[0],
422  &ModelParam.CS_cs_hjWm_ratio[0], &ModelParam.CS_gg_hjZ_ratio[0],
423  &ModelParam.CS_dd_hjZ_ratio[0], &ModelParam.CS_uu_hjZ_ratio[0],
424  &ModelParam.CS_ss_hjZ_ratio[0], &ModelParam.CS_cc_hjZ_ratio[0],
425  &ModelParam.CS_bb_hjZ_ratio[0], &ModelParam.CS_tev_vbf_ratio[0],
426  &ModelParam.CS_tev_tthj_ratio[0], &ModelParam.CS_lhc7_vbf_ratio[0],
427  &ModelParam.CS_lhc7_tthj_ratio[0], &ModelParam.CS_lhc8_vbf_ratio[0],
428  &ModelParam.CS_lhc8_tthj_ratio[0], &ModelParam.BR_hjss[0],
429  &ModelParam.BR_hjcc[0], &ModelParam.BR_hjbb[0],
430  &ModelParam.BR_hjmumu[0], &ModelParam.BR_hjtautau[0],
431  &ModelParam.BR_hjWW[0], &ModelParam.BR_hjZZ[0],
432  &ModelParam.BR_hjZga[0], &ModelParam.BR_hjgaga[0],
433  &ModelParam.BR_hjgg[0], &ModelParam.BR_hjinvisible[0], &ModelParam.BR_hjhihi[0][0]);
434 
435  // Transpose it back
436  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
437  {
438  CS_lep_hjhi_ratio(i+1,j+1) = ModelParam.CS_lep_hjhi_ratio[i][j];
439  BR_hjhihi(i+1,j+1) = ModelParam.BR_hjhihi[i][j];
440  }
441  for(int i = 0; i < 3; i++) for(int j = 0; j < 3; j++)
442  {
443  ModelParam.CS_lep_hjhi_ratio[j][i] = CS_lep_hjhi_ratio(i+1,j+1);
444  ModelParam.BR_hjhihi[j][i] = BR_hjhihi(i+1,j+1);
445  }
446 
447 
448  BEreq::HiggsBounds_charged_input_HS(&ModelParam.MHplus[0], &ModelParam.HpGammaTot[0], &ModelParam.CS_lep_HpjHmi_ratio[0],
449  &ModelParam.BR_tWpb, &ModelParam.BR_tHpjb[0], &ModelParam.BR_Hpjcs[0],
450  &ModelParam.BR_Hpjcb[0], &ModelParam.BR_Hptaunu[0]);
451 
452  BEreq::HiggsSignals_neutral_input_MassUncertainty(&ModelParam.deltaMh[0]);
453 
454  // add uncertainties to cross-sections and branching ratios
455  // double dCS[5] = {0.,0.,0.,0.,0.};
456  // double dBR[5] = {0.,0.,0.,0.,0.};
457  // BEreq::setup_rate_uncertainties(dCS,dBR);
458 
459  // run HiggsSignals
460  int mode = 1; // 1- peak-centered chi2 method (recommended)
461  double csqmu, csqmh, csqtot, Pvalue;
462  int nobs;
463  BEreq::run_HiggsSignals(mode, csqmu, csqmh, csqtot, nobs, Pvalue);
464 
465  result = -0.5*csqtot;
466 
467  #ifdef COLLIDERBIT_DEBUG
468  std::ofstream f;
469  f.open ("HB_ModelParameters_contents.dat");
470  f<<"LHC log-likleihood";
471  for (int i = 0; i < 3; i++) f<<
472  " higgs index" <<
473  " "<<i<<":CP"<<
474  " "<<i<<":Mh"<<
475  " "<<i<<":hGammaTot"<<
476  " "<<i<<":CS_lep_hjZ_ratio"<<
477  " "<<i<<":CS_tev_vbf_ratio"<<
478  " "<<i<<":CS_lep_bbhj_ratio"<<
479  " "<<i<<":CS_lep_tautauhj_ratio"<<
480  " "<<i<<":CS_gg_hj_ratio"<<
481  " "<<i<<":CS_tev_tthj_ratio"<<
482  " "<<i<<":CS_lhc7_tthj_ratio"<<
483  " "<<i<<":CS_lhc8_tthj_ratio"<<
484  " "<<i<<":CS_lep_hjhi_ratio[0]"<<
485  " "<<i<<":CS_lep_hjhi_ratio[1]"<<
486  " "<<i<<":CS_lep_hjhi_ratio[2]"<<
487  " "<<i<<":BR_ss"<<
488  " "<<i<<":BR_cc"<<
489  " "<<i<<":BR_bb"<<
490  " "<<i<<":BR_mumu"<<
491  " "<<i<<":BR_tautau"<<
492  " "<<i<<":BR_WW"<<
493  " "<<i<<":BR_ZZ"<<
494  " "<<i<<":BR_Zga"<<
495  " "<<i<<":BR_gamgam"<<
496  " "<<i<<":BR_gg"<<
497  " "<<i<<":BR_invisible"<<
498  " "<<i<<":BR_hihi[0]"<<
499  " "<<i<<":BR_hihi[1]"<<
500  " "<<i<<":BR_hihi[2]";
501  f<<
502  " higgs index" <<
503  " "<<4<<"MHplus"<<
504  " "<<4<<":HpGammaTot"<<
505  " "<<4<<":CS_lep_HpjHmi_ratio"<<
506  " "<<4<<":BR_H+->cs"<<
507  " "<<4<<":BR_H+->cb"<<
508  " "<<4<<":BR_H+->taunu"<<
509  " "<<4<<":BR_t->W+b"<<
510  " "<<4<<":BR_t->H+b";
511  f << endl << std::setw(18) << result;
512  const int w = 24;
513  for (int i = 0; i < 3; i++)
514  {
515  f << std::setw(w) << i << std::setw(w) <<
516  ModelParam.CP[i] << std::setw(w) <<
517  ModelParam.Mh[i] << std::setw(w) <<
518  ModelParam.hGammaTot[i] << std::setw(w) <<
519  ModelParam.CS_lep_hjZ_ratio[i] << std::setw(w) <<
520  ModelParam.CS_tev_vbf_ratio[i] << std::setw(w) <<
521  ModelParam.CS_lep_bbhj_ratio[i] << std::setw(w) <<
522  ModelParam.CS_lep_tautauhj_ratio[i] << std::setw(w) <<
523  ModelParam.CS_gg_hj_ratio[i] << std::setw(w) <<
524  ModelParam.CS_tev_tthj_ratio[i] << std::setw(w) <<
525  ModelParam.CS_lhc7_tthj_ratio[i] << std::setw(w) <<
526  ModelParam.CS_lhc8_tthj_ratio[i];
527  for (int j = 0; j < 3; j++) f << std::setw(w) << ModelParam.CS_lep_hjhi_ratio[i][j];
528  f << std::setw(w) <<
529  ModelParam.BR_hjss[i] << std::setw(w) <<
530  ModelParam.BR_hjcc[i] << std::setw(w) <<
531  ModelParam.BR_hjbb[i] << std::setw(w) <<
532  ModelParam.BR_hjmumu[i] << std::setw(w) <<
533  ModelParam.BR_hjtautau[i] << std::setw(w) <<
534  ModelParam.BR_hjWW[i] << std::setw(w) <<
535  ModelParam.BR_hjZZ[i] << std::setw(w) <<
536  ModelParam.BR_hjZga[i] << std::setw(w) <<
537  ModelParam.BR_hjgaga[i] << std::setw(w) <<
538  ModelParam.BR_hjgg[i] << std::setw(w) <<
539  ModelParam.BR_hjinvisible[i];
540  for (int j = 0; j < 3; j++) f << std::setw(w) << ModelParam.BR_hjhihi[i][j];
541  }
542  f << std::setw(w) << 4 << std::setw(w) <<
543  ModelParam.MHplus[0] << std::setw(w) <<
544  ModelParam.HpGammaTot[0] << std::setw(w) <<
545  ModelParam.CS_lep_HpjHmi_ratio[0] << std::setw(w) <<
546  ModelParam.BR_Hpjcs[0] << std::setw(w) <<
547  ModelParam.BR_Hpjcb[0] << std::setw(w) <<
548  ModelParam.BR_Hptaunu[0] << std::setw(w) <<
549  ModelParam.BR_tWpb << std::setw(w) <<
550  ModelParam.BR_tHpjb[0];
551  f.close();
552  #endif
553 
554  }
555 
557  void FH_HiggsProd(fh_HiggsProd &result)
558  {
559  using namespace Pipes::FH_HiggsProd;
560 
561  Farray<fh_real, 1,52> prodxs;
562 
563  fh_HiggsProd HiggsProd;
564  int error;
565  fh_real sqrts;
566 
567  // Tevatron
568  sqrts = 2.;
569  error = 1;
570  BEreq::FHHiggsProd(error, sqrts, prodxs);
571  if (error != 0)
572  {
573  std::ostringstream err;
574  err << "BEreq::FHHiggsProd raised error flag for Tevatron: " << error << ".";
575  invalid_point().raise(err.str());
576  }
577  for(int i = 0; i < 52; i++) HiggsProd.prodxs_Tev[i] = prodxs(i+1);
578  // LHC7
579  sqrts = 7.;
580  error = 1;
581  BEreq::FHHiggsProd(error, sqrts, prodxs);
582  if (error != 0)
583  {
584  std::ostringstream err;
585  err << "BEreq::FHHiggsProd raised error flag for LHC7: " << error << ".";
586  invalid_point().raise(err.str());
587  }
588  for(int i = 0; i < 52; i++) HiggsProd.prodxs_LHC7[i] = prodxs(i+1);
589  // LHC8
590  sqrts = 8.;
591  error = 1;
592  BEreq::FHHiggsProd(error, sqrts, prodxs);
593  if (error != 0)
594  {
595  std::ostringstream err;
596  err << "BEreq::FHHiggsProd raised error flag for LHC8: " << error << ".";
597  invalid_point().raise(err.str());
598  }
599  for(int i = 0; i < 52; i++) HiggsProd.prodxs_LHC8[i] = prodxs(i+1);
600 
601  // The ttbar production cross-sections for the (BSM,SM) model can be found at (prodxs_X[h+27], prodxs_X[h+30]),
602  // where h is the higgs index (0 = h0_1, 1 = h0_2, 2 = A0) and X is one of Tev, LHC7 or LHC8.
603  result = HiggsProd;
604 
605  }
606 
607 
608  }
609 }
std::vector< std::vector< double > > C_hiZ2
Rollcall header for ColliderBit module.
Array class that matches the memory structure and functionality of arrays in Fortran codes Syntax: Fa...
Definition: util_types.hpp:318
DecayTable entry class. Holds the info on all decays of a given particle.
Definition: decay_table.hpp:96
double BF(const std::vector< std::pair< int, int > > &) const
Retrieve branching fraction for decay to a given final state.
void MSSMLikeHiggs_ModelParameters(hb_ModelParameters &result)
MSSM-like (MSSM + NMSSM + ...) Higgs model parameters for HiggsBounds/Signals.
const DecayTable::Entry & get_neutral_decays(int) const
Retrieve decays of a specific neutral Higgs, by index.
void set_CS(hb_ModelParameters &result, const HiggsCouplingsTable &couplings, int n_neutral_higgses)
Helper function to set HiggsBounds/Signals parameters cross-section ratios from a GAMBIT HiggsCouplin...
void FH_HiggsProd(fh_HiggsProd &result)
Higgs production cross-sections from FeynHiggs.
#define LOCAL_INFO
Definition: local_info.hpp:34
GAMBIT native higgs coupling table class.
General small utility classes, typedefs, etc.
void SMLikeHiggs_ModelParameters(hb_ModelParameters &result)
SM-like (SM + possible invisibles) Higgs model parameters for HiggsBounds/Signals.
std::vector< std::pair< str, str > > invisibles
Particles that higgses can decay invisibly to.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum ScalarSingletDM_Z3_spectrum
GAMBIT error class.
Definition: exceptions.hpp:136
void calc_HB_LEP_LogLike(double &result)
Get a LEP chisq from HiggsBounds.
virtual void raise(const std::string &)
Raise the exception, i.e. throw it. Exact override of base method.
Definition: exceptions.cpp:422
virtual double get(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
double width_in_GeV
Total particle width (in GeV)
bool has_channel(const std::vector< std::pair< int, int > > &) const
Check if a given final state exists in this DecayTable::Entry.
An interface class for module dependencies.
std::vector< double > C_WW2
Effective couplings for neutral higgses.
Header file that includes all GAMBIT headers required for a module source file.
Virtual base class for interacting with spectrum generator output.
Definition: subspectrum.hpp:87
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum SM_spectrum
invalid_point_exception & invalid_point()
Invalid point exceptions.
std::vector< double > CP
CP of neutral higgses.
void set_SMLikeHiggs_ModelParameters(const SubSpectrum &spec, const HiggsCouplingsTable &couplings, hb_ModelParameters &result)
Helper function for populating a HiggsBounds/Signals ModelParameters object for SM-like Higgs...
TODO: see if we can use this one:
Definition: Analysis.hpp:33
virtual bool has(const Par::Tags, const str &, const SpecOverrideOptions=use_overrides, const SafeBool check_antiparticle=SafeBool(true)) const =0
Getters/Setters etc.
void calc_HS_LHC_LogLike(double &result)
Get an LHC chisq from HiggsSignals.