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