gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
FlavBit.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
42 
43 #include <string>
44 #include <iostream>
45 #include <fstream>
46 #include <map>
47 
57 #include "gambit/cmake/cmake_variables.hpp"
58 
59 //#define FLAVBIT_DEBUG
60 //#define FLAVBIT_DEBUG_LL
61 
62 namespace Gambit
63 {
64 
65  namespace FlavBit
66  {
67 
68  using namespace std;
69  namespace ublas = boost::numeric::ublas;
70 
71  const bool flav_debug =
72  #ifdef FLAVBIT_DEBUG
73  true;
74  #else
75  false;
76  #endif
77 
78  const bool flav_debug_LL =
79  #ifdef FLAVBIT_DEBUG_LL
80  true;
81  #else
82  false;
83  #endif
84 
86  void SI_fill(parameters &result)
87  {
88  using namespace Pipes::SI_fill;
89  using namespace std;
90 
91  SLHAstruct spectrum;
92  // Obtain SLHAea object from spectrum
93  if (ModelInUse("WC"))
94  {
95  spectrum = Dep::SM_spectrum->getSLHAea(2);
96  }
97  else if (ModelInUse("MSSM63atMGUT") or ModelInUse("MSSM63atQ"))
98  {
99  spectrum = Dep::MSSM_spectrum->getSLHAea(2);
100  // Add the MODSEL block if it is not provided by the spectrum object.
101  SLHAea_add(spectrum,"MODSEL",1, 0, "General MSSM", false);
102  }
103  else
104  {
105  FlavBit_error().raise(LOCAL_INFO, "Unrecognised model.");
106  }
107 
108  BEreq::Init_param(&result);
109 
110  int ie,je;
111 
112  result.model=-1;
113  if (!spectrum["MODSEL"].empty())
114  {
115  if (spectrum["MODSEL"][1].is_data_line()) result.model=SLHAea::to<int>(spectrum["MODSEL"][1][1]);
116  if (spectrum["MODSEL"][3].is_data_line()) result.NMSSM=SLHAea::to<int>(spectrum["MODSEL"][3][1]);
117  if (spectrum["MODSEL"][4].is_data_line()) result.RV=SLHAea::to<int>(spectrum["MODSEL"][4][1]);
118  if (spectrum["MODSEL"][5].is_data_line()) result.CPV=SLHAea::to<int>(spectrum["MODSEL"][5][1]);
119  if (spectrum["MODSEL"][6].is_data_line()) result.FV=SLHAea::to<int>(spectrum["MODSEL"][6][1]);
120  if (spectrum["MODSEL"][12].is_data_line()) result.Q=SLHAea::to<double>(spectrum["MODSEL"][12][1]);
121  }
122 
123  if (result.NMSSM != 0) result.model=result.NMSSM;
124  if (result.RV != 0) result.model=-2;
125  if (result.CPV != 0) result.model=-2;
126 
127  if (!spectrum["SMINPUTS"].empty())
128  {
129  if (spectrum["SMINPUTS"][1].is_data_line()) result.inv_alpha_em=SLHAea::to<double>(spectrum["SMINPUTS"][1][1]);
130  if (spectrum["SMINPUTS"][2].is_data_line()) result.Gfermi=SLHAea::to<double>(spectrum["SMINPUTS"][2][1]);
131  if (spectrum["SMINPUTS"][3].is_data_line()) result.alphas_MZ=SLHAea::to<double>(spectrum["SMINPUTS"][3][1]);
132  if (spectrum["SMINPUTS"][4].is_data_line()) result.mass_Z=SLHAea::to<double>(spectrum["SMINPUTS"][4][1]);
133  if (spectrum["SMINPUTS"][5].is_data_line()) result.mass_b=SLHAea::to<double>(spectrum["SMINPUTS"][5][1]);
134  if (spectrum["SMINPUTS"][6].is_data_line()) result.mass_top_pole=SLHAea::to<double>(spectrum["SMINPUTS"][6][1]);
135  if (spectrum["SMINPUTS"][7].is_data_line()) result.mass_tau=SLHAea::to<double>(spectrum["SMINPUTS"][7][1]);
136  if (spectrum["SMINPUTS"][8].is_data_line()) result.mass_nutau2=SLHAea::to<double>(spectrum["SMINPUTS"][8][1]);
137  if (spectrum["SMINPUTS"][11].is_data_line()) result.mass_e2=SLHAea::to<double>(spectrum["SMINPUTS"][11][1]);
138  if (spectrum["SMINPUTS"][12].is_data_line()) result.mass_nue2=SLHAea::to<double>(spectrum["SMINPUTS"][12][1]);
139  if (spectrum["SMINPUTS"][13].is_data_line()) result.mass_mu2=SLHAea::to<double>(spectrum["SMINPUTS"][13][1]);
140  if (spectrum["SMINPUTS"][14].is_data_line()) result.mass_numu2=SLHAea::to<double>(spectrum["SMINPUTS"][14][1]);
141  if (spectrum["SMINPUTS"][21].is_data_line()) result.mass_d2=SLHAea::to<double>(spectrum["SMINPUTS"][21][1]);
142  if (spectrum["SMINPUTS"][22].is_data_line()) result.mass_u2=SLHAea::to<double>(spectrum["SMINPUTS"][22][1]);
143  if (spectrum["SMINPUTS"][23].is_data_line()) result.mass_s2=SLHAea::to<double>(spectrum["SMINPUTS"][23][1]);
144  if (spectrum["SMINPUTS"][24].is_data_line()) result.mass_c2=SLHAea::to<double>(spectrum["SMINPUTS"][24][1]);
145  }
146 
147  if (!spectrum["VCKMIN"].empty())
148  {
149  if (spectrum["VCKMIN"][1].is_data_line()) result.CKM_lambda=SLHAea::to<double>(spectrum["VCKMIN"][1][1]);
150  if (spectrum["VCKMIN"][2].is_data_line()) result.CKM_A=SLHAea::to<double>(spectrum["VCKMIN"][2][1]);
151  if (spectrum["VCKMIN"][3].is_data_line()) result.CKM_rhobar=SLHAea::to<double>(spectrum["VCKMIN"][3][1]);
152  if (spectrum["VCKMIN"][4].is_data_line()) result.CKM_etabar=SLHAea::to<double>(spectrum["VCKMIN"][4][1]);
153  }
154 
155  if (!spectrum["UPMNSIN"].empty())
156  {
157  if (spectrum["UPMNSIN"][1].is_data_line()) result.PMNS_theta12=SLHAea::to<double>(spectrum["UPMNSIN"][1][1]);
158  if (spectrum["UPMNSIN"][2].is_data_line()) result.PMNS_theta23=SLHAea::to<double>(spectrum["UPMNSIN"][2][1]);
159  if (spectrum["UPMNSIN"][3].is_data_line()) result.PMNS_theta13=SLHAea::to<double>(spectrum["UPMNSIN"][3][1]);
160  if (spectrum["UPMNSIN"][4].is_data_line()) result.PMNS_delta13=SLHAea::to<double>(spectrum["UPMNSIN"][4][1]);
161  if (spectrum["UPMNSIN"][5].is_data_line()) result.PMNS_alpha1=SLHAea::to<double>(spectrum["UPMNSIN"][5][1]);
162  if (spectrum["UPMNSIN"][6].is_data_line()) result.PMNS_alpha2=SLHAea::to<double>(spectrum["UPMNSIN"][6][1]);
163  }
164 
165  if (!spectrum["MINPAR"].empty())
166  {
167  if (spectrum["MINPAR"][3].is_data_line()) result.tan_beta=SLHAea::to<double>(spectrum["MINPAR"][3][1]);
168  switch(result.model)
169  {
170  case 1:
171  {
172  if (spectrum["MINPAR"][1].is_data_line()) result.m0=SLHAea::to<double>(spectrum["MINPAR"][1][1]);
173  if (spectrum["MINPAR"][2].is_data_line()) result.m12=SLHAea::to<double>(spectrum["MINPAR"][2][1]);
174  if (spectrum["MINPAR"][4].is_data_line()) result.sign_mu=SLHAea::to<double>(spectrum["MINPAR"][4][1]);
175  if (spectrum["MINPAR"][5].is_data_line()) result.A0=SLHAea::to<double>(spectrum["MINPAR"][5][1]);
176  break;
177  }
178  case 2:
179  {
180  if (spectrum["MINPAR"][1].is_data_line()) result.Lambda=SLHAea::to<double>(spectrum["MINPAR"][1][1]);
181  if (spectrum["MINPAR"][2].is_data_line()) result.Mmess=SLHAea::to<double>(spectrum["MINPAR"][2][1]);
182  if (spectrum["MINPAR"][4].is_data_line()) result.sign_mu=SLHAea::to<double>(spectrum["MINPAR"][4][1]);
183  if (spectrum["MINPAR"][5].is_data_line()) result.N5=SLHAea::to<double>(spectrum["MINPAR"][5][1]);
184  if (spectrum["MINPAR"][6].is_data_line()) result.cgrav=SLHAea::to<double>(spectrum["MINPAR"][6][1]);
185  break;
186  }
187  case 3:
188  {
189  if (spectrum["MINPAR"][1].is_data_line()) result.m32=SLHAea::to<double>(spectrum["MINPAR"][1][1]);
190  if (spectrum["MINPAR"][2].is_data_line()) result.m0=SLHAea::to<double>(spectrum["MINPAR"][2][1]);
191  if (spectrum["MINPAR"][4].is_data_line()) result.sign_mu=SLHAea::to<double>(spectrum["MINPAR"][4][1]);
192  break;
193  }
194  }
195  }
196 
197  if (!spectrum["EXTPAR"].empty())
198  {
199  if (spectrum["EXTPAR"][0].is_data_line()) result.Min=SLHAea::to<double>(spectrum["EXTPAR"][0][1]);
200  if (spectrum["EXTPAR"][1].is_data_line()) result.M1_Min=SLHAea::to<double>(spectrum["EXTPAR"][1][1]);
201  if (spectrum["EXTPAR"][2].is_data_line()) result.M2_Min=SLHAea::to<double>(spectrum["EXTPAR"][2][1]);
202  if (spectrum["EXTPAR"][3].is_data_line()) result.M3_Min=SLHAea::to<double>(spectrum["EXTPAR"][3][1]);
203  if (spectrum["EXTPAR"][11].is_data_line()) result.At_Min=SLHAea::to<double>(spectrum["EXTPAR"][11][1]);
204  if (spectrum["EXTPAR"][12].is_data_line()) result.Ab_Min=SLHAea::to<double>(spectrum["EXTPAR"][12][1]);
205  if (spectrum["EXTPAR"][13].is_data_line()) result.Atau_Min=SLHAea::to<double>(spectrum["EXTPAR"][13][1]);
206  if (spectrum["EXTPAR"][21].is_data_line()) result.M2H1_Min=SLHAea::to<double>(spectrum["EXTPAR"][21][1]);
207  if (spectrum["EXTPAR"][22].is_data_line()) result.M2H2_Min=SLHAea::to<double>(spectrum["EXTPAR"][22][1]);
208  if (spectrum["EXTPAR"][23].is_data_line()) result.mu_Min=SLHAea::to<double>(spectrum["EXTPAR"][23][1]);
209  if (spectrum["EXTPAR"][24].is_data_line()) result.M2A_Min=SLHAea::to<double>(spectrum["EXTPAR"][24][1]);
210  if (spectrum["EXTPAR"][25].is_data_line()) result.tb_Min=SLHAea::to<double>(spectrum["EXTPAR"][25][1]);
211  if (spectrum["EXTPAR"][26].is_data_line()) result.mA_Min=SLHAea::to<double>(spectrum["EXTPAR"][26][1]);
212  if (spectrum["EXTPAR"][31].is_data_line()) result.MeL_Min=SLHAea::to<double>(spectrum["EXTPAR"][31][1]);
213  if (spectrum["EXTPAR"][32].is_data_line()) result.MmuL_Min=SLHAea::to<double>(spectrum["EXTPAR"][32][1]);
214  if (spectrum["EXTPAR"][33].is_data_line()) result.MtauL_Min=SLHAea::to<double>(spectrum["EXTPAR"][33][1]);
215  if (spectrum["EXTPAR"][34].is_data_line()) result.MeR_Min=SLHAea::to<double>(spectrum["EXTPAR"][34][1]);
216  if (spectrum["EXTPAR"][35].is_data_line()) result.MmuR_Min=SLHAea::to<double>(spectrum["EXTPAR"][35][1]);
217  if (spectrum["EXTPAR"][36].is_data_line()) result.MtauR_Min=SLHAea::to<double>(spectrum["EXTPAR"][36][1]);
218  if (spectrum["EXTPAR"][41].is_data_line()) result.MqL1_Min=SLHAea::to<double>(spectrum["EXTPAR"][41][1]);
219  if (spectrum["EXTPAR"][42].is_data_line()) result.MqL2_Min=SLHAea::to<double>(spectrum["EXTPAR"][42][1]);
220  if (spectrum["EXTPAR"][43].is_data_line()) result.MqL3_Min=SLHAea::to<double>(spectrum["EXTPAR"][43][1]);
221  if (spectrum["EXTPAR"][44].is_data_line()) result.MuR_Min=SLHAea::to<double>(spectrum["EXTPAR"][44][1]);
222  if (spectrum["EXTPAR"][45].is_data_line()) result.McR_Min=SLHAea::to<double>(spectrum["EXTPAR"][45][1]);
223  if (spectrum["EXTPAR"][46].is_data_line()) result.MtR_Min=SLHAea::to<double>(spectrum["EXTPAR"][46][1]);
224  if (spectrum["EXTPAR"][47].is_data_line()) result.MdR_Min=SLHAea::to<double>(spectrum["EXTPAR"][47][1]);
225  if (spectrum["EXTPAR"][48].is_data_line()) result.MsR_Min=SLHAea::to<double>(spectrum["EXTPAR"][48][1]);
226  if (spectrum["EXTPAR"][49].is_data_line()) result.MbR_Min=SLHAea::to<double>(spectrum["EXTPAR"][49][1]);
227  if (spectrum["EXTPAR"][51].is_data_line()) result.N51=SLHAea::to<double>(spectrum["EXTPAR"][51][1]);
228  if (spectrum["EXTPAR"][52].is_data_line()) result.N52=SLHAea::to<double>(spectrum["EXTPAR"][52][1]);
229  if (spectrum["EXTPAR"][53].is_data_line()) result.N53=SLHAea::to<double>(spectrum["EXTPAR"][53][1]);
230  if (spectrum["EXTPAR"][61].is_data_line()) result.lambdaNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][61][1]);
231  if (spectrum["EXTPAR"][62].is_data_line()) result.kappaNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][62][1]);
232  if (spectrum["EXTPAR"][63].is_data_line()) result.AlambdaNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][63][1]);
233  if (spectrum["EXTPAR"][64].is_data_line()) result.AkappaNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][64][1]);
234  if (spectrum["EXTPAR"][65].is_data_line()) result.lambdaSNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][65][1]);
235  if (spectrum["EXTPAR"][66].is_data_line()) result.xiFNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][66][1]);
236  if (spectrum["EXTPAR"][67].is_data_line()) result.xiSNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][67][1]);
237  if (spectrum["EXTPAR"][68].is_data_line()) result.mupNMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][68][1]);
238  if (spectrum["EXTPAR"][69].is_data_line()) result.mSp2NMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][69][1]);
239  if (spectrum["EXTPAR"][70].is_data_line()) result.mS2NMSSM_Min=SLHAea::to<double>(spectrum["EXTPAR"][70][1]);
240  }
241 
242  if (!spectrum["MASS"].empty())
243  {
244  if (spectrum["MASS"][1].is_data_line()) result.mass_d=SLHAea::to<double>(spectrum["MASS"][1][1]);
245  if (spectrum["MASS"][2].is_data_line()) result.mass_u=SLHAea::to<double>(spectrum["MASS"][2][1]);
246  if (spectrum["MASS"][3].is_data_line()) result.mass_s=SLHAea::to<double>(spectrum["MASS"][3][1]);
247  if (spectrum["MASS"][4].is_data_line()) result.mass_c=SLHAea::to<double>(spectrum["MASS"][4][1]);
248  if (spectrum["MASS"][6].is_data_line()) result.mass_t=SLHAea::to<double>(spectrum["MASS"][6][1]);
249  if (spectrum["MASS"][11].is_data_line()) result.mass_e=SLHAea::to<double>(spectrum["MASS"][11][1]);
250  if (spectrum["MASS"][12].is_data_line()) result.mass_nue=SLHAea::to<double>(spectrum["MASS"][12][1]);
251  if (spectrum["MASS"][13].is_data_line()) result.mass_mu=SLHAea::to<double>(spectrum["MASS"][13][1]);
252  if (spectrum["MASS"][14].is_data_line()) result.mass_num=SLHAea::to<double>(spectrum["MASS"][14][1]);
253  if (spectrum["MASS"][15].is_data_line()) result.mass_tau=result.mass_tau_pole=SLHAea::to<double>(spectrum["MASS"][15][1]);
254  if (spectrum["MASS"][16].is_data_line()) result.mass_nut=SLHAea::to<double>(spectrum["MASS"][16][1]);
255  if (spectrum["MASS"][21].is_data_line()) result.mass_gluon=SLHAea::to<double>(spectrum["MASS"][21][1]);
256  if (spectrum["MASS"][22].is_data_line()) result.mass_photon=SLHAea::to<double>(spectrum["MASS"][22][1]);
257  if (spectrum["MASS"][23].is_data_line()) result.mass_Z0=SLHAea::to<double>(spectrum["MASS"][23][1]);
258  if (spectrum["MASS"][24].is_data_line()) result.mass_W=SLHAea::to<double>(spectrum["MASS"][24][1]);
259  if (spectrum["MASS"][25].is_data_line()) result.mass_h0=SLHAea::to<double>(spectrum["MASS"][25][1]);
260  if (spectrum["MASS"][35].is_data_line()) result.mass_H0=SLHAea::to<double>(spectrum["MASS"][35][1]);
261  if (spectrum["MASS"][36].is_data_line()) result.mass_A0=SLHAea::to<double>(spectrum["MASS"][36][1]);
262  if (spectrum["MASS"][37].is_data_line()) result.mass_H=SLHAea::to<double>(spectrum["MASS"][37][1]);
263  if (spectrum["MASS"][39].is_data_line()) result.mass_graviton=SLHAea::to<double>(spectrum["MASS"][39][1]);
264  if (spectrum["MASS"][45].is_data_line()) result.mass_H03=SLHAea::to<double>(spectrum["MASS"][45][1]);
265  if (spectrum["MASS"][46].is_data_line()) result.mass_A02=SLHAea::to<double>(spectrum["MASS"][46][1]);
266  if (spectrum["MASS"][1000001].is_data_line()) result.mass_dnl=SLHAea::to<double>(spectrum["MASS"][1000001][1]);
267  if (spectrum["MASS"][1000002].is_data_line()) result.mass_upl=SLHAea::to<double>(spectrum["MASS"][1000002][1]);
268  if (spectrum["MASS"][1000003].is_data_line()) result.mass_stl=SLHAea::to<double>(spectrum["MASS"][1000003][1]);
269  if (spectrum["MASS"][1000004].is_data_line()) result.mass_chl=SLHAea::to<double>(spectrum["MASS"][1000004][1]);
270  if (spectrum["MASS"][1000005].is_data_line()) result.mass_b1=SLHAea::to<double>(spectrum["MASS"][1000005][1]);
271  if (spectrum["MASS"][1000006].is_data_line()) result.mass_t1=SLHAea::to<double>(spectrum["MASS"][1000006][1]);
272  if (spectrum["MASS"][1000011].is_data_line()) result.mass_el=SLHAea::to<double>(spectrum["MASS"][1000011][1]);
273  if (spectrum["MASS"][1000012].is_data_line()) result.mass_nuel=SLHAea::to<double>(spectrum["MASS"][1000012][1]);
274  if (spectrum["MASS"][1000013].is_data_line()) result.mass_mul=SLHAea::to<double>(spectrum["MASS"][1000013][1]);
275  if (spectrum["MASS"][1000014].is_data_line()) result.mass_numl=SLHAea::to<double>(spectrum["MASS"][1000014][1]);
276  if (spectrum["MASS"][1000015].is_data_line()) result.mass_tau1=SLHAea::to<double>(spectrum["MASS"][1000015][1]);
277  if (spectrum["MASS"][1000016].is_data_line()) result.mass_nutl=SLHAea::to<double>(spectrum["MASS"][1000016][1]);
278  if (spectrum["MASS"][1000021].is_data_line()) result.mass_gluino=SLHAea::to<double>(spectrum["MASS"][1000021][1]);
279  if (spectrum["MASS"][1000022].is_data_line()) result.mass_neut[1]=SLHAea::to<double>(spectrum["MASS"][1000022][1]);
280  if (spectrum["MASS"][1000023].is_data_line()) result.mass_neut[2]=SLHAea::to<double>(spectrum["MASS"][1000023][1]);
281  if (spectrum["MASS"][1000024].is_data_line()) result.mass_cha1=SLHAea::to<double>(spectrum["MASS"][1000024][1]);
282  if (spectrum["MASS"][1000025].is_data_line()) result.mass_neut[3]=SLHAea::to<double>(spectrum["MASS"][1000025][1]);
283  if (spectrum["MASS"][1000035].is_data_line()) result.mass_neut[4]=SLHAea::to<double>(spectrum["MASS"][1000035][1]);
284  if (spectrum["MASS"][1000037].is_data_line()) result.mass_cha2=SLHAea::to<double>(spectrum["MASS"][1000037][1]);
285  if (spectrum["MASS"][1000039].is_data_line()) result.mass_gravitino=SLHAea::to<double>(spectrum["MASS"][1000039][1]);
286  if (spectrum["MASS"][1000045].is_data_line()) result.mass_neut[5]=SLHAea::to<double>(spectrum["MASS"][1000045][1]);
287  if (spectrum["MASS"][2000001].is_data_line()) result.mass_dnr=SLHAea::to<double>(spectrum["MASS"][2000001][1]);
288  if (spectrum["MASS"][2000002].is_data_line()) result.mass_upr=SLHAea::to<double>(spectrum["MASS"][2000002][1]);
289  if (spectrum["MASS"][2000003].is_data_line()) result.mass_str=SLHAea::to<double>(spectrum["MASS"][2000003][1]);
290  if (spectrum["MASS"][2000004].is_data_line()) result.mass_chr=SLHAea::to<double>(spectrum["MASS"][2000004][1]);
291  if (spectrum["MASS"][2000005].is_data_line()) result.mass_b2=SLHAea::to<double>(spectrum["MASS"][2000005][1]);
292  if (spectrum["MASS"][2000006].is_data_line()) result.mass_t2=SLHAea::to<double>(spectrum["MASS"][2000006][1]);
293  if (spectrum["MASS"][2000011].is_data_line()) result.mass_er=SLHAea::to<double>(spectrum["MASS"][2000011][1]);
294  if (spectrum["MASS"][2000012].is_data_line()) result.mass_nuer=SLHAea::to<double>(spectrum["MASS"][2000012][1]);
295  if (spectrum["MASS"][2000013].is_data_line()) result.mass_mur=SLHAea::to<double>(spectrum["MASS"][2000013][1]);
296  if (spectrum["MASS"][2000014].is_data_line()) result.mass_numr=SLHAea::to<double>(spectrum["MASS"][2000014][1]);
297  if (spectrum["MASS"][2000015].is_data_line()) result.mass_tau2=SLHAea::to<double>(spectrum["MASS"][2000015][1]);
298  if (spectrum["MASS"][2000016].is_data_line()) result.mass_nutr=SLHAea::to<double>(spectrum["MASS"][2000016][1]);
299  }
300 
301  // The following blocks will only appear for SUSY models so let's not waste time checking them if we're not scanning one of those.
302  if (ModelInUse("MSSM63atMGUT") or ModelInUse("MSSM63atQ"))
303  {
304  // The scale doesn't come through in MODSEL with all spectrum generators
305  result.Q = Dep::MSSM_spectrum->get_HE().GetScale();
306 
307  if (!spectrum["ALPHA"].empty()) if (spectrum["ALPHA"].back().is_data_line()) result.alpha=SLHAea::to<double>(spectrum["ALPHA"].back().at(0));
308 
309  if (!spectrum["STOPMIX"].empty()) for (ie=1;ie<=2;ie++) for (je=1;je<=2;je++)
310  if (spectrum["STOPMIX"][max(ie,je)].is_data_line()) result.stop_mix[ie][je]=SLHAea::to<double>(spectrum["STOPMIX"].at(ie,je)[2]);
311  if (!spectrum["SBOTMIX"].empty()) for (ie=1;ie<=2;ie++) for (je=1;je<=2;je++)
312  if (spectrum["SBOTMIX"][max(ie,je)].is_data_line()) result.sbot_mix[ie][je]=SLHAea::to<double>(spectrum["SBOTMIX"].at(ie,je)[2]);
313  if (!spectrum["STAUMIX"].empty()) for (ie=1;ie<=2;ie++) for (je=1;je<=2;je++)
314  if (spectrum["STAUMIX"][max(ie,je)].is_data_line()) result.stau_mix[ie][je]=SLHAea::to<double>(spectrum["STAUMIX"].at(ie,je)[2]);
315  if (!spectrum["NMIX"].empty()) for (ie=1;ie<=4;ie++) for (je=1;je<=4;je++)
316  if (spectrum["NMIX"][max(ie,je)].is_data_line()) result.neut_mix[ie][je]=SLHAea::to<double>(spectrum["NMIX"].at(ie,je)[2]);
317  if (!spectrum["NMNMIX"].empty()) for (ie=1;ie<=5;ie++) for (je=1;je<=5;je++)
318  if (spectrum["NMNMIX"][max(ie,je)].is_data_line()) result.neut_mix[ie][je]=SLHAea::to<double>(spectrum["NMNMIX"].at(ie,je)[2]);
319  if (!spectrum["UMIX"].empty()) for (ie=1;ie<=2;ie++) for (je=1;je<=2;je++)
320  if (spectrum["UMIX"][max(ie,je)].is_data_line()) result.charg_Umix[ie][je]=SLHAea::to<double>(spectrum["UMIX"].at(ie,je)[2]);
321  if (!spectrum["VMIX"].empty()) for (ie=1;ie<=2;ie++) for (je=1;je<=2;je++)
322  if (spectrum["VMIX"][max(ie,je)].is_data_line()) result.charg_Vmix[ie][je]=SLHAea::to<double>(spectrum["VMIX"].at(ie,je)[2]);
323 
324  if (!spectrum["GAUGE"].empty())
325  {
326  if (spectrum["GAUGE"][1].is_data_line()) result.gp_Q=SLHAea::to<double>(spectrum["GAUGE"][1][1]);
327  if (spectrum["GAUGE"][2].is_data_line()) result.g2_Q=SLHAea::to<double>(spectrum["GAUGE"][2][1]);
328  if (spectrum["GAUGE"][3].is_data_line()) result.g3_Q=SLHAea::to<double>(spectrum["GAUGE"][3][1]);
329  }
330 
331  if (!spectrum["YU"].empty()) for (ie=1;ie<=3;ie++) if (spectrum["YU"][ie].is_data_line()) result.yut[ie]=SLHAea::to<double>(spectrum["YU"].at(ie,ie)[2]);
332  if (!spectrum["YD"].empty()) for (ie=1;ie<=3;ie++) if (spectrum["YD"][ie].is_data_line()) result.yub[ie]=SLHAea::to<double>(spectrum["YD"].at(ie,ie)[2]);
333  if (!spectrum["YE"].empty()) for (ie=1;ie<=3;ie++) if (spectrum["YE"][ie].is_data_line()) result.yutau[ie]=SLHAea::to<double>(spectrum["YE"].at(ie,ie)[2]);
334 
335  if (!spectrum["HMIX"].empty())
336  {
337  if (spectrum["HMIX"][1].is_data_line()) result.mu_Q=SLHAea::to<double>(spectrum["HMIX"][1][1]);
338  if (spectrum["HMIX"][2].is_data_line()) result.tanb_GUT=SLHAea::to<double>(spectrum["HMIX"][2][1]);
339  if (spectrum["HMIX"][3].is_data_line()) result.Higgs_VEV=SLHAea::to<double>(spectrum["HMIX"][3][1]);
340  if (spectrum["HMIX"][4].is_data_line()) result.mA2_Q=SLHAea::to<double>(spectrum["HMIX"][4][1]);
341  }
342 
343  if (!spectrum["NMHMIX"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
344  if (spectrum["NMHMIX"][max(ie,je)].is_data_line()) result.H0_mix[ie][je]=SLHAea::to<double>(spectrum["NMHMIX"].at(ie,je)[2]);
345 
346  if (!spectrum["NMAMIX"].empty()) for (ie=1;ie<=2;ie++) for (je=1;je<=2;je++)
347  if (spectrum["NMAMIX"][max(ie,je)].is_data_line()) result.A0_mix[ie][je]=SLHAea::to<double>(spectrum["NMAMIX"].at(ie,je)[2]);
348 
349  if (!spectrum["MSOFT"].empty())
350  {
351  if (!spectrum["MSOFT"].front().empty()) result.MSOFT_Q=SLHAea::to<double>(spectrum["MSOFT"].front().at(3));
352  if (spectrum["MSOFT"][1].is_data_line()) result.M1_Q=SLHAea::to<double>(spectrum["MSOFT"][1][1]);
353  if (spectrum["MSOFT"][2].is_data_line()) result.M2_Q=SLHAea::to<double>(spectrum["MSOFT"][2][1]);
354  if (spectrum["MSOFT"][3].is_data_line()) result.M3_Q=SLHAea::to<double>(spectrum["MSOFT"][3][1]);
355  if (spectrum["MSOFT"][21].is_data_line()) result.M2H1_Q=SLHAea::to<double>(spectrum["MSOFT"][21][1]);
356  if (spectrum["MSOFT"][22].is_data_line()) result.M2H2_Q=SLHAea::to<double>(spectrum["MSOFT"][22][1]);
357  if (spectrum["MSOFT"][31].is_data_line()) result.MeL_Q=SLHAea::to<double>(spectrum["MSOFT"][31][1]);
358  if (spectrum["MSOFT"][32].is_data_line()) result.MmuL_Q=SLHAea::to<double>(spectrum["MSOFT"][32][1]);
359  if (spectrum["MSOFT"][33].is_data_line()) result.MtauL_Q=SLHAea::to<double>(spectrum["MSOFT"][33][1]);
360  if (spectrum["MSOFT"][34].is_data_line()) result.MeR_Q=SLHAea::to<double>(spectrum["MSOFT"][34][1]);
361  if (spectrum["MSOFT"][35].is_data_line()) result.MmuR_Q=SLHAea::to<double>(spectrum["MSOFT"][35][1]);
362  if (spectrum["MSOFT"][36].is_data_line()) result.MtauR_Q=SLHAea::to<double>(spectrum["MSOFT"][36][1]);
363  if (spectrum["MSOFT"][41].is_data_line()) result.MqL1_Q=SLHAea::to<double>(spectrum["MSOFT"][41][1]);
364  if (spectrum["MSOFT"][42].is_data_line()) result.MqL2_Q=SLHAea::to<double>(spectrum["MSOFT"][42][1]);
365  if (spectrum["MSOFT"][43].is_data_line()) result.MqL3_Q=SLHAea::to<double>(spectrum["MSOFT"][43][1]);
366  if (spectrum["MSOFT"][44].is_data_line()) result.MuR_Q=SLHAea::to<double>(spectrum["MSOFT"][44][1]);
367  if (spectrum["MSOFT"][45].is_data_line()) result.McR_Q=SLHAea::to<double>(spectrum["MSOFT"][45][1]);
368  if (spectrum["MSOFT"][46].is_data_line()) result.MtR_Q=SLHAea::to<double>(spectrum["MSOFT"][46][1]);
369  if (spectrum["MSOFT"][47].is_data_line()) result.MdR_Q=SLHAea::to<double>(spectrum["MSOFT"][47][1]);
370  if (spectrum["MSOFT"][48].is_data_line()) result.MsR_Q=SLHAea::to<double>(spectrum["MSOFT"][48][1]);
371  if (spectrum["MSOFT"][49].is_data_line()) result.MbR_Q=SLHAea::to<double>(spectrum["MSOFT"][49][1]);
372  }
373 
374  if (!spectrum["AU"].empty())
375  {
376  if (spectrum["AU"][1].is_data_line()) result.A_u=SLHAea::to<double>(spectrum["AU"].at(1,1)[2]);
377  if (spectrum["AU"][2].is_data_line()) result.A_c=SLHAea::to<double>(spectrum["AU"].at(2,2)[2]);
378  if (spectrum["AU"][3].is_data_line()) result.A_t=SLHAea::to<double>(spectrum["AU"].at(3,3)[2]);
379  }
380 
381  if (!spectrum["AD"].empty())
382  {
383  if (spectrum["AD"][1].is_data_line()) result.A_d=SLHAea::to<double>(spectrum["AD"].at(1,1)[2]);
384  if (spectrum["AD"][2].is_data_line()) result.A_s=SLHAea::to<double>(spectrum["AD"].at(2,2)[2]);
385  if (spectrum["AD"][3].is_data_line()) result.A_b=SLHAea::to<double>(spectrum["AD"].at(3,3)[2]);
386  }
387 
388  if (!spectrum["AE"].empty())
389  {
390  if (spectrum["AE"][1].is_data_line()) result.A_e=SLHAea::to<double>(spectrum["AE"].at(1,1)[2]);
391  if (spectrum["AE"][2].is_data_line()) result.A_mu=SLHAea::to<double>(spectrum["AE"].at(2,2)[2]);
392  if (spectrum["AE"][3].is_data_line()) result.A_tau=SLHAea::to<double>(spectrum["AE"].at(3,3)[2]);
393  }
394 
395  if (!spectrum["NMSSMRUN"].empty())
396  {
397  if (spectrum["NMSSMRUN"][1].is_data_line()) result.lambdaNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][1][1]);
398  if (spectrum["NMSSMRUN"][2].is_data_line()) result.kappaNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][2][1]);
399  if (spectrum["NMSSMRUN"][3].is_data_line()) result.AlambdaNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][3][1]);
400  if (spectrum["NMSSMRUN"][4].is_data_line()) result.AkappaNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][4][1]);
401  if (spectrum["NMSSMRUN"][5].is_data_line()) result.lambdaSNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][5][1]);
402  if (spectrum["NMSSMRUN"][6].is_data_line()) result.xiFNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][6][1]);
403  if (spectrum["NMSSMRUN"][7].is_data_line()) result.xiSNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][7][1]);
404  if (spectrum["NMSSMRUN"][8].is_data_line()) result.mupNMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][8][1]);
405  if (spectrum["NMSSMRUN"][9].is_data_line()) result.mSp2NMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][9][1]);
406  if (spectrum["NMSSMRUN"][10].is_data_line()) result.mS2NMSSM=SLHAea::to<double>(spectrum["NMSSMRUN"][10][1]);
407  }
408 
409  if (!spectrum["USQMIX"].empty()) for (ie=1;ie<=6;ie++) for (je=1;je<=6;je++)
410  if (spectrum["USQMIX"][max(ie,je)].is_data_line()) result.sU_mix[ie][je]=SLHAea::to<double>(spectrum["USQMIX"].at(ie,je)[2]);
411  if (!spectrum["DSQMIX"].empty()) for (ie=1;ie<=6;ie++) for (je=1;je<=6;je++)
412  if (spectrum["DSQMIX"][max(ie,je)].is_data_line()) result.sD_mix[ie][je]=SLHAea::to<double>(spectrum["DSQMIX"].at(ie,je)[2]);
413  if (!spectrum["SELMIX"].empty()) for (ie=1;ie<=6;ie++) for (je=1;je<=6;je++)
414  if (spectrum["SELMIX"][max(ie,je)].is_data_line()) result.sE_mix[ie][je]=SLHAea::to<double>(spectrum["SELMIX"].at(ie,je)[2]);
415  if (!spectrum["SNUMIX"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
416  if (spectrum["SNUMIX"][max(ie,je)].is_data_line()) result.sNU_mix[ie][je]=SLHAea::to<double>(spectrum["SNUMIX"].at(ie,je)[2]);
417 
418  if (!spectrum["MSQ2"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
419  if (spectrum["MSQ2"][max(ie,je)].is_data_line()) result.sCKM_msq2[ie][je]=SLHAea::to<double>(spectrum["MSQ2"].at(ie,je)[2]);
420  if (!spectrum["MSL2"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
421  if (spectrum["MSL2"][max(ie,je)].is_data_line()) result.sCKM_msl2[ie][je]=SLHAea::to<double>(spectrum["MSL2"].at(ie,je)[2]);
422  if (!spectrum["MSD2"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
423  if (spectrum["MSD2"][max(ie,je)].is_data_line()) result.sCKM_msd2[ie][je]=SLHAea::to<double>(spectrum["MSD2"].at(ie,je)[2]);
424  if (!spectrum["MSU2"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
425  if (spectrum["MSU2"][max(ie,je)].is_data_line()) result.sCKM_msu2[ie][je]=SLHAea::to<double>(spectrum["MSU2"].at(ie,je)[2]);
426  if (!spectrum["MSE2"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
427  if (spectrum["MSE2"][max(ie,je)].is_data_line()) result.sCKM_mse2[ie][je]=SLHAea::to<double>(spectrum["MSE2"].at(ie,je)[2]);
428 
429  if (!spectrum["IMVCKM"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
430  if (spectrum["IMVCKM"][max(ie,je)].is_data_line()) result.IMCKM[ie][je]=SLHAea::to<double>(spectrum["IMVCKM"].at(ie,je)[2]);
431  if (!spectrum["IMVCKM"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
432  if (spectrum["IMVCKM"][max(ie,je)].is_data_line()) result.IMCKM[ie][je]=SLHAea::to<double>(spectrum["IMVCKM"].at(ie,je)[2]);
433 
434  if (!spectrum["UPMNS"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
435  if (spectrum["UPMNS"][max(ie,je)].is_data_line()) result.PMNS_U[ie][je]=SLHAea::to<double>(spectrum["UPMNS"].at(ie,je)[2]);
436 
437  if (!spectrum["TU"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
438  if (spectrum["TU"][max(ie,je)].is_data_line()) result.TU[ie][je]=SLHAea::to<double>(spectrum["TU"].at(ie,je)[2]);
439  if (!spectrum["TD"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
440  if (spectrum["TD"][max(ie,je)].is_data_line()) result.TD[ie][je]=SLHAea::to<double>(spectrum["TD"].at(ie,je)[2]);
441  if (!spectrum["TE"].empty()) for (ie=1;ie<=3;ie++) for (je=1;je<=3;je++)
442  if (spectrum["TE"][max(ie,je)].is_data_line()) result.TE[ie][je]=SLHAea::to<double>(spectrum["TE"].at(ie,je)[2]);
443  }
444 
445  else if (ModelInUse("WC"))
446  {
447  // The Higgs mass doesn't come through in the SLHAea object, as that's only for SLHA2 SM inputs.
448  result.mass_h0 = Dep::SM_spectrum->get(Par::Pole_Mass, "h0_1");
449  // Set the scale.
450  result.Q = result.mass_Z;
451  }
452 
453  BEreq::slha_adjust(&result);
454 
455  // Set the Z and W widths
456  result.width_Z = Dep::Z_decay_rates->width_in_GeV;
457  result.width_W = Dep::W_plus_decay_rates->width_in_GeV;
458 
459  // If requested, override the SuperIso b pole mass with the SpecBit value and recompute the 1S b mass.
460  if (runOptions->getValueOrDef<bool>(false, "take_b_pole_mass_from_spectrum"))
461  {
462  if (ModelInUse("MSSM63atMGUT") or ModelInUse("MSSM63atQ"))
463  {
464  result.mass_h0 = Dep::MSSM_spectrum->get(Par::Pole_Mass, "h0_1");
465  }
466  else if (ModelInUse("WC"))
467  {
468  result.mass_h0 = Dep::SM_spectrum->get(Par::Pole_Mass, "h0_1");
469  }
470  result.mass_b_1S = BEreq::mb_1S(&result);
471  }
472 
473  if (ModelInUse("WC"))
474  {
475 
476  // Tell SuperIso to do its Wilson coefficient calculations for the SM.
477  // We will adjust them with our BSM deviations in backend convenience
478  // functions before we send them to SuperIso's observable calculation functions.
479  result.SM = 1;
480 
481  // So far our model only deals with 5 operators: O_7, O_9, O_10, Q_1 and Q_2.
482  // SuperIso can actually only handle real O_7, O_9 and O_10 too, so the imaginary
483  // parts of those operators get ignored in subsequent calculations.
484  result.Re_DeltaC7 = *Param["Re_DeltaC7"];
485  result.Im_DeltaC7 = *Param["Im_DeltaC7"];
486  result.Re_DeltaC9 = *Param["Re_DeltaC9"];
487  result.Im_DeltaC9 = *Param["Im_DeltaC9"];
488  result.Re_DeltaC10 = *Param["Re_DeltaC10"];
489  result.Im_DeltaC10 = *Param["Im_DeltaC10"];
490  result.Re_DeltaCQ1 = *Param["Re_DeltaCQ1"];
491  result.Im_DeltaCQ1 = *Param["Im_DeltaCQ1"];
492  result.Re_DeltaCQ2 = *Param["Re_DeltaCQ2"];
493  result.Im_DeltaCQ2 = *Param["Im_DeltaCQ2"];
494 
495  }
496 
497  if (flav_debug) cout<<"Finished SI_fill"<<endl;
498  }
499 
500 
502  void SI_bsgamma(double &result)
503  {
504  using namespace Pipes::SI_bsgamma;
505  if (flav_debug) cout<<"Starting SI_bsgamma"<<endl;
506 
507  parameters const& param = *Dep::SuperIso_modelinfo;
508  double E_cut=1.6;
509  result=BEreq::bsgamma_CONV(&param, byVal(E_cut));
510 
511  if (flav_debug) printf("BR(b->s gamma)=%.3e\n",result);
512  if (flav_debug) cout<<"Finished SI_bsgamma"<<endl;
513  }
514 
515 
517  void SI_Bsmumu_untag(double &result)
518  {
519  using namespace Pipes::SI_Bsmumu_untag;
520  if (flav_debug) cout<<"Starting SI_Bsmumu_untag"<<endl;
521 
522  parameters const& param = *Dep::SuperIso_modelinfo;
523  int flav=2;
524  result=BEreq::Bsll_untag_CONV(&param, byVal(flav));
525 
526  if (flav_debug) printf("BR(Bs->mumu)_untag=%.3e\n",result);
527  if (flav_debug) cout<<"Finished SI_Bsmumu_untag"<<endl;
528  }
529 
530 
532  void SI_Bsee_untag(double &result)
533  {
534  using namespace Pipes::SI_Bsee_untag;
535  if (flav_debug) cout<<"Starting SI_Bsee_untag"<<endl;
536 
537  parameters const& param = *Dep::SuperIso_modelinfo;
538  int flav=1;
539  result=BEreq::Bsll_untag_CONV(&param, byVal(flav));
540 
541  if (flav_debug) printf("BR(Bs->ee)_untag=%.3e\n",result);
542  if (flav_debug) cout<<"Finished SI_Bsee_untag"<<endl;
543  }
544 
545 
547  void SI_Bmumu(double &result)
548  {
549  using namespace Pipes::SI_Bmumu;
550  if (flav_debug) cout<<"Starting SI_Bmumu"<<endl;
551 
552  parameters const& param = *Dep::SuperIso_modelinfo;
553  int flav=2;
554  result=BEreq::Bll_CONV(&param, byVal(flav));
555 
556  if (flav_debug) printf("BR(B->mumu)=%.3e\n",result);
557  if (flav_debug) cout<<"Finished SI_Bmumu"<<endl;
558  }
559 
560 
562  void SI_Btaunu(double &result)
563  {
564  using namespace Pipes::SI_Btaunu;
565  if (flav_debug) cout<<"Starting SI_Btaunu"<<endl;
566 
567  parameters const& param = *Dep::SuperIso_modelinfo;
568  result = BEreq::Btaunu(&param);
569 
570  if (flav_debug) printf("BR(B->tau nu)=%.3e\n",result);
571  if (flav_debug) cout<<"Finished SI_Btaunu"<<endl;
572  }
573 
574 
576  void SI_Dstaunu(double &result)
577  {
578  using namespace Pipes::SI_Dstaunu;
579  if (flav_debug) cout<<"Starting SI_Dstaunu"<<endl;
580 
581  parameters const& param = *Dep::SuperIso_modelinfo;
582  result = BEreq::Dstaunu(&param);
583 
584  if (flav_debug) printf("BR(Ds->tau nu)=%.3e\n",result);
585  if (flav_debug) cout<<"Finished SI_Dstaunu"<<endl;
586  }
587 
588 
590  void SI_Dsmunu(double &result)
591  {
592  using namespace Pipes::SI_Dsmunu;
593  if (flav_debug) cout<<"Starting SI_Dsmunu"<<endl;
594 
595  parameters const& param = *Dep::SuperIso_modelinfo;
596  result = BEreq::Dsmunu(&param);
597 
598  if (flav_debug) printf("BR(Ds->mu nu)=%.3e\n",result);
599  if (flav_debug) cout<<"Finished SI_Dsmunu"<<endl;
600  }
601 
602 
604  void SI_Dmunu(double &result)
605  {
606  using namespace Pipes::SI_Dmunu;
607  if (flav_debug) cout<<"Starting SI_Dmunu"<<endl;
608 
609  parameters const& param = *Dep::SuperIso_modelinfo;
610  result = BEreq::Dmunu(&param);
611 
612  if (flav_debug) printf("BR(D->mu nu)=%.3e\n",result);
613  if (flav_debug) cout<<"Finished SI_Dmunu"<<endl;
614  }
615 
616 
618  void SI_BDtaunu(double &result)
619  {
620  using namespace Pipes::SI_BDtaunu;
621  if (flav_debug) cout<<"Starting SI_BDtaunu"<<endl;
622 
623  parameters const& param = *Dep::SuperIso_modelinfo;
624  if (param.model < 0) FlavBit_error().raise(LOCAL_INFO, "Unsupported model.");
625 
626  double q2_min_tau_D = 3.16; // 1.776**2
627  double q2_max_tau_D = 11.6; // (5.28-1.869)**2
628  int gen_tau_D = 3;
629  int charge_tau_D = 0;// D* is the charged version
630  double obs_tau_D[3];
631 
632  result=BEreq::BRBDlnu(byVal(gen_tau_D), byVal( charge_tau_D), byVal(q2_min_tau_D), byVal(q2_max_tau_D), byVal(obs_tau_D), &param);
633 
634  if (flav_debug) printf("BR(B-> D tau nu )=%.3e\n",result);
635  if (flav_debug) cout<<"Finished SI_BDtaunu"<<endl;
636  }
637 
638 
640  void SI_BDmunu(double &result)
641  {
642  using namespace Pipes::SI_BDmunu;
643  if (flav_debug) cout<<"Starting SI_BDmunu"<<endl;
644 
645  parameters const& param = *Dep::SuperIso_modelinfo;
646  if (param.model < 0) FlavBit_error().raise(LOCAL_INFO, "Unsupported model.");
647 
648  double q2_min_mu_D= 0.012; // 0.105*0.105
649  double q2_max_mu_D= 11.6; // (5.28-1.869)**2
650  int gen_mu_D =2;
651  int charge_mu_D =0;// D* is the charged version
652  double obs_mu_D[3];
653 
654  result= BEreq::BRBDlnu(byVal(gen_mu_D), byVal( charge_mu_D), byVal(q2_min_mu_D), byVal(q2_max_mu_D), byVal(obs_mu_D), &param);
655 
656  if (flav_debug) printf("BR(B->D mu nu)=%.3e\n",result);
657  if (flav_debug) cout<<"Finished SI_BDmunu"<<endl;
658  }
659 
660 
662  void SI_BDstartaunu(double &result)
663  {
664  using namespace Pipes::SI_BDstartaunu;
665  if (flav_debug) cout<<"Starting SI_BDstartaunu"<<endl;
666 
667  parameters const& param = *Dep::SuperIso_modelinfo;
668  if (param.model < 0) FlavBit_error().raise(LOCAL_INFO, "Unsupported model.");
669 
670  double q2_min_tau_Dstar = 3.16; // 1.776**2
671  double q2_max_tau_Dstar = 10.67; //(5.279-2.01027)*(5.279-2.01027);
672  int gen_tau_Dstar =3;
673  int charge_tau_Dstar =1;// D* is the charged version
674  double obs_tau_Dstar[3];
675 
676  result= BEreq::BRBDstarlnu(byVal(gen_tau_Dstar), byVal( charge_tau_Dstar), byVal(q2_min_tau_Dstar), byVal(q2_max_tau_Dstar), byVal(obs_tau_Dstar), &param);
677 
678  if (flav_debug) printf("BR(B->Dstar tau nu)=%.3e\n",result);
679  if (flav_debug) cout<<"Finished SI_BDstartaunu"<<endl;
680  }
681 
682 
684  void SI_BDstarmunu(double &result)
685  {
686  using namespace Pipes::SI_BDstarmunu;
687  if (flav_debug) cout<<"Starting SI_BDstarmunu"<<endl;
688 
689  parameters const& param = *Dep::SuperIso_modelinfo;
690  if (param.model < 0) FlavBit_error().raise(LOCAL_INFO, "Unsupported model.");
691 
692  double q2_min_mu_Dstar = 0.012; // 0.105*0.105
693  double q2_max_mu_Dstar = 10.67; //(5.279-2.01027)*(5.279-2.01027);
694  int gen_mu_Dstar =2;
695  int charge_mu_Dstar =1;// D* is the charged version
696  double obs_mu_Dstar[3];
697 
698  result=BEreq::BRBDstarlnu(byVal(gen_mu_Dstar), byVal( charge_mu_Dstar), byVal(q2_min_mu_Dstar), byVal(q2_max_mu_Dstar), byVal(obs_mu_Dstar), &param);
699 
700  if (flav_debug) printf("BR(B->Dstar mu nu)=%.3e\n",result);
701  if (flav_debug) cout<<"Finished SI_BDstarmunu"<<endl;
702  }
703 
704 
706  void SI_RD(double &result)
707  {
708  using namespace Pipes::SI_RD;
709  if (flav_debug) cout<<"Starting SI_RD"<<endl;
710 
711  parameters const& param = *Dep::SuperIso_modelinfo;
712  result = BEreq::BDtaunu_BDenu(&param);
713 
714  if (flav_debug) printf("BR(B->D tau nu)/BR(B->D e nu)=%.3e\n",result);
715  if (flav_debug) cout<<"Finished SI_RD"<<endl;
716  }
717 
718 
720  void SI_RDstar(double &result)
721  {
722  using namespace Pipes::SI_RDstar;
723  if (flav_debug) cout<<"Starting SI_RDstart"<<endl;
724 
725  parameters const& param = *Dep::SuperIso_modelinfo;
726  result = BEreq::BDstartaunu_BDstarenu(&param);
727 
728  if (flav_debug) printf("BR(B->D* tau nu)/BR(B->D* e nu)=%.3e\n",result);
729  if (flav_debug) cout<<"Finished SI_RD*"<<endl;
730  }
731 
732 
734  void SI_Rmu(double &result)
735  {
736  using namespace Pipes::SI_Rmu;
737  if (flav_debug) cout<<"Starting SI_Rmu"<<endl;
738 
739  parameters const& param = *Dep::SuperIso_modelinfo;
740  result = BEreq::Kmunu_pimunu(&param);
741 
742  if (flav_debug) printf("R_mu=BR(K->mu nu)/BR(pi->mu nu)=%.3e\n",result);
743  if (flav_debug) cout<<"Finished SI_Rmu"<<endl;
744  }
745 
746 
748  void SI_Rmu23(double &result)
749  {
750  using namespace Pipes::SI_Rmu23;
751  if (flav_debug) cout<<"Starting SI_Rmu23"<<endl;
752 
753  parameters const& param = *Dep::SuperIso_modelinfo;
754  result = BEreq::Rmu23(&param);
755 
756  if (flav_debug) printf("Rmu23=%.3e\n",result);
757  if (flav_debug) cout<<"Finished SI_Rmu23"<<endl;
758  }
759 
760 
762  void SI_delta0(double &result)
763  {
764  using namespace Pipes::SI_delta0;
765  if (flav_debug) cout<<"Starting SI_delta0"<<endl;
766 
767  parameters const& param = *Dep::SuperIso_modelinfo;
768  result=BEreq::delta0_CONV(&param);
769 
770  if (flav_debug) printf("Delta0(B->K* gamma)=%.3e\n",result);
771  if (flav_debug) cout<<"Finished SI_delta0"<<endl;
772  }
773 
774 
776  void SI_BRBXsmumu_lowq2(double &result)
777  {
778  using namespace Pipes::SI_BRBXsmumu_lowq2;
779  if (flav_debug) cout<<"Starting SI_BRBXsmumu_lowq2"<<endl;
780 
781  parameters const& param = *Dep::SuperIso_modelinfo;
782  result=BEreq::BRBXsmumu_lowq2_CONV(&param);
783 
784  if (flav_debug) printf("BR(B->Xs mu mu)_lowq2=%.3e\n",result);
785  if (flav_debug) cout<<"Finished SI_BRBXsmumu_lowq2"<<endl;
786  }
787 
788 
790  void SI_BRBXsmumu_highq2(double &result)
791  {
792  using namespace Pipes::SI_BRBXsmumu_highq2;
793  if (flav_debug) cout<<"Starting SI_BRBXsmumu_highq2"<<endl;
794 
795  parameters const& param = *Dep::SuperIso_modelinfo;
796  result=BEreq::BRBXsmumu_highq2_CONV(&param);
797 
798  if (flav_debug) printf("BR(B->Xs mu mu)_highq2=%.3e\n",result);
799  if (flav_debug) cout<<"Finished SI_BRBXsmumu_highq2"<<endl;
800  }
801 
802 
804  void SI_A_BXsmumu_lowq2(double &result)
805  {
806  using namespace Pipes::SI_A_BXsmumu_lowq2;
807  if (flav_debug) cout<<"Starting SI_A_BXsmumu_lowq2"<<endl;
808 
809  parameters const& param = *Dep::SuperIso_modelinfo;
810  result=BEreq::A_BXsmumu_lowq2_CONV(&param);
811 
812  if (flav_debug) printf("AFB(B->Xs mu mu)_lowq2=%.3e\n",result);
813  if (flav_debug) cout<<"Finished SI_A_BXsmumu_lowq2"<<endl;
814  }
815 
816 
818  void SI_A_BXsmumu_highq2(double &result)
819  {
820  using namespace Pipes::SI_A_BXsmumu_highq2;
821  if (flav_debug) cout<<"Starting SI_A_BXsmumu_highq2"<<endl;
822 
823  parameters const& param = *Dep::SuperIso_modelinfo;
824  result=BEreq::A_BXsmumu_highq2_CONV(&param);
825 
826  if (flav_debug) printf("AFB(B->Xs mu mu)_highq2=%.3e\n",result);
827  if (flav_debug) cout<<"Finished SI_A_BXsmumu_highq2"<<endl;
828  }
829 
830 
832  void SI_A_BXsmumu_zero(double &result)
833  {
834  using namespace Pipes::SI_A_BXsmumu_zero;
835  if (flav_debug) cout<<"Starting SI_A_BXsmumu_zero"<<endl;
836 
837  parameters const& param = *Dep::SuperIso_modelinfo;
838  result=BEreq::A_BXsmumu_zero_CONV(&param);
839 
840  if (flav_debug) printf("AFB(B->Xs mu mu)_zero=%.3e\n",result);
841  if (flav_debug) cout<<"Finished SI_A_BXsmumu_zero"<<endl;
842  }
843 
844 
846  void SI_BRBXstautau_highq2(double &result)
847  {
848  using namespace Pipes::SI_BRBXstautau_highq2;
849  if (flav_debug) cout<<"Starting SI_BRBXstautau_highq2"<<endl;
850 
851  parameters const& param = *Dep::SuperIso_modelinfo;
852  result=BEreq::BRBXstautau_highq2_CONV(&param);
853 
854  if (flav_debug) printf("BR(B->Xs tau tau)_highq2=%.3e\n",result);
855  if (flav_debug) cout<<"Finished SI_BRBXstautau_highq2"<<endl;
856  }
857 
858 
860  void SI_A_BXstautau_highq2(double &result)
861  {
862  using namespace Pipes::SI_A_BXstautau_highq2;
863  if (flav_debug) cout<<"Starting SI_A_BXstautau_highq2"<<endl;
864 
865  parameters const& param = *Dep::SuperIso_modelinfo;
866  result=BEreq::A_BXstautau_highq2_CONV(&param);
867 
868  if (flav_debug) printf("AFB(B->Xs tau tau)_highq2=%.3e\n",result);
869  if (flav_debug) cout<<"Finished SI_A_BXstautau_highq2"<<endl;
870  }
871 
872 
875  #define DEFINE_BKSTARMUMU(Q2MIN, Q2MAX, Q2MIN_TAG, Q2MAX_TAG) \
876  void CAT_4(SI_BKstarmumu_,Q2MIN_TAG,_,Q2MAX_TAG)(Flav_KstarMuMu_obs &result) \
877  { \
878  using namespace Pipes::CAT_4(SI_BKstarmumu_,Q2MIN_TAG,_,Q2MAX_TAG); \
879  if (flav_debug) cout<<"Starting " STRINGIFY(CAT_4(SI_BKstarmumu_,Q2MIN_TAG,_,Q2MAX_TAG))<<endl; \
880  parameters const& param = *Dep::SuperIso_modelinfo; \
881  result=BEreq::BKstarmumu_CONV(&param, Q2MIN, Q2MAX); \
882  if (flav_debug) cout<<"Finished " STRINGIFY(CAT_4(SI_BKstarmumu_,Q2MIN_TAG,_,Q2MAX_TAG))<<endl; \
883  }
884  DEFINE_BKSTARMUMU(1.1, 2.5, 11, 25)
885  DEFINE_BKSTARMUMU(2.5, 4.0, 25, 40)
886  DEFINE_BKSTARMUMU(4.0, 6.0, 40, 60)
887  DEFINE_BKSTARMUMU(6.0, 8.0, 60, 80)
888  DEFINE_BKSTARMUMU(15., 17., 15, 17)
889  DEFINE_BKSTARMUMU(17., 19., 17, 19)
891  #undef DEFINE_BKSTARMUMU
892 
894  void SI_RKstar_0045_11(double &result)
895  {
896  using namespace Pipes::SI_RKstar_0045_11;
897  if (flav_debug) cout<<"Starting SI_RKstar_0045_11"<<endl;
898 
899  parameters const& param = *Dep::SuperIso_modelinfo;
900  result=BEreq::RKstar_CONV(&param,0.045,1.1);
901 
902  if (flav_debug) printf("RK*_lowq2=%.3e\n",result);
903  if (flav_debug) cout<<"Finished SI_RKstar_0045_11"<<endl;
904  }
905 
907  void SI_RKstar_11_60(double &result)
908  {
909  using namespace Pipes::SI_RKstar_11_60;
910  if (flav_debug) cout<<"Starting SI_RKstar_11_60"<<endl;
911 
912  parameters const& param = *Dep::SuperIso_modelinfo;
913  result=BEreq::RKstar_CONV(&param,1.1,6.0);
914 
915  if (flav_debug) printf("RK*_intermq2=%.3e\n",result);
916  if (flav_debug) cout<<"Finished SI_RKstar_11_60"<<endl;
917  }
918 
919  // RK* for RHN, using same approximations as RK, low q^2
920  void RHN_RKstar_0045_11(double &result)
921  {
922  using namespace Pipes::RHN_RKstar_0045_11;
923  SMInputs sminputs = *Dep::SMINPUTS;
924  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
925  std::vector<double> mN = {*Param["M_1"],*Param["M_2"],*Param["M_3"]};
926  double mt = *Param["mT"];
927 
928  if (flav_debug) cout << "Starting RHN_RKstar_0045_11" << endl;
929 
930  const double mW = sminputs.mW;
931  const double sinW2 = sqrt(1.0 - pow(sminputs.mW/sminputs.mZ,2));
932 
933  // NNLL calculation of SM Wilson coefficients from 1712.01593 and 0811.1214
934  const double C10_SM = -4.103;
935  const double C9_SM = 4.211;
936 
937  // Wilson coefficients for the RHN model, from 1706.07570
938  std::complex<double> C10_mu = {0.0, 0.0}, C10_e = {0.0, 0.0};
939  for(int i=0; i<3; i++)
940  {
941  C10_mu += 1.0/(4.0*sinW2)*Theta.adjoint()(i,1)*Theta(1,i) * LoopFunctions::E(pow(mt/mW,2),pow(mN[i]/mW,2));
942  C10_e += 1.0/(4.0*sinW2)*Theta.adjoint()(i,0)*Theta(0,i) * LoopFunctions::E(pow(mt/mW,2),pow(mN[i]/mW,2));
943  }
944  std::complex<double> C9_mu = - C10_mu, C9_e = -C10_e;
945 
946  // Aproximated solution from eq A.3 in 1408.4097
947  result = std::norm(C10_SM + C10_mu) + std::norm(C9_SM + C9_mu);
948  result /= std::norm(C10_SM + C10_e) + std::norm(C9_SM + C9_e);
949 
950  if (flav_debug) cout << "RK = " << result << endl;
951  if (flav_debug) cout << "Finished RHN_RKstar_0045_11" << endl;
952 
953  }
954 
955  // RK* for RHN, using same approximations as RK, intermediate q^2
956  void RHN_RKstar_11_60(double &result)
957  {
958  using namespace Pipes::RHN_RKstar_11_60;
959  SMInputs sminputs = *Dep::SMINPUTS;
960  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
961  std::vector<double> mN = {*Param["M_1"],*Param["M_2"],*Param["M_3"]};
962  double mt = *Param["mT"];
963 
964  if (flav_debug) cout << "Starting RHN_RKstar_11_60" << endl;
965 
966  const double mW = sminputs.mW;
967  const double sinW2 = sqrt(1.0 - pow(sminputs.mW/sminputs.mZ,2));
968 
969  // NNLL calculation of SM Wilson coefficients from 1712.01593 and 0811.1214
970  const double C10_SM = -4.103;
971  const double C9_SM = 4.211;
972 
973  // Wilson coefficients for the RHN model, from 1706.07570
974  std::complex<double> C10_mu = {0.0, 0.0}, C10_e = {0.0, 0.0};
975  for(int i=0; i<3; i++)
976  {
977  C10_mu += 1.0/(4.0*sinW2)*Theta.adjoint()(i,1)*Theta(1,i) * LoopFunctions::E(pow(mt/mW,2),pow(mN[i]/mW,2));
978  C10_e += 1.0/(4.0*sinW2)*Theta.adjoint()(i,0)*Theta(0,i) * LoopFunctions::E(pow(mt/mW,2),pow(mN[i]/mW,2));
979  }
980  std::complex<double> C9_mu = - C10_mu, C9_e = -C10_e;
981 
982  // Aproximated solution from eq A.3 in 1408.4097
983  result = std::norm(C10_SM + C10_mu) + std::norm(C9_SM + C9_mu);
984  result /= std::norm(C10_SM + C10_e) + std::norm(C9_SM + C9_e);
985 
986  if (flav_debug) cout << "RK = " << result << endl;
987  if (flav_debug) cout << "Finished RHN_RKstar_11_60" << endl;
988 
989  }
990 
992  void SI_RK(double &result)
993  {
994  using namespace Pipes::SI_RK;
995  if (flav_debug) cout<<"Starting SI_RK"<<endl;
996 
997  parameters const& param = *Dep::SuperIso_modelinfo;
998  result=BEreq::RK_CONV(&param,1.0,6.0);
999 
1000  if (flav_debug) printf("RK=%.3e\n",result);
1001  if (flav_debug) cout<<"Finished SI_RK"<<endl;
1002  }
1003 
1005  void RHN_RK(double &result)
1006  {
1007  using namespace Pipes::RHN_RK;
1008  SMInputs sminputs = *Dep::SMINPUTS;
1009  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1010  std::vector<double> mN = {*Param["M_1"],*Param["M_2"],*Param["M_3"]};
1011  double mt = *Param["mT"];
1012 
1013  if (flav_debug) cout << "Starting RHN_RK" << endl;
1014 
1015  const double mW = sminputs.mW;
1016  const double sinW2 = sqrt(1.0 - pow(sminputs.mW/sminputs.mZ,2));
1017 
1018  // NNLL calculation of SM Wilson coefficients from 1712.01593 and 0811.1214
1019  const double C10_SM = -4.103;
1020  const double C9_SM = 4.211;
1021 
1022  // Wilson coefficients for the RHN model, from 1706.07570
1023  std::complex<double> C10_mu = {0.0, 0.0}, C10_e = {0.0, 0.0};
1024  for(int i=0; i<3; i++)
1025  {
1026  C10_mu += 1.0/(4.0*sinW2)*Theta.adjoint()(i,1)*Theta(1,i) * LoopFunctions::E(pow(mt/mW,2),pow(mN[i]/mW,2));
1027  C10_e += 1.0/(4.0*sinW2)*Theta.adjoint()(i,0)*Theta(0,i) * LoopFunctions::E(pow(mt/mW,2),pow(mN[i]/mW,2));
1028  }
1029  std::complex<double> C9_mu = - C10_mu, C9_e = -C10_e;
1030 
1031  // Aproximated solution from eq A.3 in 1408.4097
1032  result = std::norm(C10_SM + C10_mu) + std::norm(C9_SM + C9_mu);
1033  result /= std::norm(C10_SM + C10_e) + std::norm(C9_SM + C9_e);
1034 
1035  if (flav_debug) cout << "RK = " << result << endl;
1036  if (flav_debug) cout << "Finished RHN_RK" << endl;
1037  }
1038 
1040  void SI_AI_BKstarmumu(double &result)
1041  {
1042  using namespace Pipes::SI_AI_BKstarmumu;
1043  if (flav_debug) cout<<"Starting SI_AI_BKstarmumu"<<endl;
1044 
1045  parameters const& param = *Dep::SuperIso_modelinfo;
1046  result=BEreq::SI_AI_BKstarmumu_CONV(&param);
1047 
1048  if (flav_debug) printf("A_I(B->K* mu mu)_lowq2=%.3e\n",result);
1049  if (flav_debug) cout<<"Finished SI_AI_BKstarmumu"<<endl;
1050  }
1051 
1052 
1054  void SI_AI_BKstarmumu_zero(double &result)
1055  {
1056  using namespace Pipes::SI_AI_BKstarmumu_zero;
1057 
1058  if (flav_debug) cout<<"Starting SI_AI_BKstarmumu_zero"<<endl;
1059 
1060  parameters const& param = *Dep::SuperIso_modelinfo;
1061  result=BEreq::SI_AI_BKstarmumu_zero_CONV(&param);
1062 
1063  if (flav_debug) printf("A_I(B->K* mu mu)_zero=%.3e\n",result);
1064  if (flav_debug) cout<<"Finished SI_AI_BKstarmumu_zero"<<endl;
1065  }
1066 
1067 
1069  void FH_FlavourObs(fh_FlavourObs &result)
1070  {
1071  using namespace Pipes::FH_FlavourObs;
1072 
1073  if (flav_debug) cout<<"Starting FH_FlavourObs"<<endl;
1074 
1075  fh_real bsgMSSM; // B -> Xs gamma in MSSM
1076  fh_real bsgSM; // B -> Xs gamma in SM
1077  fh_real deltaMsMSSM; // delta Ms in MSSM
1078  fh_real deltaMsSM; // delta Ms in SM
1079  fh_real bsmumuMSSM; // Bs -> mu mu in MSSM
1080  fh_real bsmumuSM; // Bs -> mu mu in SM
1081 
1082  int error = 1;
1083  BEreq::FHFlavour(error, bsgMSSM, bsgSM,
1084  deltaMsMSSM, deltaMsSM,
1085  bsmumuMSSM, bsmumuSM);
1086 
1087  fh_FlavourObs FlavourObs;
1088  FlavourObs.Bsg_MSSM = bsgMSSM;
1089  FlavourObs.Bsg_SM = bsgSM;
1090  FlavourObs.deltaMs_MSSM = deltaMsMSSM;
1091  FlavourObs.deltaMs_SM = deltaMsSM;
1092  FlavourObs.Bsmumu_MSSM = bsmumuMSSM;
1093  FlavourObs.Bsmumu_SM = bsmumuSM;
1094 
1095  result = FlavourObs;
1096  if (flav_debug) cout<<"Finished FH_FlavourObs"<<endl;
1097  }
1098 
1099 
1102  void FH_bsgamma(double &result)
1103  {
1104  result = Pipes::FH_bsgamma::Dep::FH_FlavourObs->Bsg_MSSM;
1105  }
1106  void FH_Bsmumu (double &result)
1107  {
1108  result = Pipes::FH_Bsmumu::Dep::FH_FlavourObs->Bsmumu_MSSM;
1109  }
1110  void FH_DeltaMs(double &result)
1111  {
1112  result = Pipes::FH_DeltaMs::Dep::FH_FlavourObs->deltaMs_MSSM;
1113  }
1115 
1116 
1119  {
1120  using namespace Pipes::b2sll_measurements;
1121 
1122  static bool first = true;
1123  static int n_experiments;
1124 
1125  if (flav_debug) cout<<"Starting b2sll_measurements function"<<endl;
1126  if (flav_debug and first) cout <<"Initialising Flav Reader in b2sll_measurements"<<endl;
1127 
1128  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
1129  if (first)
1130  {
1131  pmc.LL_name="b2sll_likelihood";
1132 
1133  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
1134  fread.debug_mode(flav_debug);
1135 
1136  const vector<string> observablesn = {"FL", "AFB", "S3", "S4", "S5", "S7", "S8", "S9"};
1137  const vector<string> observablesq = {"1.1-2.5", "2.5-4", "4-6", "6-8", "15-17", "17-19"};
1138  vector<string> observables;
1139  for (unsigned i=0;i<observablesq.size();++i)
1140  {
1141  for (unsigned j=0;j<observablesn.size();++j)
1142  {
1143  observables.push_back(observablesn[j]+"_B0Kstar0mumu_"+observablesq[i]);
1144  }
1145  }
1146 
1147  for (unsigned i=0;i<observables.size();++i)
1148  {
1149  fread.read_yaml_measurement("flav_data.yaml", observables[i]);
1150  }
1151 
1152  fread.initialise_matrices();
1153  pmc.cov_exp = fread.get_exp_cov();
1154  pmc.value_exp = fread.get_exp_value();
1155  pmc.cov_th = Kstarmumu_theory_err().get_th_cov(observables);
1156  n_experiments = pmc.cov_th.size1();
1157  pmc.value_th.resize(n_experiments,1);
1158  pmc.dim=n_experiments;
1159 
1160  // We assert that the experiments and the observables are the same size
1161  assert(pmc.value_exp.size1() == observables.size());
1162 
1163  // Init out.
1164  first = false;
1165  }
1166 
1167  pmc.value_th(0,0)=Dep::BKstarmumu_11_25->FL;
1168  pmc.value_th(1,0)=Dep::BKstarmumu_11_25->AFB;
1169  pmc.value_th(2,0)=Dep::BKstarmumu_11_25->S3;
1170  pmc.value_th(3,0)=Dep::BKstarmumu_11_25->S4;
1171  pmc.value_th(4,0)=Dep::BKstarmumu_11_25->S5;
1172  pmc.value_th(5,0)=Dep::BKstarmumu_11_25->S7;
1173  pmc.value_th(6,0)=Dep::BKstarmumu_11_25->S8;
1174  pmc.value_th(7,0)=Dep::BKstarmumu_11_25->S9;
1175 
1176  pmc.value_th(8,0)=Dep::BKstarmumu_25_40->FL;
1177  pmc.value_th(9,0)=Dep::BKstarmumu_25_40->AFB;
1178  pmc.value_th(10,0)=Dep::BKstarmumu_25_40->S3;
1179  pmc.value_th(11,0)=Dep::BKstarmumu_25_40->S4;
1180  pmc.value_th(12,0)=Dep::BKstarmumu_25_40->S5;
1181  pmc.value_th(13,0)=Dep::BKstarmumu_25_40->S7;
1182  pmc.value_th(14,0)=Dep::BKstarmumu_25_40->S8;
1183  pmc.value_th(15,0)=Dep::BKstarmumu_25_40->S9;
1184 
1185  pmc.value_th(16,0)=Dep::BKstarmumu_40_60->FL;
1186  pmc.value_th(17,0)=Dep::BKstarmumu_40_60->AFB;
1187  pmc.value_th(18,0)=Dep::BKstarmumu_40_60->S3;
1188  pmc.value_th(19,0)=Dep::BKstarmumu_40_60->S4;
1189  pmc.value_th(20,0)=Dep::BKstarmumu_40_60->S5;
1190  pmc.value_th(21,0)=Dep::BKstarmumu_40_60->S7;
1191  pmc.value_th(22,0)=Dep::BKstarmumu_40_60->S8;
1192  pmc.value_th(23,0)=Dep::BKstarmumu_40_60->S9;
1193 
1194  pmc.value_th(24,0)=Dep::BKstarmumu_60_80->FL;
1195  pmc.value_th(25,0)=Dep::BKstarmumu_60_80->AFB;
1196  pmc.value_th(26,0)=Dep::BKstarmumu_60_80->S3;
1197  pmc.value_th(27,0)=Dep::BKstarmumu_60_80->S4;
1198  pmc.value_th(28,0)=Dep::BKstarmumu_60_80->S5;
1199  pmc.value_th(29,0)=Dep::BKstarmumu_60_80->S7;
1200  pmc.value_th(30,0)=Dep::BKstarmumu_60_80->S8;
1201  pmc.value_th(31,0)=Dep::BKstarmumu_60_80->S9;
1202 
1203  pmc.value_th(32,0)=Dep::BKstarmumu_15_17->FL;
1204  pmc.value_th(33,0)=Dep::BKstarmumu_15_17->AFB;
1205  pmc.value_th(34,0)=Dep::BKstarmumu_15_17->S3;
1206  pmc.value_th(35,0)=Dep::BKstarmumu_15_17->S4;
1207  pmc.value_th(36,0)=Dep::BKstarmumu_15_17->S5;
1208  pmc.value_th(37,0)=Dep::BKstarmumu_15_17->S7;
1209  pmc.value_th(38,0)=Dep::BKstarmumu_15_17->S8;
1210  pmc.value_th(39,0)=Dep::BKstarmumu_15_17->S9;
1211 
1212  pmc.value_th(40,0)=Dep::BKstarmumu_17_19->FL;
1213  pmc.value_th(41,0)=Dep::BKstarmumu_17_19->AFB;
1214  pmc.value_th(42,0)=Dep::BKstarmumu_17_19->S3;
1215  pmc.value_th(43,0)=Dep::BKstarmumu_17_19->S4;
1216  pmc.value_th(44,0)=Dep::BKstarmumu_17_19->S5;
1217  pmc.value_th(45,0)=Dep::BKstarmumu_17_19->S7;
1218  pmc.value_th(46,0)=Dep::BKstarmumu_17_19->S8;
1219  pmc.value_th(47,0)=Dep::BKstarmumu_17_19->S9;
1220 
1221  pmc.diff.clear();
1222  for (int i=0;i<n_experiments;++i)
1223  {
1224  pmc.diff.push_back(pmc.value_exp(i,0)-pmc.value_th(i,0));
1225  }
1226 
1227  if (flav_debug) cout<<"Finished b2sll_measurements function"<<endl;
1228  }
1229 
1230 
1232  void b2sll_likelihood(double &result)
1233  {
1234  using namespace Pipes::b2sll_likelihood;
1235 
1236  if (flav_debug) cout<<"Starting b2sll_likelihood"<<endl;
1237 
1238  // Get experimental measurements
1239  predictions_measurements_covariances pmc=*Dep::b2sll_M;
1240 
1241  // Get experimental covariance
1242  boost::numeric::ublas::matrix<double> cov=pmc.cov_exp;
1243 
1244  // adding theory and experimenta covariance
1245  cov+=pmc.cov_th;
1246 
1247  //calculating a diff
1248  vector<double> diff;
1249  diff=pmc.diff;
1250  boost::numeric::ublas::matrix<double> cov_inv(pmc.dim, pmc.dim);
1251  InvertMatrix(cov, cov_inv);
1252 
1253  double Chi2=0;
1254 
1255  for (int i=0; i < pmc.dim; ++i)
1256  {
1257  for (int j=0; j<pmc.dim; ++j)
1258  {
1259  Chi2+= diff[i] * cov_inv(i,j)*diff[j] ;
1260  }
1261  }
1262 
1263  result=-0.5*Chi2;
1264 
1265  if (flav_debug) cout<<"Finished b2sll_likelihood"<<endl;
1266  if (flav_debug_LL) cout<<"Likelihood result b2sll_likelihood : "<< result<<endl;
1267 
1268  }
1269 
1270 
1272  void deltaMB_likelihood(double &result)
1273  {
1274  using namespace Pipes::deltaMB_likelihood;
1275  static bool th_err_absolute, first = true;
1276  static double exp_meas, exp_DeltaMs_err, th_err;
1277 
1278  if (flav_debug) cout << "Starting Delta_Ms_likelihood"<<endl;
1279 
1280  if (first)
1281  {
1282  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
1283  fread.debug_mode(flav_debug);
1284  if (flav_debug) cout<<"Initialised Flav reader in Delta_Ms_likelihood"<<endl;
1285  fread.read_yaml_measurement("flav_data.yaml", "DeltaMs");
1286  fread.initialise_matrices(); // here we have a single measurement ;) so let's be sneaky:
1287  exp_meas = fread.get_exp_value()(0,0);
1288  exp_DeltaMs_err = sqrt(fread.get_exp_cov()(0,0));
1289  th_err = fread.get_th_err()(0,0).first;
1290  th_err_absolute = fread.get_th_err()(0,0).second;
1291  first = false;
1292  }
1293 
1294  if (flav_debug) cout << "Experiment: " << exp_meas << " " << exp_DeltaMs_err << " " << th_err << endl;
1295 
1296  // Now we do the stuff that actually depends on the parameters
1297  double theory_prediction = *Dep::DeltaMs;
1298  double theory_DeltaMs_err = th_err * (th_err_absolute ? 1.0 : std::abs(theory_prediction));
1299  if (flav_debug) cout<<"Theory prediction: "<<theory_prediction<<" +/- "<<theory_DeltaMs_err<<endl;
1300 
1302  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
1303 
1304  result = Stats::gaussian_loglikelihood(theory_prediction, exp_meas, theory_DeltaMs_err, exp_DeltaMs_err, profile);
1305  }
1306 
1307 
1309  void b2sgamma_likelihood(double &result)
1310  {
1311  using namespace Pipes::b2sgamma_likelihood;
1312 
1313  static bool th_err_absolute, first = true;
1314  static double exp_meas, exp_b2sgamma_err, th_err;
1315 
1316  if (flav_debug) cout << "Starting b2sgamma_measurements"<<endl;
1317 
1318  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
1319  if (first)
1320  {
1321  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
1322  fread.debug_mode(flav_debug);
1323  if (flav_debug) cout<<"Initialised Flav reader in b2sgamma_measurements"<<endl;
1324  fread.read_yaml_measurement("flav_data.yaml", "BR_b2sgamma");
1325  fread.initialise_matrices(); // here we have a single measurement ;) so let's be sneaky:
1326  exp_meas = fread.get_exp_value()(0,0);
1327  exp_b2sgamma_err = sqrt(fread.get_exp_cov()(0,0));
1328  th_err = fread.get_th_err()(0,0).first;
1329  th_err_absolute = fread.get_th_err()(0,0).second;
1330  first = false;
1331  }
1332 
1333  if (flav_debug) cout << "Experiment: " << exp_meas << " " << exp_b2sgamma_err << " " << th_err << endl;
1334 
1335  // Now we do the stuff that actually depends on the parameters
1336  double theory_prediction = *Dep::bsgamma;
1337  double theory_b2sgamma_err = th_err * (th_err_absolute ? 1.0 : std::abs(theory_prediction));
1338  if (flav_debug) cout<<"Theory prediction: "<<theory_prediction<<" +/- "<<theory_b2sgamma_err<<endl;
1339 
1341  bool profile = runOptions->getValueOrDef<bool>(false, "profile_systematics");
1342 
1343  result = Stats::gaussian_loglikelihood(theory_prediction, exp_meas, theory_b2sgamma_err, exp_b2sgamma_err, profile);
1344  }
1345 
1346 
1349  {
1350  using namespace Pipes::b2ll_measurements;
1351 
1352  static bool bs2mumu_err_absolute, b2mumu_err_absolute, first = true;
1353  static double theory_bs2mumu_err, theory_b2mumu_err;
1354 
1355  if (flav_debug) cout<<"Starting b2ll_measurements"<<endl;
1356 
1357  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
1358  if (first)
1359  {
1360  pmc.LL_name="b2ll_likelihood";
1361 
1362  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
1363  fread.debug_mode(flav_debug);
1364 
1365  if (flav_debug) cout<<"Initiated Flav reader in b2ll_measurements"<<endl;
1366  fread.read_yaml_measurement("flav_data.yaml", "BR_Bs2mumu");
1367  fread.read_yaml_measurement("flav_data.yaml", "BR_B02mumu");
1368  if (flav_debug) cout<<"Finished reading b->mumu data"<<endl;
1369 
1370  fread.initialise_matrices();
1371 
1372  theory_bs2mumu_err = fread.get_th_err()(0,0).first;
1373  theory_b2mumu_err = fread.get_th_err()(1,0).first;
1374  bs2mumu_err_absolute = fread.get_th_err()(0,0).second;
1375  b2mumu_err_absolute = fread.get_th_err()(1,0).second;
1376 
1377  pmc.value_exp=fread.get_exp_value();
1378  pmc.cov_exp=fread.get_exp_cov();
1379 
1380  pmc.value_th.resize(2,1);
1381  pmc.cov_th.resize(2,2);
1382 
1383  pmc.dim=2;
1384 
1385  // Init over and out.
1386  first = false;
1387  }
1388 
1389  // Get theory prediction
1390  pmc.value_th(0,0)=*Dep::Bsmumu_untag;
1391  pmc.value_th(1,0)=*Dep::Bmumu;
1392 
1393  // Compute error on theory prediction and populate the covariance matrix
1394  double theory_bs2mumu_error = theory_bs2mumu_err * (bs2mumu_err_absolute ? 1.0 : *Dep::Bsmumu_untag);
1395  double theory_b2mumu_error = theory_b2mumu_err * (b2mumu_err_absolute ? 1.0 : *Dep::Bmumu);
1396  pmc.cov_th(0,0)=theory_bs2mumu_error*theory_bs2mumu_error;
1397  pmc.cov_th(0,1)=0.;
1398  pmc.cov_th(1,0)=0.;
1399  pmc.cov_th(1,1)=theory_b2mumu_error*theory_b2mumu_error;
1400 
1401  // Save the differences between theory and experiment
1402  pmc.diff.clear();
1403  for (int i=0;i<2;++i)
1404  {
1405  pmc.diff.push_back(pmc.value_exp(i,0)-pmc.value_th(i,0));
1406  }
1407 
1408  if (flav_debug) cout<<"Finished b2ll_measurements"<<endl;
1409 
1410  }
1411 
1412 
1414  void b2ll_likelihood(double &result)
1415  {
1416  using namespace Pipes::b2ll_likelihood;
1417 
1418  if (flav_debug) cout<<"Starting b2ll_likelihood"<<endl;
1419 
1420  predictions_measurements_covariances pmc = *Dep::b2ll_M;
1421 
1422  boost::numeric::ublas::matrix<double> cov=pmc.cov_exp;
1423 
1424  // adding theory and experimental covariance
1425  cov+=pmc.cov_th;
1426 
1427  //calculating a diff
1428  vector<double> diff;
1429  diff=pmc.diff;
1430 
1431  boost::numeric::ublas::matrix<double> cov_inv(pmc.dim, pmc.dim);
1432  InvertMatrix(cov, cov_inv);
1433 
1434  // calculating the chi2
1435  double Chi2=0;
1436 
1437  for (int i=0; i < pmc.dim; ++i)
1438  {
1439  for (int j=0; j<pmc.dim; ++j)
1440  {
1441  Chi2+= diff[i] * cov_inv(i,j)*diff[j];
1442  }
1443  }
1444 
1445  result=-0.5*Chi2;
1446 
1447  if (flav_debug) cout<<"Finished b2ll_likelihood"<<endl;
1448  if (flav_debug_LL) cout<<"Likelihood result b2ll_likelihood : "<< result<<endl;
1449 
1450  }
1451 
1452 
1455  {
1456  using namespace Pipes::SL_measurements;
1457 
1458  const int n_experiments=8;
1459  static bool th_err_absolute[n_experiments], first = true;
1460  static double th_err[n_experiments];
1461 
1462  if (flav_debug) cout<<"Starting SL_measurements"<<endl;
1463 
1464  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
1465  if (first)
1466  {
1467  pmc.LL_name="SL_likelihood";
1468 
1469  // Read in experimental measuremens
1470  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
1471  fread.debug_mode(flav_debug);
1472  if (flav_debug) cout<<"Initialised Flav reader in SL_measurements"<<endl;
1473 
1474  // B-> tau nu
1475  fread.read_yaml_measurement("flav_data.yaml", "BR_Btaunu");
1476  // B-> D mu nu
1477  fread.read_yaml_measurement("flav_data.yaml", "BR_BDmunu");
1478  // B-> D* mu nu
1479  fread.read_yaml_measurement("flav_data.yaml", "BR_BDstarmunu");
1480  // RD
1481  fread.read_yaml_measurement("flav_data.yaml", "RD");
1482  // RDstar
1483  fread.read_yaml_measurement("flav_data.yaml", "RDstar");
1484  // Ds-> tau nu
1485  fread.read_yaml_measurement("flav_data.yaml", "BR_Dstaunu");
1486  // Ds -> mu nu
1487  fread.read_yaml_measurement("flav_data.yaml", "BR_Dsmunu");
1488  // D -> mu nu
1489  fread.read_yaml_measurement("flav_data.yaml", "BR_Dmunu");
1490 
1491  fread.initialise_matrices();
1492  pmc.cov_exp=fread.get_exp_cov();
1493  pmc.value_exp=fread.get_exp_value();
1494 
1495  pmc.value_th.resize(n_experiments,1);
1496  // Set all entries in the covariance matrix explicitly to zero, as we will only write the diagonal ones later.
1497  pmc.cov_th = boost::numeric::ublas::zero_matrix<double>(n_experiments,n_experiments);
1498  for (int i = 0; i < n_experiments; ++i)
1499  {
1500  th_err[i] = fread.get_th_err()(i,0).first;
1501  th_err_absolute[i] = fread.get_th_err()(i,0).second;
1502  }
1503 
1504  pmc.dim=n_experiments;
1505 
1506  // Init over.
1507  first = false;
1508  }
1509 
1510  // R(D) is calculated assuming isospin symmetry
1511  double theory[8];
1512  // B-> tau nu SI
1513  theory[0] = *Dep::Btaunu;
1514  // B-> D mu nu
1515  theory[1] = *Dep::BDmunu;
1516  // B-> D* mu nu
1517  theory[2] = *Dep::BDstarmunu;
1518  // RD
1519  theory[3] = *Dep::RD;
1520  // RDstar
1521  theory[4] = *Dep::RDstar;
1522  // Ds-> tau nu
1523  theory[5] = *Dep::Dstaunu;
1524  // Ds -> mu nu
1525  theory[6] = *Dep::Dsmunu;
1526  // D -> mu nu
1527  theory[7] =*Dep::Dmunu;
1528 
1529  for (int i = 0; i < n_experiments; ++i)
1530  {
1531  pmc.value_th(i,0) = theory[i];
1532  pmc.cov_th(i,i) = th_err[i]*th_err[i] * (th_err_absolute[i] ? 1.0 : theory[i]*theory[i]);
1533  }
1534  // Add in the correlations between B-> D mu nu and RD
1535  pmc.cov_th(1,3) = pmc.cov_th(3,1) = -0.55 * th_err[1]*th_err[3] * (th_err_absolute[1] ? 1.0 : theory[1]) * (th_err_absolute[3] ? 1.0 : theory[3]);
1536  // Add in the correlations between B-> D* mu nu and RD*
1537  pmc.cov_th(2,4) = pmc.cov_th(4,2) = -0.62 * th_err[2]*th_err[4] * (th_err_absolute[2] ? 1.0 : theory[2]) * (th_err_absolute[4] ? 1.0 : theory[4]);
1538 
1539  pmc.diff.clear();
1540  for (int i=0;i<n_experiments;++i)
1541  {
1542  pmc.diff.push_back(pmc.value_exp(i,0)-pmc.value_th(i,0));
1543  }
1544 
1545  if (flav_debug) cout<<"Finished SL_measurements"<<endl;
1546 
1547  }
1548 
1549 
1551  void SL_likelihood(double &result)
1552  {
1553  using namespace Pipes::SL_likelihood;
1554 
1555  if (flav_debug) cout<<"Starting SL_likelihood"<<endl;
1556 
1557  predictions_measurements_covariances pmc = *Dep::SL_M;
1558 
1559  boost::numeric::ublas::matrix<double> cov=pmc.cov_exp;
1560 
1561  // adding theory and experimental covariance
1562  cov+=pmc.cov_th;
1563 
1564  //calculating a diff
1565  vector<double> diff;
1566  diff=pmc.diff;
1567 
1568  boost::numeric::ublas::matrix<double> cov_inv(pmc.dim, pmc.dim);
1569  InvertMatrix(cov, cov_inv);
1570 
1571  double Chi2=0;
1572  for (int i=0; i < pmc.dim; ++i)
1573  {
1574  for (int j=0; j<pmc.dim; ++j)
1575  {
1576  Chi2+= diff[i] * cov_inv(i,j)*diff[j];
1577  }
1578  }
1579 
1580  result=-0.5*Chi2;
1581 
1582  if (flav_debug) cout<<"Finished SL_likelihood"<<endl;
1583 
1584  if (flav_debug_LL) cout<<"Likelihood result SL_likelihood : "<< result<<endl;
1585 
1586  }
1587 
1588  // Helper function
1589  double G(const double x)
1590  {
1591  if(x)
1592  return (10.0 - 43.0*x + 78.0*pow(x,2) - 49.0*pow(x,3) + 4.0*pow(x,4) + 18.0*pow(x,3)*log(x)) / (3.0*pow(x - 1,4));
1593  else
1594  return 10.0/3;
1595  }
1596 
1597  // Contribution to mu -> e gamma from RHNs
1598  void RHN_muegamma(double &result)
1599  {
1600  using namespace Pipes::RHN_muegamma;
1601  SMInputs sminputs = *Dep::SMINPUTS;
1602 
1603  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1604  vector<double> ml = {sminputs.mE, sminputs.mMu, sminputs.mTau};
1605  vector<double> mnu = {real(m_nu(0,0)), real(m_nu(1,1)), real(m_nu(2,2)), *Param["M_1"], *Param["M_2"], *Param["M_3"]};
1606 
1607  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1608  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1609  Eigen::Matrix<complex<double>,3,6> U;
1610 
1611  for(int i=0; i<3; i++)
1612  for(int j=0; j<3; j++)
1613  {
1614  U(i,j) = Vnu(i,j);
1615  U(i,j+3) = Theta(i,j);
1616  }
1617 
1618  result = pow(sminputs.mMu,5)/(4 * sminputs.alphainv);
1619 
1620  // Form factors
1621  int e = 0, mu = 1;
1622  complex<double> k2l = FormFactors::K2L(mu, e, sminputs, U, ml, mnu);
1623  complex<double> k2r = FormFactors::K2R(mu, e, sminputs, U, ml, mnu);
1624 
1625  result *= (norm(k2l) + norm(k2r));
1626 
1627  result /= Dep::mu_minus_decay_rates->width_in_GeV;
1628 
1629  }
1630 
1631  // Contribution to tau -> e gamma from RHNs
1632  void RHN_tauegamma(double &result)
1633  {
1634  using namespace Pipes::RHN_tauegamma;
1635  SMInputs sminputs = *Dep::SMINPUTS;
1636 
1637  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1638  vector<double> ml = {sminputs.mE, sminputs.mMu, sminputs.mTau};
1639  vector<double> mnu = {real(m_nu(0,0)), real(m_nu(1,1)), real(m_nu(2,2)), *Param["M_1"], *Param["M_2"], *Param["M_3"]};
1640 
1641  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1642  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1643  Eigen::Matrix<complex<double>,3,6> U;
1644 
1645  for(int i=0; i<3; i++)
1646  for(int j=0; j<3; j++)
1647  {
1648  U(i,j) = Vnu(i,j);
1649  U(i,j+3) = Theta(i,j);
1650  }
1651 
1652  result = pow(sminputs.mTau,5)/(4*sminputs.alphainv);
1653 
1654  // Form factors
1655  int e = 0, tau = 2;
1656  complex<double> k2l = FormFactors::K2L(tau, e, sminputs, U, ml, mnu);
1657  complex<double> k2r = FormFactors::K2R(tau, e, sminputs, U, ml, mnu);
1658 
1659  result *= (norm(k2l) + norm(k2r));
1660 
1661  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1662 
1663  }
1664 
1665  // Contribution to tau -> mu gamma from RHNs
1666  void RHN_taumugamma(double &result)
1667  {
1668  using namespace Pipes::RHN_taumugamma;
1669  SMInputs sminputs = *Dep::SMINPUTS;
1670 
1671  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1672  vector<double> ml = {sminputs.mE, sminputs.mMu, sminputs.mTau};
1673  vector<double> mnu = {real(m_nu(0,0)), real(m_nu(1,1)), real(m_nu(2,2)), *Param["M_1"], *Param["M_2"], *Param["M_3"]};
1674 
1675  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1676  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1677  Eigen::Matrix<complex<double>,3,6> U;
1678 
1679  for(int i=0; i<3; i++)
1680  for(int j=0; j<3; j++)
1681  {
1682  U(i,j) = Vnu(i,j);
1683  U(i,j+3) = Theta(i,j);
1684  }
1685 
1686  result = pow(sminputs.mTau,5)/(4 * sminputs.alphainv);
1687 
1688  // Form factors
1689  int mu = 1, tau = 2;
1690  complex<double> k2l = FormFactors::K2L(tau, mu, sminputs, U, ml, mnu);
1691  complex<double> k2r = FormFactors::K2R(tau, mu, sminputs, U, ml, mnu);
1692 
1693  result *= (norm(k2l) + norm(k2r));
1694 
1695  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1696  }
1697 
1698  // General contribution to l_\alpha^- -> l_\beta^- l_\gamma^- l_\delta^+ from RHNs
1699  double RHN_l2lll(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix3cd Vnu, Eigen::Matrix3cd Theta, Eigen::Matrix3cd m_nu, double M1, double M2, double M3, double mH)
1700  {
1701  vector<double> ml = {sminputs.mE, sminputs.mMu, sminputs.mTau};
1702  vector<double> mnu = {real(m_nu(0,0)), real(m_nu(1,1)), real(m_nu(2,2)), M1, M2, M3};
1703 
1704  Eigen::Matrix<complex<double>,3,6> U;
1705 
1706  for(int i=0; i<3; i++)
1707  for(int j=0; j<3; j++)
1708  {
1709  U(i,j) = Vnu(i,j);
1710  U(i,j+3) = Theta(i,j);
1711  }
1712 
1713  // Form factors
1714  complex<double> k2l = FormFactors::K2L(alpha, beta, sminputs, U, ml, mnu);
1715  complex<double> k2r = FormFactors::K2R(alpha, beta, sminputs, U, ml, mnu);
1716  complex<double> k1r = FormFactors::K1R(alpha, beta, sminputs, U, mnu);
1717  complex<double> asll = FormFactors::ASLL(alpha, beta, gamma, delta, sminputs, U, ml, mnu, mH);
1718  complex<double> aslr = FormFactors::ASLR(alpha, beta, gamma, delta, sminputs, U, ml, mnu, mH);
1719  complex<double> asrl = FormFactors::ASRL(alpha, beta, gamma, delta, sminputs, U, ml, mnu, mH);
1720  complex<double> asrr = FormFactors::ASRR(alpha, beta, gamma, delta, sminputs, U, ml, mnu, mH);
1721  complex<double> avll = FormFactors::AVLL(alpha, beta, gamma, delta, sminputs, U, ml, mnu);
1722  complex<double> avlr = FormFactors::AVLR(alpha, beta, gamma, delta, sminputs, U, ml, mnu);
1723  complex<double> avrl = FormFactors::AVLL(alpha, beta, gamma, delta, sminputs, U, ml, mnu);
1724  complex<double> avrr = FormFactors::AVRR(alpha, beta, gamma, delta, sminputs, U, ml, mnu);
1725 
1726  complex<double> avhatll = avll;
1727  complex<double> avhatlr = avlr;
1728  complex<double> avhatrl = avrl + 4. * pi / sminputs.alphainv * k1r;
1729  complex<double> avhatrr = avrr + 4. * pi / sminputs.alphainv * k1r;
1730 
1731  double l2lll = 0;
1732  if(beta == gamma and gamma == delta) // l(alpha)- -> l(beta)- l(beta)- l(beta)+
1733  {
1734  l2lll = real(16. * pow(pi,2) / pow(sminputs.alphainv,2) * (norm(k2l) + norm(k2r)) * (16./3.*log(ml[alpha]/ml[beta]) - 22./3.) + 1./24. * (norm(asll) + norm(asrr) + 2.*norm(aslr) + 2.*norm(asrl)) + 1./3. * (2.*norm(avhatll) + 2.*norm(avhatrr) + norm(avhatlr) + norm(avhatrl)) + 4.*pi/(3.*sminputs.alphainv)*(k2l*conj(asrl - 2.*avhatrl - 4.*avhatrr) + conj(k2l)*(asrl - 2.*avhatrl - 4.*avhatrr) + k2r*conj(aslr - 2.*avhatlr - 4.*avhatll) + conj(k2r)*(aslr - 2.*avhatlr - 4.*avhatll)) - 1./6. * (aslr*conj(avhatlr) + asrl*conj(avhatrl) + conj(aslr)*avhatlr + conj(asrl)*avhatrl));
1735  }
1736  else if(gamma == delta) // l(alpha)- -> l(beta)- l(gamma)- l(gamma)+
1737  {
1738  l2lll = real(16. *pow(pi,2) / pow(sminputs.alphainv,2) * (norm(k2l) + norm(k2r)) * (16./3.*log(ml[alpha]/ml[gamma]) - 8.) + 1./12. *(norm(asll) + norm(asrr) + norm(aslr) + norm(asrl)) + 1./3. * (norm(avhatll) + norm(avhatrr) + norm(avhatlr) + norm(avhatrl)) + 8.*pi/(3.*sminputs.alphainv) * (k2l*conj(avhatrl + avhatrr) + k2r*conj(avhatlr + avhatll) + conj(k2l)*(avhatrl + avhatrr) + conj(k2r)*(avhatlr + avhatll)));
1739  }
1740  else if(beta == gamma) // l(alpha)- -> l(beta)- l(beta)- l(delta)+
1741  {
1742  l2lll = real(1./24. * (norm(asll) + norm(asrr) + 2.*norm(aslr) + 2.*norm(asrl)) + 1./3.*(2.*norm(avhatll) + 2.*norm(avhatrr) + norm(avhatlr) + norm(avhatrl)) - 1./6.*(aslr*conj(avhatlr) + asrl*conj(avhatrl) + conj(aslr)*avhatlr + conj(asrl)*avhatrl));
1743  }
1744  return l2lll;
1745 
1746  }
1747 
1748  // Contribution to mu -> e e e from RHNs
1749  void RHN_mueee(double &result)
1750  {
1751  using namespace Pipes::RHN_mueee;
1752  SMInputs sminputs = *Dep::SMINPUTS;
1753 
1754  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1755  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1756  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1757 
1758  result = pow(sminputs.mMu,5)/(512*pow(pi,3));
1759 
1760  int e = 0, mu = 1;
1761  result *= RHN_l2lll(mu, e, e, e, sminputs, Vnu, Theta, m_nu, *Param["M_1"], *Param["M_2"], *Param["M_3"], *Param["mH"]);
1762 
1763  result /= Dep::mu_minus_decay_rates->width_in_GeV;
1764 
1765  }
1766 
1767  // Contribution to tau -> e e e from RHNs
1768  void RHN_taueee(double &result)
1769  {
1770  using namespace Pipes::RHN_taueee;
1771  SMInputs sminputs = *Dep::SMINPUTS;
1772 
1773  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1774  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1775  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1776 
1777  result = pow(sminputs.mTau,5)/(512*pow(pi,3));
1778 
1779  int e = 0, tau = 2;
1780  result *= RHN_l2lll(tau, e, e, e, sminputs, Vnu, Theta, m_nu, *Param["M_1"], *Param["M_2"], *Param["M_3"], *Param["mH"]);
1781 
1782  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1783 
1784  }
1785 
1786  // Contribution to tau -> mu mu mu from RHNs
1787  void RHN_taumumumu(double &result)
1788  {
1789  using namespace Pipes::RHN_taumumumu;
1790  SMInputs sminputs = *Dep::SMINPUTS;
1791 
1792  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1793  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1794  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1795 
1796  result = pow(sminputs.mTau,5)/(512*pow(pi,3));
1797 
1798  int mu = 1, tau = 2;
1799  result *= RHN_l2lll(tau, mu, mu, mu, sminputs, Vnu, Theta, m_nu, *Param["M_1"], *Param["M_2"], *Param["M_3"], *Param["mH"]);
1800 
1801  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1802 
1803  }
1804 
1805  // Contribution to tau^- -> mu^- e^- e^+ from RHNs
1806  void RHN_taumuee(double &result)
1807  {
1808  using namespace Pipes::RHN_taumuee;
1809  SMInputs sminputs = *Dep::SMINPUTS;
1810 
1811  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1812  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1813  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1814 
1815  result = pow(sminputs.mTau,5)/(512*pow(pi,3));
1816 
1817  int e = 0, mu = 1, tau = 2;
1818  result *= RHN_l2lll(tau, mu, e, e, sminputs, Vnu, Theta, m_nu, *Param["M_1"], *Param["M_2"], *Param["M_3"], *Param["mH"]);
1819 
1820  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1821  }
1822 
1823  // Contribution to tau^- -> e^- e^- mu^+ from RHNs
1824  void RHN_taueemu(double &result)
1825  {
1826  using namespace Pipes::RHN_taueemu;
1827  SMInputs sminputs = *Dep::SMINPUTS;
1828 
1829  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1830  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1831  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1832 
1833  result = pow(sminputs.mTau,5)/(512*pow(pi,3));
1834 
1835  int e = 0, mu = 1, tau = 2;
1836  result *= RHN_l2lll(tau, e, e, mu, sminputs, Vnu, Theta, m_nu, *Param["M_1"], *Param["M_2"], *Param["M_3"], *Param["mH"]);
1837 
1838  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1839  }
1840 
1841  // Contribution to tau^- -> e^- mu^- mu^+ from RHNs
1842  void RHN_tauemumu(double &result)
1843  {
1844  using namespace Pipes::RHN_tauemumu;
1845  SMInputs sminputs = *Dep::SMINPUTS;
1846 
1847  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1848  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1849  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1850 
1851  result = pow(sminputs.mTau,5)/(512*pow(pi,3));
1852 
1853  int e = 0, mu = 1, tau = 2;
1854  result *= RHN_l2lll(tau, e, mu, mu, sminputs, Vnu, Theta, m_nu, *Param["M_1"], *Param["M_2"], *Param["M_3"], *Param["mH"]);
1855 
1856  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1857  }
1858 
1859  // Contribution to tau^- -> mu^- mu^- e^+ from RHNs
1860  void RHN_taumumue(double &result)
1861  {
1862  using namespace Pipes::RHN_taumumue;
1863  SMInputs sminputs = *Dep::SMINPUTS;
1864 
1865  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1866  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1867  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1868 
1869  result = pow(sminputs.mTau,5)/(512*pow(pi,3));
1870 
1871  int e = 0, mu = 1, tau = 2;
1872  result *= RHN_l2lll(tau, mu, mu, e, sminputs, Vnu, Theta, m_nu, *Param["M_1"], *Param["M_2"], *Param["M_3"], *Param["mH"]);
1873 
1874  result /= Dep::tau_minus_decay_rates->width_in_GeV;
1875  }
1876 
1877  // Form factors for to mu - e conversion
1878  void RHN_mue_FF(const SMInputs sminputs, std::vector<double> &mnu, Eigen::Matrix<complex<double>,3,6> &U, const double mH, complex<double> &g0SL, complex<double> &g0SR, complex<double> &g0VL, complex<double> &g0VR, complex<double> &g1SL, complex<double> &g1SR, complex<double> &g1VL, complex<double> &g1VR)
1879  {
1880  vector<double> ml = {sminputs.mE, sminputs.mMu, sminputs.mTau};
1881 
1882  int e = 0, mu = 1;
1883  complex<double> k1r = FormFactors::K1R(mu, e, sminputs, U, mnu);
1884  complex<double> k2l = FormFactors::K2L(mu, e, sminputs, U, ml, mnu);
1885  complex<double> k2r = FormFactors::K2R(mu, e, sminputs, U, ml, mnu);
1886 
1887  int u = 0, d =0, s = 1;
1888  complex<double> CVLLu = FormFactors::CVLL(mu, e, u, u, sminputs, U, ml, mnu);
1889  complex<double> CVLLd = FormFactors::BVLL(mu, e, d, d, sminputs, U, ml, mnu);
1890  complex<double> CVLLs = FormFactors::BVLL(mu, e, s, s, sminputs, U, ml, mnu);
1891  complex<double> CVLRu = FormFactors::CVLR(mu, e, u, u, sminputs, U, ml, mnu);
1892  complex<double> CVLRd = FormFactors::BVLR(mu, e, d, d, sminputs, U, ml, mnu);
1893  complex<double> CVLRs = FormFactors::BVLR(mu, e, s, s, sminputs, U, ml, mnu);
1894  complex<double> CVRLu = FormFactors::CVRL(mu, e, u, u, sminputs, U, ml, mnu);
1895  complex<double> CVRLd = FormFactors::BVRL(mu, e, d, d, sminputs, U, ml, mnu);
1896  complex<double> CVRLs = FormFactors::BVRL(mu, e, s, s, sminputs, U, ml, mnu);
1897  complex<double> CVRRu = FormFactors::CVRR(mu, e, u, u, sminputs, U, ml, mnu);
1898  complex<double> CVRRd = FormFactors::BVRR(mu, e, d, d, sminputs, U, ml, mnu);
1899  complex<double> CVRRs = FormFactors::BVRR(mu, e, s, s, sminputs, U, ml, mnu);
1900 
1901  complex<double> CSLLu = FormFactors::CSLL(mu, e, u, u, sminputs, U, ml, mnu, mH);
1902  complex<double> CSLLd = FormFactors::BSLL(mu, e, d, d, sminputs, U, ml, mnu, mH);
1903  complex<double> CSLLs = FormFactors::BSLL(mu, e, s, s, sminputs, U, ml, mnu, mH);
1904  complex<double> CSLRu = FormFactors::CSLL(mu, e, u, u, sminputs, U, ml, mnu, mH);
1905  complex<double> CSLRd = FormFactors::BSLL(mu, e, d, d, sminputs, U, ml, mnu, mH);
1906  complex<double> CSLRs = FormFactors::BSLL(mu, e, s, s, sminputs, U, ml, mnu, mH);
1907  complex<double> CSRLu = FormFactors::CSLL(mu, e, u, u, sminputs, U, ml ,mnu, mH);
1908  complex<double> CSRLd = FormFactors::BSLL(mu, e, d, d, sminputs, U, ml ,mnu, mH);
1909  complex<double> CSRLs = FormFactors::BSLL(mu, e, s, s, sminputs, U, ml ,mnu, mH);
1910  complex<double> CSRRu = FormFactors::CSLL(mu, e, u, u, sminputs, U, ml ,mnu, mH);
1911  complex<double> CSRRd = FormFactors::BSLL(mu, e, d, d, sminputs, U, ml, mnu, mH);
1912  complex<double> CSRRs = FormFactors::BSLL(mu, e, s, s, sminputs, U, ml ,mnu, mH);
1913 
1914  double Qu = 2./3.;
1915  complex<double> gVLu = sqrt(2)/sminputs.GF * (4.*pi / sminputs.alphainv * Qu * (0. - k2r) - 0.5*(CVLLu + CVLRu));
1916  complex<double> gSLu = -1./(sqrt(2)*sminputs.GF)*(CSLLu + CSLRu);
1917  complex<double> gVRu = sqrt(2)/sminputs.GF * (4.*pi / sminputs.alphainv * Qu * (k1r - k2l) - 0.5*(CVRRu + CVRLu));
1918  complex<double> gSRu = -1./(sqrt(2)*sminputs.GF)*(CSRRu + CSRLu);
1919 
1920  double Qd = -1./3.;
1921  complex<double> gVLd = sqrt(2)/sminputs.GF * (4.*pi / sminputs.alphainv * Qd * (0. - k2r) - 0.5*(CVLLd + CVLRd));
1922  complex<double> gSLd = -1./(sqrt(2)*sminputs.GF)*(CSLLd + CSLRd);
1923  complex<double> gVRd = sqrt(2)/sminputs.GF * (4.*pi / sminputs.alphainv * Qd * (k1r - k2l) - 0.5*(CVRRd + CVRLd));
1924  complex<double> gSRd = -1./(sqrt(2)*sminputs.GF)*(CSRRd + CSRLd);
1925 
1926  double Qs = -1./3.;
1927  complex<double> gVLs = sqrt(2)/sminputs.GF * (4.*pi / sminputs.alphainv * Qs * (0. - k2r) - 0.5*(CVLLs + CVLRs));
1928  complex<double> gSLs = -1./(sqrt(2)*sminputs.GF)*(CSLLs + CSLRs);
1929  complex<double> gVRs = sqrt(2)/sminputs.GF * (4.*pi / sminputs.alphainv * Qs * (k1r - k2l) - 0.5*(CVRRs + CVRLs));
1930  complex<double> gSRs = -1./(sqrt(2)*sminputs.GF)*(CSRRs + CSRLs);
1931 
1932  double GVup = 2, GVdn = 2, GVdp = 1, GVun = 1, GVsp = 0, GVsn = 0;
1933  double GSup = 5.1, GSdn = 5.1, GSdp = 4.3, GSun = 4.3, GSsp = 2.5, GSsn = 2.5;
1934 
1935  g0SL = 0.5*(gSLu*(GSup + GSun) + gSLd*(GSdp + GSdn) + gSLs*(GSsp + GSsn));
1936  g0SR = 0.5*(gSRu*(GSup + GSun) + gSRd*(GSdp + GSdn) + gSRs*(GSsp + GSsn));
1937  g0VL = 0.5*(gVLu*(GVup + GVun) + gVLd*(GVdp + GVdn) + gVLs*(GVsp + GVsn));
1938  g0VR = 0.5*(gVRu*(GVup + GVun) + gVRd*(GVdp + GVdn) + gVRs*(GVsp + GVsn));
1939  g1SL = 0.5*(gSLu*(GSup - GSun) + gSLd*(GSdp - GSdn) + gSLs*(GSsp - GSsn));
1940  g1SR = 0.5*(gSRu*(GSup - GSun) + gSRd*(GSdp - GSdn) + gSRs*(GSsp - GSsn));
1941  g1VL = 0.5*(gVLu*(GVup - GVun) + gVLd*(GVdp - GVdn) + gVLs*(GVsp - GVsn));
1942  g1VR = 0.5*(gVRu*(GVup - GVun) + gVRd*(GVdp - GVdn) + gVRs*(GVsp - GVsn));
1943 
1944  }
1945 
1946  // Contribution to mu - e conversion in Ti nuclei from RHNs
1947  void RHN_mueTi(double &result)
1948  {
1949  using namespace Pipes::RHN_mueTi;
1950  const SMInputs sminputs = *Dep::SMINPUTS;
1951  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1952  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1953  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1954 
1955  vector<double> mnu = {real(m_nu(0,0)), real(m_nu(1,1)), real(m_nu(2,2)), *Param["M_1"], *Param["M_2"], *Param["M_3"]};
1956  Eigen::Matrix<complex<double>,3,6> U;
1957 
1958  for(int i=0; i<3; i++)
1959  for(int j=0; j<3; j++)
1960  {
1961  U(i,j) = Vnu(i,j);
1962  U(i,j+3) = Theta(i,j);
1963  }
1964 
1965  complex<double> g0SL, g0SR, g0VL, g0VR, g1SL, g1SR, g1VL, g1VR;
1966  RHN_mue_FF(sminputs, mnu, U, *Param["mH"], g0SL, g0SR, g0VL, g0VR, g1SL, g1SR, g1VL, g1VR);
1967 
1968  // Parameters for Ti, from Table 1 in 1209.2679 for Ti
1969  double Z = 22, N = 26;
1970  double Zeff = 17.6, Fp = 0.54;
1971  double hbar = 6.582119514e-25; // GeV * s
1972  double GammaCapt = 2.59e6 * hbar;
1973 
1974  result = (pow(sminputs.GF,2)*pow(sminputs.mMu,5)*pow(Zeff,4)*pow(Fp,2)) / (8.*pow(pi,4)*pow(sminputs.alphainv,3)*Z*GammaCapt) * (norm((Z+N)*(g0VL + g0SL) + (Z-N)*(g1VL + g1SL)) + norm((Z+N)*(g0VR + g0SR) + (Z-N)*(g1VR + g1SR)));
1975 
1976  }
1977 
1978  // Contribution to mu - e conversion in Au nuclei from RHNs
1979  void RHN_mueAu(double &result)
1980  {
1981  using namespace Pipes::RHN_mueAu;
1982  const SMInputs sminputs = *Dep::SMINPUTS;
1983  Eigen::Matrix3cd m_nu = *Dep::m_nu;
1984  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
1985  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
1986 
1987  vector<double> mnu = {real(m_nu(0,0)), real(m_nu(1,1)), real(m_nu(2,2)), *Param["M_1"], *Param["M_2"], *Param["M_3"]};
1988  Eigen::Matrix<complex<double>,3,6> U;
1989 
1990  for(int i=0; i<3; i++)
1991  for(int j=0; j<3; j++)
1992  {
1993  U(i,j) = Vnu(i,j);
1994  U(i,j+3) = Theta(i,j);
1995  }
1996 
1997  complex<double> g0SL, g0SR, g0VL, g0VR, g1SL, g1SR, g1VL, g1VR;
1998  RHN_mue_FF(sminputs, mnu, U, *Param["mH"], g0SL, g0SR, g0VL, g0VR, g1SL, g1SR, g1VL, g1VR);
1999 
2000 
2001  // Parameters for Au, from Table 1 in 1209.2679 for Au
2002  double Z = 79, N = 118;
2003  double Zeff = 33.5, Fp = 0.16;
2004  double hbar = 6.582119514e-25; // GeV * s
2005  double GammaCapt = 13.07e6 * hbar;
2006 
2007  result = (pow(sminputs.GF,2)*pow(sminputs.mMu,5)*pow(Zeff,4)*pow(Fp,2)) / (8.*pow(pi,4)*pow(sminputs.alphainv,3)*Z*GammaCapt) * (norm((Z+N)*(g0VL + g0SL) + (Z-N)*(g1VL + g1SL)) + norm((Z+N)*(g0VR + g0SR) + (Z-N)*(g1VR + g1SR)));
2008 
2009  }
2010 
2011 
2012  // Contribution to mu - e conversion in Pb nuclei from RHNs
2013  void RHN_muePb(double &result)
2014  {
2015  using namespace Pipes::RHN_muePb;
2016  const SMInputs sminputs = *Dep::SMINPUTS;
2017  Eigen::Matrix3cd m_nu = *Dep::m_nu;
2018  Eigen::Matrix3cd Vnu = *Dep::SeesawI_Vnu;
2019  Eigen::Matrix3cd Theta = *Dep::SeesawI_Theta;
2020 
2021  vector<double> mnu = {real(m_nu(0,0)), real(m_nu(1,1)), real(m_nu(2,2)), *Param["M_1"], *Param["M_2"], *Param["M_3"]};
2022  Eigen::Matrix<complex<double>,3,6> U;
2023 
2024  for(int i=0; i<3; i++)
2025  for(int j=0; j<3; j++)
2026  {
2027  U(i,j) = Vnu(i,j);
2028  U(i,j+3) = Theta(i,j);
2029  }
2030 
2031  complex<double> g0SL, g0SR, g0VL, g0VR, g1SL, g1SR, g1VL, g1VR;
2032  RHN_mue_FF(sminputs, mnu, U, *Param["mH"], g0SL, g0SR, g0VL, g0VR, g1SL, g1SR, g1VL, g1VR);
2033 
2034 
2035  // Parameters for Pb, from Table 1 in 1209.2679 for Pb
2036  double Z = 82, N = 126;
2037  double Zeff = 34., Fp = 0.15;
2038  double hbar = 6.582119514e-25; // GeV * s
2039  double GammaCapt = 13.45e6 * hbar;
2040 
2041  result = (pow(sminputs.GF,2)*pow(sminputs.mMu,5)*pow(Zeff,4)*pow(Fp,2)) / (8.*pow(pi,4)*pow(sminputs.alphainv,3)*Z*GammaCapt) * (norm((Z+N)*(g0VL + g0SL) + (Z-N)*(g1VL + g1SL)) + norm((Z+N)*(g0VR + g0SR) + (Z-N)*(g1VR + g1SR)));
2042 
2043  }
2044 
2046  void l2lgamma_likelihood(double &result)
2047  {
2048  using namespace Pipes::l2lgamma_likelihood;
2049 
2050  static bool first = true;
2051  static boost::numeric::ublas::matrix<double> cov_exp, value_exp;
2052  static double th_err[3];
2053  double theory[3];
2054 
2055  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
2056  if (first)
2057  {
2058  // Read in experimental measuremens
2059  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
2060  fread.debug_mode(flav_debug);
2061 
2062  // mu -> e gamma
2063  fread.read_yaml_measurement("flav_data.yaml", "BR_muegamma");
2064  // tau -> e gamma
2065  fread.read_yaml_measurement("flav_data.yaml", "BR_tauegamma");
2066  // tau -> mu gamma
2067  fread.read_yaml_measurement("flav_data.yaml", "BR_taumugamma");
2068 
2069  fread.initialise_matrices();
2070  cov_exp=fread.get_exp_cov();
2071  value_exp=fread.get_exp_value();
2072 
2073  for (int i = 0; i < 3; ++i)
2074  th_err[i] = fread.get_th_err()(i,0).first;
2075 
2076  // Init over.
2077  first = false;
2078  }
2079 
2080  theory[0] = *Dep::muegamma;
2081  if(flav_debug) cout << "mu- -> e- gamma = " << theory[0] << endl;
2082  theory[1] = *Dep::tauegamma;
2083  if(flav_debug) cout << "tau- -> e- gamma = " << theory[1] << endl;
2084  theory[2] = *Dep::taumugamma;
2085  if(flav_debug) cout << "tau- -> mu- gamma = " << theory[2] << endl;
2086 
2087  result = 0;
2088  for (int i = 0; i < 3; ++i)
2089  result += Stats::gaussian_upper_limit(theory[i], value_exp(i,0), th_err[i], sqrt(cov_exp(i,i)), false);
2090 
2091  }
2092 
2094  void l2lll_likelihood(double &result)
2095  {
2096  using namespace Pipes::l2lll_likelihood;
2097 
2098  static bool first = true;
2099  static boost::numeric::ublas::matrix<double> cov_exp, value_exp;
2100  static double th_err[7];
2101  double theory[7];
2102 
2103 
2104  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
2105  if (first)
2106  {
2107  // Read in experimental measuremens
2108  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
2109  fread.debug_mode(flav_debug);
2110 
2111  // mu- -> e- e- e+
2112  fread.read_yaml_measurement("flav_data.yaml", "BR_mueee");
2113  // tau- -> e- e- e+
2114  fread.read_yaml_measurement("flav_data.yaml", "BR_taueee");
2115  // tau- -> mu- mu- mu+
2116  fread.read_yaml_measurement("flav_data.yaml", "BR_taumumumu");
2117  // tau- -> mu- e- e+
2118  fread.read_yaml_measurement("flav_data.yaml", "BR_taumuee");
2119  // tau- -> e- e- mu+
2120  fread.read_yaml_measurement("flav_data.yaml", "BR_taueemu");
2121  // tau- -> e- mu- mu+
2122  fread.read_yaml_measurement("flav_data.yaml", "BR_tauemumu");
2123  // tau- -> mu- mu- e+
2124  fread.read_yaml_measurement("flav_data.yaml", "BR_taumumue");
2125 
2126  fread.initialise_matrices();
2127  cov_exp=fread.get_exp_cov();
2128  value_exp=fread.get_exp_value();
2129 
2130  for (int i = 0; i < 7; ++i)
2131  th_err[i] = fread.get_th_err()(i,0).first;
2132 
2133  // Init over.
2134  first = false;
2135  }
2136 
2137  theory[0] = *Dep::mueee;
2138  if(flav_debug) cout << "mu- -> e- e- e+ = " << theory[0] << endl;
2139  theory[1] = *Dep::taueee;
2140  if(flav_debug) cout << "tau- -> e- e- e+ = " << theory[1] << endl;
2141  theory[2] = *Dep::taumumumu;
2142  if(flav_debug) cout << "tau- -> mu- mu- mu+ = " << theory[2] << endl;
2143  theory[3] = *Dep::taumuee;
2144  if(flav_debug) cout << "tau- -> mu- e- e- = " << theory[3] << endl;
2145  theory[4] = *Dep::taueemu;
2146  if(flav_debug) cout << "tau- -> e- e- mu+ = " << theory[4] << endl;
2147  theory[5] = *Dep::tauemumu;
2148  if(flav_debug) cout << "tau- -> e- mu- mu+ = " << theory[5] << endl;
2149  theory[6] = *Dep::taumumue;
2150  if(flav_debug) cout << "tau- -> mu- mu- e+ = " << theory[6] << endl;
2151 
2152  result = 0;
2153  for (int i = 0; i < 7; ++i)
2154  result += Stats::gaussian_upper_limit(theory[i], value_exp(i,0), th_err[i], sqrt(cov_exp(i,i)), false);
2155 
2156  }
2157 
2159  void mu2e_likelihood(double &result)
2160  {
2161  using namespace Pipes::mu2e_likelihood;
2162 
2163  static bool first = true;
2164  static boost::numeric::ublas::matrix<double> cov_exp, value_exp;
2165  static int n_measurements = 3;
2166  static double th_err[3];
2167  double theory[3];
2168 
2169 
2170  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
2171  if (first)
2172  {
2173  // Read in experimental measuremens
2174  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
2175  fread.debug_mode(flav_debug);
2176 
2177  // mu - e (Ti)
2178  fread.read_yaml_measurement("flav_data.yaml", "R_mueTi");
2179  // mu - e (Au)
2180  fread.read_yaml_measurement("flav_data.yaml", "R_mueAu");
2181  // mu - e (Pb)
2182  fread.read_yaml_measurement("flav_data.yaml", "R_muePb");
2183 
2184  fread.initialise_matrices();
2185  cov_exp=fread.get_exp_cov();
2186  value_exp=fread.get_exp_value();
2187 
2188  for (int i = 0; i < n_measurements; ++i)
2189  th_err[i] = fread.get_th_err()(i,0).first;
2190 
2191  // Init over.
2192  first = false;
2193  }
2194 
2195  theory[0] = *Dep::mueTi;
2196  if(flav_debug) cout << "mu - e (Ti) = " << theory[0] << endl;
2197  theory[1] = *Dep::mueAu;
2198  if(flav_debug) cout << "mu - e (Au) = " << theory[1] << endl;
2199  theory[2] = *Dep::muePb;
2200  if(flav_debug) cout << "mu - e (Pb) = " << theory[2] << endl;
2201 
2202  result = 0;
2203  for (int i = 0; i < n_measurements; ++i)
2204  result += Stats::gaussian_upper_limit(theory[i], value_exp(i,0), th_err[i], sqrt(cov_exp(i,i)), false);
2205 
2206  }
2207 
2210  {
2211  using namespace Pipes::LUV_measurements;
2212  static bool first = true;
2213 
2214  static double theory_RKstar_0045_11_err, theory_RKstar_11_60_err, theory_RK_err;
2215  if (flav_debug) cout<<"Starting LUV_measurements"<<endl;
2216 
2217  // Read and calculate things based on the observed data only the first time through, as none of it depends on the model parameters.
2218  if (first)
2219  {
2220  pmc.LL_name="LUV_likelihood";
2221 
2222  Flav_reader fread(GAMBIT_DIR "/FlavBit/data");
2223  fread.debug_mode(flav_debug);
2224 
2225  if (flav_debug) cout<<"Initiated Flav reader in LUV_measurements"<<endl;
2226  fread.read_yaml_measurement("flav_data.yaml", "RKstar_0045_11");
2227  fread.read_yaml_measurement("flav_data.yaml", "RKstar_11_60");
2228  fread.read_yaml_measurement("flav_data.yaml", "RK");
2229 
2230  if (flav_debug) cout<<"Finished reading LUV data"<<endl;
2231 
2232  fread.initialise_matrices();
2233 
2234  theory_RKstar_0045_11_err = fread.get_th_err()(0,0).first;
2235  theory_RKstar_11_60_err = fread.get_th_err()(1,0).first;
2236  theory_RK_err = fread.get_th_err()(2,0).first;
2237 
2238  pmc.value_exp=fread.get_exp_value();
2239  pmc.cov_exp=fread.get_exp_cov();
2240 
2241  pmc.value_th.resize(3,1);
2242  pmc.cov_th.resize(3,3);
2243 
2244  pmc.dim=3;
2245 
2246  // Init over and out.
2247  first = false;
2248  }
2249 
2250  // Get theory prediction
2251  pmc.value_th(0,0)=*Dep::RKstar_0045_11;
2252  pmc.value_th(1,0)=*Dep::RKstar_11_60;
2253  pmc.value_th(2,0)=*Dep::RK;
2254 
2255  // Compute error on theory prediction and populate the covariance matrix
2256  pmc.cov_th(0,0)=theory_RKstar_0045_11_err;
2257  pmc.cov_th(0,1)=0.;
2258  pmc.cov_th(0,2)=0.;
2259  pmc.cov_th(1,0)=0.;
2260  pmc.cov_th(1,1)=theory_RKstar_11_60_err;
2261  pmc.cov_th(1,2)=0.;
2262  pmc.cov_th(2,0)=0.;
2263  pmc.cov_th(2,1)=0.;
2264  pmc.cov_th(2,2)=theory_RK_err;
2265 
2266 
2267 
2268  // Save the differences between theory and experiment
2269  pmc.diff.clear();
2270  for (int i=0;i<3;++i)
2271  {
2272  pmc.diff.push_back(pmc.value_exp(i,0)-pmc.value_th(i,0));
2273  }
2274 
2275  if (flav_debug) cout<<"Finished LUV_measurements"<<endl;
2276 
2277 
2278  }
2280  void LUV_likelihood(double &result)
2281  {
2282  using namespace Pipes::LUV_likelihood;
2283 
2284  if (flav_debug) cout<<"Starting LUV_likelihood"<<endl;
2285 
2286  predictions_measurements_covariances pmc = *Dep::LUV_M;
2287 
2288  boost::numeric::ublas::matrix<double> cov=pmc.cov_exp;
2289 
2290  // adding theory and experimental covariance
2291  cov+=pmc.cov_th;
2292 
2293  //calculating a diff
2294  vector<double> diff;
2295  diff=pmc.diff;
2296 
2297  boost::numeric::ublas::matrix<double> cov_inv(pmc.dim, pmc.dim);
2298  InvertMatrix(cov, cov_inv);
2299 
2300  double Chi2=0;
2301  for (int i=0; i < pmc.dim; ++i)
2302  {
2303  for (int j=0; j<pmc.dim; ++j)
2304  {
2305  Chi2+= diff[i] * cov_inv(i,j)*diff[j];
2306  }
2307  }
2308 
2309  result=-0.5*Chi2;
2310 
2311  if (flav_debug) cout<<"Finished LUV_likelihood"<<endl;
2312 
2313  if (flav_debug_LL) cout<<"Likelihood result LUV_likelihood : "<< result<<endl;
2314 
2315  }
2316 
2317 
2318 
2319  }
2320 }
void SI_Bsee_untag(double &result)
Br Bs->ee decays for the untagged case (CP-averaged)
Definition: FlavBit.cpp:532
void SI_A_BXsmumu_zero(double &result)
Zero crossing of the forward-backward asymmetry of B -> X_s mu mu.
Definition: FlavBit.cpp:832
DEFINE_BKSTARMUMU(1.1, 2.5, 11, 25) DEFINE_BKSTARMUMU(2.5
void deltaMB_likelihood(double &result)
Likelihood for Delta Ms.
Definition: FlavBit.cpp:1272
complex< double > BVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void RHN_taumugamma(double &result)
Definition: FlavBit.cpp:1666
void RHN_taumumumu(double &result)
Definition: FlavBit.cpp:1787
This class is used to deliver both information defined in the Standard Model (or potentially just QED...
Structure for holding predicted and observed values of multiple observables, and experimental and the...
boost::numeric::ublas::matrix< double > value_th
void SI_BRBXsmumu_highq2(double &result)
Inclusive branching fraction B -> X_s mu mu at high q^2.
Definition: FlavBit.cpp:790
Holder for theory errors for B->K* mumu observables.
void RHN_mueAu(double &result)
Definition: FlavBit.cpp:1979
void b2ll_likelihood(double &result)
Likelihood for rare purely leptonic B decays.
Definition: FlavBit.cpp:1414
boost::numeric::ublas::matrix< double > get_exp_cov()
Return the covariance matrix.
Definition: Flav_reader.hpp:90
boost::numeric::ublas::matrix< double > get_th_cov(std::vector< str > observables)
Return theory error covariance matrix for selected observables.
START_MODEL beta
Definition: Axions.hpp:35
complex< double > AVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void SI_BDstarmunu(double &result)
Br B -> D* mu nu.
Definition: FlavBit.cpp:684
boost::numeric::ublas::matrix< double > value_exp
Type definition header for module FlavBit.
boost::numeric::ublas::matrix< double > cov_th
void SI_RDstar(double &result)
B->D* tau nu / B-> D* e nu decays.
Definition: FlavBit.cpp:720
Funk delta(std::string arg, double pos, double width)
Definition: daFunk.hpp:902
#define LOCAL_INFO
Definition: local_info.hpp:34
const double mt
Definition: topness.h:39
void SI_RKstar_11_60(double &result)
RK* in intermediate q^2.
Definition: FlavBit.cpp:907
STL namespace.
void SI_Bsmumu_untag(double &result)
Br Bs->mumu decays for the untagged case (CP-averaged)
Definition: FlavBit.cpp:517
double gaussian_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is Gaussian...
Definition: statistics.cpp:133
void SI_Rmu23(double &result)
2-to-3-body decay ratio for semileptonic K and pi decays
Definition: FlavBit.cpp:748
complex< double > AVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void RHN_tauegamma(double &result)
Definition: FlavBit.cpp:1632
void SI_BDtaunu(double &result)
Br B -> D tau nu.
Definition: FlavBit.cpp:618
void b2sgamma_likelihood(double &result)
Likelihood for b->s gamma.
Definition: FlavBit.cpp:1309
void RHN_taumuee(double &result)
Definition: FlavBit.cpp:1806
void SI_RD(double &result)
B-> D tau nu / B-> D e nu decays.
Definition: FlavBit.cpp:706
void LUV_likelihood(double &result)
Likelihood for LUV in b->sll.
Definition: FlavBit.cpp:2280
void RHN_taumumue(double &result)
Definition: FlavBit.cpp:1860
void FH_bsgamma(double &result)
These functions extract observables from a FeynHiggs flavour result.
Definition: FlavBit.cpp:1102
const double mW
Definition: topness.h:40
void RHN_RKstar_0045_11(double &result)
Definition: FlavBit.cpp:920
void SLHAea_add(SLHAstruct &slha, const str &block, const int index, const double value, const str &comment="", const bool overwrite=false)
Add an entry to an SLHAea object (if overwrite=false, only if it doesn&#39;t already exist) ...
Helper utilities for FlavBit.
void SI_Dmunu(double &result)
Br D -> mu nu.
Definition: FlavBit.cpp:604
void SI_fill(parameters &result)
Fill SuperIso model info structure.
Definition: FlavBit.cpp:86
complex< double > K1R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > mnu)
double RHN_l2lll(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix3cd Vnu, Eigen::Matrix3cd Theta, Eigen::Matrix3cd m_nu, double M1, double M2, double M3, double mH)
Definition: FlavBit.cpp:1699
Declarations of B->Kstar mumu theory error class.
complex< double > ASRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
complex< double > CSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
const double hbar
void SI_AI_BKstarmumu_zero(double &result)
Zero crossing of isospin asymmetry of B-> K* mu mu.
Definition: FlavBit.cpp:1054
void RHN_muegamma(double &result)
Definition: FlavBit.cpp:1598
void FH_DeltaMs(double &result)
Definition: FlavBit.cpp:1110
void b2ll_measurements(predictions_measurements_covariances &pmc)
Measurements for rare purely leptonic B decays.
Definition: FlavBit.cpp:1348
START_MODEL alpha
void SI_A_BXsmumu_highq2(double &result)
Forward-backward asymmetry of B -> X_s mu mu at high q^2.
Definition: FlavBit.cpp:818
Spectrum W_plus_decay_rates
void SI_AI_BKstarmumu(double &result)
Isospin asymmetry of B-> K* mu mu.
Definition: FlavBit.cpp:1040
complex< double > CVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
SLHAea::Coll SLHAstruct
Less confusing name for SLHAea container class.
void SL_likelihood(double &result)
Likelihood for tree-level leptonic and semileptonic B decays.
Definition: FlavBit.cpp:1551
void mu2e_likelihood(double &result)
Likelihood for mu - e conversion in nuclei.
Definition: FlavBit.cpp:2159
GAMBIT error class.
Definition: exceptions.hpp:136
Declarations of statistical utilities.
void SI_bsgamma(double &result)
Br b-> s gamma decays.
Definition: FlavBit.cpp:502
void b2sll_measurements(predictions_measurements_covariances &pmc)
Measurements for electroweak penguin decays.
Definition: FlavBit.cpp:1118
const double pi
void b2sll_likelihood(double &result)
Likelihood for electroweak penguin decays.
Definition: FlavBit.cpp:1232
double gaussian_upper_limit(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a gaussian log-likelihood for an upper limit.
Definition: statistics.cpp:184
const double mu
Definition: SM_Z.hpp:42
Rollcall header for module FlavBit.
void RHN_taueemu(double &result)
Definition: FlavBit.cpp:1824
double G(const double x)
Definition: FlavBit.cpp:1589
Loop functions for flavour violating decays of charged leptons (from hep-ph/9403398) And for RK from ...
complex< double > ASLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
void initialise_matrices()
Compute the covariance matrix and populate the measurement and theory error vectors.
complex< double > ASRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
void RHN_tauemumu(double &result)
Definition: FlavBit.cpp:1842
complex< double > K2L(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
Header file that includes all GAMBIT headers required for a module source file.
void LUV_measurements(predictions_measurements_covariances &pmc)
Measurements for LUV in b->sll.
Definition: FlavBit.cpp:2209
complex< double > AVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void SI_BRBXstautau_highq2(double &result)
Inclusive branching fraction B -> X_s tau tau at high q^2.
Definition: FlavBit.cpp:846
void SI_delta0(double &result)
Delta_0 (CP-averaged isospin asymmetry of B -> K* gamma)
Definition: FlavBit.cpp:762
const bool flav_debug
Definition: FlavBit.cpp:71
void SI_A_BXsmumu_lowq2(double &result)
Forward-backward asymmetry of B -> X_s mu mu at low q^2.
Definition: FlavBit.cpp:804
void SI_Dstaunu(double &result)
Br B->D_s tau nu.
Definition: FlavBit.cpp:576
void RHN_mue_FF(const SMInputs sminputs, std::vector< double > &mnu, Eigen::Matrix< complex< double >, 3, 6 > &U, const double mH, complex< double > &g0SL, complex< double > &g0SR, complex< double > &g0VL, complex< double > &g0VR, complex< double > &g1SL, complex< double > &g1SR, complex< double > &g1VL, complex< double > &g1VR)
Definition: FlavBit.cpp:1878
Reader class for FlavBit YAML database.
Definition: Flav_reader.hpp:40
boost::numeric::ublas::matrix< double > get_exp_value()
Return the experimental central values.
Definition: Flav_reader.hpp:93
complex< double > CVLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void RHN_muePb(double &result)
Definition: FlavBit.cpp:2013
void read_yaml_measurement(str name, str measurement_name)
Read a single measurement from the database into memory.
void l2lll_likelihood(double &result)
Likelihood for l -> l l l processes.
Definition: FlavBit.cpp:2094
void RHN_RKstar_11_60(double &result)
Definition: FlavBit.cpp:956
void SI_BDmunu(double &result)
Br B -> D mu nu.
Definition: FlavBit.cpp:640
void FH_Bsmumu(double &result)
Definition: FlavBit.cpp:1106
void l2lgamma_likelihood(double &result)
Likelihood for l -> l gamma processes.
Definition: FlavBit.cpp:2046
void debug_mode(bool k)
Set debug mode for reader.
Definition: Flav_reader.hpp:87
void SL_measurements(predictions_measurements_covariances &pmc)
Measurements for tree-level leptonic and semileptonic B decays.
Definition: FlavBit.cpp:1454
void SI_Bmumu(double &result)
Br B0->mumu decays.
Definition: FlavBit.cpp:547
complex< double > K2R(int alpha, int beta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void SI_RKstar_0045_11(double &result)
RK* in low q^2.
Definition: FlavBit.cpp:894
double pow(const double &a)
Outputs a^i.
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum SM_spectrum
bool InvertMatrix(const ublas::matrix< T > &input, ublas::matrix< T > &inverse)
Matrix inversion routine using Boost.
Definition: flav_utils.hpp:36
void SI_Rmu(double &result)
B->K mu nu / B-> pi mu nu.
Definition: FlavBit.cpp:734
Declaration of reader class for FlavBit YAML database.
boost::numeric::ublas::matrix< double > cov_exp
void SI_BRBXsmumu_lowq2(double &result)
Inclusive branching fraction B -> X_s mu mu at low q^2.
Definition: FlavBit.cpp:776
void RHN_RK(double &result)
RK for RHN.
Definition: FlavBit.cpp:1005
void SI_BDstartaunu(double &result)
Br B -> D* tau nu.
Definition: FlavBit.cpp:662
T byVal(T t)
Redirection function to turn an lvalue into an rvalue, so that it is correctly passed by value when d...
def profile(file_name, frac_error=0.1, min_=0., max_=1., log_normal=True)
void RHN_mueTi(double &result)
Definition: FlavBit.cpp:1947
void SI_A_BXstautau_highq2(double &result)
Forward-backward asymmetry of B -> X_s tau tau at high q^2.
Definition: FlavBit.cpp:860
complex< double > BVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum Spectrum SMINPUTS
complex< double > BVRL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > ASLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)
void SI_RK(double &result)
RK between 1 and 6 GeV^2.
Definition: FlavBit.cpp:992
void FH_FlavourObs(fh_FlavourObs &result)
Flavour observables from FeynHiggs: B_s mass asymmetry, Br B_s -> mu mu, Br B -> X_s gamma...
Definition: FlavBit.cpp:1069
void SI_Btaunu(double &result)
Br B->tau nu_tau decays.
Definition: FlavBit.cpp:562
complex< double > BVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void SI_Dsmunu(double &result)
Br B->D_s mu nu.
Definition: FlavBit.cpp:590
TODO: see if we can use this one:
Definition: Analysis.hpp:33
double E(const double x, const double y)
boost::numeric::ublas::matrix< std::pair< double, bool > > get_th_err()
Return the (uncorrelated) theory errors.
Definition: Flav_reader.hpp:96
const bool flav_debug_LL
Definition: FlavBit.cpp:78
void RHN_mueee(double &result)
Definition: FlavBit.cpp:1749
complex< double > CVRR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
void RHN_taueee(double &result)
Definition: FlavBit.cpp:1768
Container class for Standard Model input information (defined as in SLHA2)
Definition: sminputs.hpp:29
complex< double > CVLR(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu)
complex< double > BSLL(int alpha, int beta, int gamma, int delta, SMInputs sminputs, Eigen::Matrix< complex< double >, 3, 6 > U, vector< double > ml, vector< double > mnu, double mh)