SUSYMatching Class Reference

A class for the matching in the MSSM. More...

#include <SUSYMatching.h>

Inheritance diagram for SUSYMatching:
[legend]
Collaboration diagram for SUSYMatching:
[legend]

Detailed Description

A class for the matching in the MSSM.

Author
HEPfit Collaboration

Definition at line 31 of file SUSYMatching.h.

Public Member Functions

virtual gslpp::vector< gslpp::complexAFunctions (int n)
 Calculates gamma penguin amplitudes for the process \( \ell_j \to \ell_i \gamma \) from [90] [14]. More...
 
virtual gslpp::vector< gslpp::complexBFunctions (int n)
 Calculates box diagram amplitudes for the process \( \ell_j \to \ell_i \ell_i \ell_i \) from from [90] [14]. More...
 
virtual gslpp::vector< gslpp::complexBHFunctions (int n)
 Calculates Higgs penguin amplitudes for the process \( \ell_j \to \ell_i \ell_i \ell_i \) from [14]. More...
 
virtual gslpp::vector< gslpp::complexC10_Lepton (int n)
 Calculates \( C_{10} \) and \( C^{\prime}_{10} \) Wilson coefficients for the process \( \ell_j \to \ell_i \ell_i \ell_i \). More...
 
virtual gslpp::vector< gslpp::complexC7_Lepton (int n)
 Calculates \( C_{7} \) and \( C^{\prime}_{7} \) Wilson coefficients for the process \( \ell_j \to \ell_i \gamma \). More...
 
virtual gslpp::vector< gslpp::complexC9_Lepton (int n)
 Calculates \( C_{9} \) and \( C^{\prime}_{9} \) Wilson coefficients for the process \( \ell_j \to \ell_i \ell_i \ell_i \). More...
 
virtual std::vector< WilsonCoefficient > & CMDLi3j (int li_lj)
 Wilson coefficient for the process \( \ell_j \to \ell_i \ell_i \ell_i \). More...
 
virtual std::vector< WilsonCoefficient > & CMDLij (int li_lj)
 Wilson coefficient for the process \( \ell_j \to \ell_i \gamma \). More...
 
virtual std::vector< WilsonCoefficient > & CMgminus2mu ()
 Wilson coefficient for \( (g-2)_{\mu} \) at one-loop. More...
 
virtual std::vector< WilsonCoefficient > & CMmueconv ()
 Wilson coefficient for the process \( \mu \to e \) conversion in Nuclei. More...
 
virtual gslpp::vector< gslpp::complexCP_Lepton (int n)
 Calculates \( C_{P} \) and \( C^{\prime}_{P} \) Wilson coefficients for the process \( \ell_j \to \ell_i \ell_i \ell_i \). More...
 
virtual gslpp::vector< gslpp::complexCS_Lepton (int n)
 Calculates \( C_{S} \) and \( C^{\prime}_{S} \) Wilson coefficients for the process \( \ell_j \to \ell_i \ell_i \ell_i \). More...
 
virtual gslpp::vector< gslpp::complexCT_Lepton (int n)
 Calculates \( C_{T} \) and \( C_{T5} \) Wilson coefficients for the process \( \ell_j \to \ell_i \ell_i \ell_i \). More...
 
virtual gslpp::vector< gslpp::complexDFunctions ()
 Calculates box diagram amplitudes for the process \( \mu \to e \) conversion in Nuclei from from [90]. More...
 
virtual gslpp::vector< gslpp::complexFFunctions (int n)
 Calculates Z penguin amplitudes for the process \( \ell_j \to \ell_i \ell_i \ell_i \) from [90] [14]. More...
 
virtual gslpp::vector< gslpp::complexgminus2mu ()
 Calculates amplitudes for \( (g-2)_{\mu} \) at one-loop from [90]. More...
 
 SUSYMatching (const SUSY &SUSY_i)
 
void updateSUSYParameters ()
 Updates to new SUSY parameter sets. More...
 
- Public Member Functions inherited from StandardModelMatching
double A0t (double x) const
 loop function which appear in the Wilson coefficient for the magnetic operator in the non-effective Misiak basis, Bobeth et al hep-ph/9910220 More...
 
double A1t (double x, double mu) const
 loop function which appear in the Wilson coefficient for the semileptonic operator in the non-effective Misiak basis, Misiak and Urban hep-ph/9901278v1 More...
 
double B0b (double x) const
 loop functions for non-leptonic B decays, Buiras Basis Buras et al, hep-ph/9512380v1 More...
 
double B0t (double x) const
 loop function which appear in the Wilson coefficient for the semileptonic operator in the non-effective Misiak basis, Bobeth et al hep-ph/9910220 More...
 
double B1t (double x, double mu) const
 loop function which appear in the Wilson coefficient for the semileptonic operator in the non-effective Misiak basis, Misiak and Urban hep-ph/9901278v1 More...
 
double C0b (double x) const
 loop functions for non-leptonic B decays, Buiras Basis Buras et al, hep-ph/9512380v1 More...
 
double C0t (double x) const
 loop function which appear in the Wilson coefficient for the magnetic operator in the non-effective Misiak basis, Bobeth et al hep-ph/9910220 More...
 
double C1t (double x, double mu) const
 loop function which appear in the Wilson coefficient for the magnetic operator in the non-effective Misiak basis, Misiak and Urban hep-ph/9901278v1 More...
 
double C7LOeff (double x) const
 loop function which appear in the Wilson coefficient for the magnetic operator in the effective Misiak basis, LO term, Chetyrkin et al hep-ph/9612313 More...
 
double C7NLOeff (double x) const
 loop function which appear in the Wilson coefficient for the magnetic operator in the effective Misiak basis, NLO term, Chetyrkin et al hep-ph/9612313 More...
 
double C8LOeff (double x) const
 loop function which appear in the Wilson coefficient for the chromomagnetic operator in the effective Misiak basis, LO term, Chetyrkin et al hep-ph/9612313 More...
 
double C8NLOeff (double x) const
 loop function which appear in the Wilson coefficient for the chromomagnetic operator in the effective Misiak basis, LO term, Chetyrkin et al hep-ph/9612313 More...
 
virtual std::vector< WilsonCoefficient > & CMbdmm ()
 
virtual std::vector< WilsonCoefficient > & CMBMll ()
 operator basis: current current; qcd penguins; magnetic and chromomagnetic penguins; semileptonic More...
 
virtual std::vector< WilsonCoefficient > & CMbnlep (int a)
 operator basis: More...
 
virtual std::vector< WilsonCoefficient > & CMbnlepCC (int a)
 operator basis: - current current opertors More...
 
virtual std::vector< WilsonCoefficient > & CMbsg ()
 operator basis: current current; qcd penguins; magnetic and chromomagnetic penguins; semileptonic More...
 
virtual std::vector< WilsonCoefficient > & CMbsmm ()
 
virtual std::vector< WilsonCoefficient > & CMbtaunu ()
 
virtual std::vector< WilsonCoefficient > & CMBXdnn ()
 
virtual std::vector< WilsonCoefficient > & CMBXsnn ()
 
virtual std::vector< WilsonCoefficient > & CMd1 ()
 current-current oerators, Misiak basis More...
 
virtual std::vector< WilsonCoefficient > & CMd1Buras ()
 current-current oerators, Buras basis More...
 
virtual std::vector< WilsonCoefficient > & CMdbd2 ()
 \( \Delta B = 2 \), \( B_{d} \) More...
 
virtual std::vector< WilsonCoefficient > & CMdbs2 ()
 \( \Delta B = 2 \), \( B_{s} \) More...
 
virtual std::vector< WilsonCoefficient > & CMdd2 ()
 \( \Delta C = 2 \), More...
 
virtual std::vector< WilsonCoefficient > & CMdk2 ()
 \( \Delta S = 2 \) More...
 
virtual std::vector< WilsonCoefficient > & CMK ()
 operator basis: More...
 
virtual std::vector< WilsonCoefficient > & CMKCC ()
 operator basis: More...
 
virtual std::vector< WilsonCoefficient > & CMkmm ()
 
virtual std::vector< WilsonCoefficient > & CMkpnn ()
 
virtual std::vector< WilsonCoefficient > & CMprimeBMll ()
 operator basis: current current; qcd penguins; magnetic and chromomagnetic penguins; semileptonic More...
 
virtual std::vector< WilsonCoefficient > & CMprimebsg ()
 operator basis: current current; qcd penguins; magnetic and chromomagnetic penguins; semileptonic More...
 
double D0b (double x) const
 loop functions for non-leptonic B decays, Buiras Basis Buras et al, hep-ph/9512380v1 More...
 
double D0t (double x) const
 loop function which appear in the Wilson coefficient for the magnetic operator in the non-effective Misiak basis, Bobeth et al hep-ph/9910220 More...
 
double D1t (double x, double mu) const
 loop function which appear in the Wilson coefficient for the magnetic operator in the non-effective Misiak basis, Misiak and Urban hep-ph/9901278v1 More...
 
double E0b (double x) const
 loop functions for non-leptonic B decays, Buiras Basis Buras et al, hep-ph/9512380v1 More...
 
double E0t (double x) const
 loop function which appear in the Wilson coefficient for the chromomagnetic operator in the Misiak basis, Chetyrkin et al hep-ph/9612313 More...
 
double E1t (double x, double mu) const
 loop function which appear in the Wilson coefficient in the non-effective Misiak basis, Misiak and Urban hep-ph/9910220 More...
 
double Eet (double x) const
 loop function which appear in the Wilson coefficient in the non-effective Misiak basis, Misiak and Urban hep-ph/0512066 More...
 
double F0t (double x) const
 loop function which appear in the Wilson coefficient for the chromomagnetic operator in the non-effective Misiak basis, Bobeth et al hep-ph/9910220 More...
 
double F1t (double x, double mu) const
 loop function which appear in the Wilson coefficient for the semileptonic operator in the non-effective Misiak basis, Misiak and Urban hep-ph/9901278v1 More...
 
double G1t (double x, double mu) const
 loop function which appear in the Wilson coefficient in the non-effective Misiak basis, Misiak and Urban hep-ph/9910220 More...
 
double mt2omh2 (const double mu, const orders order=FULLNNLO) const
 
double Rest (double x, double mu) const
 approximation of two-loops EW correction for Q_10 operator in the non-effective Misiak basis, Misiak and Urban hep-ph/1311.1348 More...
 
double S0 (double, double) const
 
gslpp::complex S0c () const
 hep-ph/9512380 More...
 
gslpp::complex S0ct () const
 hep-ph/9512380 More...
 
gslpp::complex S0tt () const
 hep-ph/9512380v1 More...
 
double S1 (double x) const
 
 StandardModelMatching (const StandardModel &SM_i)
 
double Tt (double x) const
 loop function which appear in the Wilson coefficient in the non-effective Misiak basis, Misiak and Urban hep-ph/9910220 More...
 
void updateSMParameters ()
 Updates to new Standard Model parameter sets. More...
 
double Wt (double x) const
 loop function which appear in the Wilson coefficient in the non-effective Misiak basis, Misiak and Urban hep-ph/0512066 More...
 
double X0t (double x) const
 hep-ph/9512380v1 More...
 
double X1t (double x) const
 hep-ph/1009.0947v2 More...
 
double x_c (const double mu, const orders order=FULLNNLO) const
 
double x_t (const double mu, const orders order=FULLNNLO) const
 
double Xewt (double x, double a, double mu) const
 hep-ph/1009.0947v2 More...
 
double Y0 (double x) const
 
double Y1 (double x, double mu) const
 
- Public Member Functions inherited from ModelMatching
virtual ~ModelMatching ()
 

Private Member Functions

int delta_ab (int a, int b)
 Kronecker delta. More...
 
void NeutralinoRemixing ()
 

Private Attributes

gslpp::matrix< gslpp::complexAmpA1LC
 
gslpp::matrix< gslpp::complexAmpA1LN
 Amplitudes of Chargino and Neutralino contribution to various LFV observables. More...
 
gslpp::matrix< gslpp::complexAmpA1RC
 
gslpp::matrix< gslpp::complexAmpA1RN
 
gslpp::matrix< gslpp::complexAmpALC
 
gslpp::matrix< gslpp::complexAmpALN
 
gslpp::matrix< gslpp::complexAmpARC
 
gslpp::matrix< gslpp::complexAmpARN
 
gslpp::matrix< gslpp::complexAmpTauA1LC
 
gslpp::matrix< gslpp::complexAmpTauA1LN
 
gslpp::matrix< gslpp::complexAmpTauA1RC
 
gslpp::matrix< gslpp::complexAmpTauA1RN
 
gslpp::matrix< gslpp::complexAmpTauALC
 
gslpp::matrix< gslpp::complexAmpTauALN
 
gslpp::matrix< gslpp::complexAmpTauARC
 
gslpp::matrix< gslpp::complexAmpTauARN
 
gslpp::matrix< gslpp::complexAmpTEA1LC
 
gslpp::matrix< gslpp::complexAmpTEA1LN
 
gslpp::matrix< gslpp::complexAmpTEA1RC
 
gslpp::matrix< gslpp::complexAmpTEA1RN
 
gslpp::matrix< gslpp::complexAmpTEALC
 
gslpp::matrix< gslpp::complexAmpTEALN
 
gslpp::matrix< gslpp::complexAmpTEARC
 
gslpp::matrix< gslpp::complexAmpTEARN
 
gslpp::matrix< gslpp::complexCLlE
 
gslpp::matrix< gslpp::complexCLlMU
 
gslpp::matrix< gslpp::complexCLlTAU
 
gslpp::matrix< gslpp::complexCLqDOWN
 
gslpp::matrix< gslpp::complexCLqUP
 
double cosb
 
gslpp::matrix< gslpp::complexCRlE
 Chargino and Neutralino couplings to sfermions. More...
 
gslpp::matrix< gslpp::complexCRlMU
 
gslpp::matrix< gslpp::complexCRlTAU
 
gslpp::matrix< gslpp::complexCRqDOWN
 
gslpp::matrix< gslpp::complexCRqUP
 
double gW
 
gslpp::matrix< double > Leptf1
 
gslpp::matrix< double > Leptf2
 
gslpp::matrix< double > Leptf3
 
gslpp::matrix< double > Leptf4
 
gslpp::matrix< double > Leptfa1
 
gslpp::matrix< double > Leptfa2
 
std::complex< double > Leptfzc [3][2][2]
 
std::complex< double > Leptfzn [6][4][4]
 
std::complex< double > Leptgzc [3][2][2]
 
std::complex< double > Leptgzn [6][4][4]
 
gslpp::matrix< double > Lepty
 Functions needed to calculate various LFV observables. More...
 
gslpp::matrix< double > Leptz
 
WilsonCoefficient mcDLi3j
 
WilsonCoefficient mcDLij
 
WilsonCoefficient mcgminus2mu
 
gslpp::vector< double > MChi
 Chargino mass-eigenvalue. More...
 
gslpp::vector< double > MChi0
 Neutralino mass-eigenvalue. More...
 
WilsonCoefficient mcmueconv
 
gslpp::vector< double > MNeig
 
gslpp::vector< double > mym_sd_sq
 
gslpp::vector< double > mym_se_sq
 
gslpp::vector< double > mym_sn_sq
 
gslpp::vector< double > mym_su_sq
 Sfermion mass-eigenvalue squared. More...
 
gslpp::matrix< gslpp::complexmyN
 Neutralino mixing matrix. More...
 
gslpp::matrix< gslpp::complexmyRd
 
gslpp::matrix< gslpp::complexmyRl
 
gslpp::matrix< gslpp::complexmyRn
 
gslpp::matrix< gslpp::complexmyRu
 Sfermion mixing matrices. More...
 
const SUSYmySUSY
 
gslpp::matrix< gslpp::complexmyU
 
gslpp::matrix< gslpp::complexmyV
 Chargino mixing matrices. More...
 
gslpp::matrix< gslpp::complexNLlE
 
gslpp::matrix< gslpp::complexNLlMU
 
gslpp::matrix< gslpp::complexNLlTAU
 
gslpp::matrix< gslpp::complexNLqDOWN
 
gslpp::matrix< gslpp::complexNLqUP
 
gslpp::matrix< gslpp::complexNRlE
 
gslpp::matrix< gslpp::complexNRlMU
 
gslpp::matrix< gslpp::complexNRlTAU
 
gslpp::matrix< gslpp::complexNRqDOWN
 
gslpp::matrix< gslpp::complexNRqUP
 
gslpp::matrix< gslpp::complexON
 
const PVfunctions PV
 
double sinb
 
double tanb
 
gslpp::matrix< gslpp::complexTEhat
 Slepton tri-linear coupling matrix. More...
 
double v
 
double v1
 
double v2
 

Additional Inherited Members

- Protected Attributes inherited from StandardModelMatching
std::vector< WilsonCoefficientvmcbdmm
 
std::vector< WilsonCoefficientvmcbdnn
 
std::vector< WilsonCoefficientvmcBMll
 
std::vector< WilsonCoefficientvmcbnlep
 
std::vector< WilsonCoefficientvmcbnlepCC
 
std::vector< WilsonCoefficientvmcbsg
 
std::vector< WilsonCoefficientvmcbsmm
 
std::vector< WilsonCoefficientvmcbsnn
 
std::vector< WilsonCoefficientvmcbtaunu
 
std::vector< WilsonCoefficientvmcd1
 
std::vector< WilsonCoefficientvmcd1Buras
 
std::vector< WilsonCoefficientvmcd2
 
std::vector< WilsonCoefficientvmcdb
 
std::vector< WilsonCoefficientvmcDLi3j
 
std::vector< WilsonCoefficientvmcDLij
 
std::vector< WilsonCoefficientvmcds
 
std::vector< WilsonCoefficientvmcgminus2mu
 
std::vector< WilsonCoefficientvmck
 
std::vector< WilsonCoefficientvmck2
 
std::vector< WilsonCoefficientvmckcc
 
std::vector< WilsonCoefficientvmckmm
 
std::vector< WilsonCoefficientvmckpnn
 
std::vector< WilsonCoefficientvmcmueconv
 
std::vector< WilsonCoefficientvmcprimeBMll
 
std::vector< WilsonCoefficientvmcprimebsg
 

Constructor & Destructor Documentation

SUSYMatching::SUSYMatching ( const SUSY SUSY_i)

Definition at line 13 of file SUSYMatching.cpp.

13  :
14 
15  StandardModelMatching(SUSY_i),
16  mySUSY(SUSY_i),
17  PV(true),
18  mcDLij(2, NDR, LO),
19  mcDLi3j(20, NDR, LO),
20  mcmueconv(8, NDR, LO),
21  mcgminus2mu(2, NDR, LO),
22 
23  mym_su_sq(6, 0.),
24  mym_sd_sq(6, 0.),
25  mym_se_sq(6, 0.),
26  mym_sn_sq(6, 0.),
27 
28  myRu(6, 6, 0.),
29  myRd(6, 6, 0.),
30  myRl(6, 6, 0.),
31  myRn(6, 6, 0.),
32 
33  MChi(2, 0.),
34  myV(2, 2, 0.),
35  myU(2, 2, 0.),
36  MChi0(4, 0.),
37  MNeig(4, 0.),
38  myN(4, 4, 0.),
39  ON(4, 4, 0.),
40 
41  Lepty(4, 6, 0.),
42  Leptz(2, 3, 0.),
43  Leptfa1(4, 6, 0.),
44  Leptfa2(2, 3, 0.),
45  Leptf1(4, 6, 0.),
46  Leptf2(4, 6, 0.),
47  Leptf3(2, 3, 0.),
48  Leptf4(2, 3, 0.),
49  CRlE(2, 3, 0.),
50  CRlMU(2, 3, 0.),
51  CRlTAU(2, 3, 0.),
52  CRqUP(2, 6, 0.),
53  CRqDOWN(2, 6, 0.),
54  CLlE(2, 3, 0.),
55  CLlMU(2, 3, 0.),
56  CLlTAU(2, 3, 0.),
57  CLqUP(2, 6, 0.),
58  CLqDOWN(2, 6, 0.),
59  NRlE(4, 6, 0.),
60  NRlMU(4, 6, 0.),
61  NRlTAU(4, 6, 0.),
62  NRqUP(4, 6, 0.),
63  NRqDOWN(4, 6, 0.),
64  NLlE(4, 6, 0.),
65  NLlMU(4, 6, 0.),
66  NLlTAU(4, 6, 0.),
67  NLqUP(4, 6, 0.),
68  NLqDOWN(4, 6, 0.),
69  AmpA1LN(4, 6, 0.),
70  AmpA1RN(4, 6, 0.),
71  AmpA1LC(2, 3, 0.),
72  AmpA1RC(2, 3, 0.),
73  AmpTauA1LN(4, 6, 0.),
74  AmpTauA1RN(4, 6, 0.),
75  AmpTauA1LC(2, 3, 0.),
76  AmpTauA1RC(2, 3, 0.),
77  AmpTEA1LN(4, 6, 0.),
78  AmpTEA1RN(4, 6, 0.),
79  AmpTEA1LC(2, 3, 0.),
80  AmpTEA1RC(2, 3, 0.),
81  AmpALN(4, 6, 0.),
82  AmpARN(4, 6, 0.),
83  AmpALC(2, 3, 0.),
84  AmpARC(2, 3, 0.),
85  AmpTauALN(4, 6, 0.),
86  AmpTauARN(4, 6, 0.),
87  AmpTauALC(2, 3, 0.),
88  AmpTauARC(2, 3, 0.),
89  AmpTEALN(4, 6, 0.),
90  AmpTEARN(4, 6, 0.),
91  AmpTEALC(2, 3, 0.),
92  AmpTEARC(2, 3, 0.),
93 
94  TEhat(3, 3, 0.)
95 {
96 }
WilsonCoefficient mcgminus2mu
Definition: SUSYMatching.h:180
gslpp::matrix< gslpp::complex > AmpTauA1LN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > myRl
Definition: SUSYMatching.h:190
gslpp::matrix< gslpp::complex > AmpTauALN
Definition: SUSYMatching.h:225
WilsonCoefficient mcmueconv
Definition: SUSYMatching.h:180
gslpp::matrix< gslpp::complex > AmpTEARN
Definition: SUSYMatching.h:225
gslpp::vector< double > mym_se_sq
Definition: SUSYMatching.h:185
gslpp::matrix< gslpp::complex > NLlE
Definition: SUSYMatching.h:220
gslpp::vector< double > MChi
Chargino mass-eigenvalue.
Definition: SUSYMatching.h:195
gslpp::vector< double > mym_sn_sq
Definition: SUSYMatching.h:185
gslpp::vector< double > mym_sd_sq
Definition: SUSYMatching.h:185
gslpp::matrix< gslpp::complex > NLqUP
Definition: SUSYMatching.h:220
gslpp::matrix< double > Lepty
Functions needed to calculate various LFV observables.
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > ON
Definition: SUSYMatching.h:209
gslpp::matrix< double > Leptf4
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > CRlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< double > Leptfa1
Definition: SUSYMatching.h:214
WilsonCoefficient mcDLij
Definition: SUSYMatching.h:180
gslpp::matrix< double > Leptf2
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > AmpTEA1LC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CLlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpA1LC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CLqDOWN
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpTauALC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > NRlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpA1LN
Amplitudes of Chargino and Neutralino contribution to various LFV observables.
Definition: SUSYMatching.h:225
gslpp::vector< double > MNeig
Definition: SUSYMatching.h:205
gslpp::matrix< gslpp::complex > myRd
Definition: SUSYMatching.h:190
gslpp::matrix< gslpp::complex > CLqUP
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpTauA1RN
Definition: SUSYMatching.h:225
gslpp::matrix< double > Leptz
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > AmpTEALN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpARN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpA1RN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CRlMU
Definition: SUSYMatching.h:220
const PVfunctions PV
Definition: SUSYMatching.h:178
gslpp::matrix< gslpp::complex > myRn
Definition: SUSYMatching.h:190
gslpp::matrix< gslpp::complex > AmpTEARC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > NRlE
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NRqDOWN
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > myRu
Sfermion mixing matrices.
Definition: SUSYMatching.h:190
gslpp::matrix< gslpp::complex > AmpTEA1RC
Definition: SUSYMatching.h:225
WilsonCoefficient mcDLi3j
Definition: SUSYMatching.h:180
gslpp::matrix< gslpp::complex > myV
Chargino mixing matrices.
Definition: SUSYMatching.h:200
gslpp::matrix< double > Leptf3
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > AmpTauARN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpARC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > NRqUP
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > CRlE
Chargino and Neutralino couplings to sfermions.
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpTEA1RN
Definition: SUSYMatching.h:225
gslpp::vector< double > mym_su_sq
Sfermion mass-eigenvalue squared.
Definition: SUSYMatching.h:185
Definition: OrderScheme.h:33
gslpp::matrix< gslpp::complex > CRqDOWN
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > CLlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > CRqUP
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NLlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NRlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > TEhat
Slepton tri-linear coupling matrix.
Definition: SUSYMatching.h:230
gslpp::matrix< double > Leptf1
Definition: SUSYMatching.h:214
gslpp::matrix< double > Leptfa2
Definition: SUSYMatching.h:214
StandardModelMatching(const StandardModel &SM_i)
gslpp::matrix< gslpp::complex > AmpTEALC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpA1RC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpTauA1RC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CLlE
Definition: SUSYMatching.h:220
gslpp::vector< double > MChi0
Neutralino mass-eigenvalue.
Definition: SUSYMatching.h:205
gslpp::matrix< gslpp::complex > NLqDOWN
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NLlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpALN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpALC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpTauARC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > myU
Definition: SUSYMatching.h:200
gslpp::matrix< gslpp::complex > myN
Neutralino mixing matrix.
Definition: SUSYMatching.h:209
gslpp::matrix< gslpp::complex > AmpTauA1LC
Definition: SUSYMatching.h:225
const SUSY & mySUSY
Definition: SUSYMatching.h:177
gslpp::matrix< gslpp::complex > AmpTEA1LN
Definition: SUSYMatching.h:225

Member Function Documentation

gslpp::vector< gslpp::complex > SUSYMatching::AFunctions ( int  n)
virtual

Calculates gamma penguin amplitudes for the process \( \ell_j \to \ell_i \gamma \) from [90] [14].

Calculates gamma penguin amplitudes for m->(3)e (1), t->(3)m (2) and t->(3)e (3)

Parameters
[in]ndetermines the process, e.g., 1 = \( \mu \to e \gamma \), 2 = \( \tau \to \mu \gamma \), 3 = \( \tau \to e \gamma \)
Returns
returns the vector of gamma penguin amplitude

Definition at line 147 of file SUSYMatching.cpp.

147  {
148  //gamma penguin contributions
149 
151 
152  double MW = mySUSY.Mw_tree();
153  double pi = M_PI;
154  double piconst = 1.0/(32.0 * pi * pi);
155  double sw2 = mySUSY.StandardModel::sW2(MW);
156  double stw = sqrt(sw2);
157  double ctw = sqrt(1.0 - sw2);
158  double ttw = stw/ctw;
160  double mMU = mySUSY.getLeptons(StandardModel::MU).getMass();
161  double mTAU = mySUSY.getLeptons(StandardModel::TAU).getMass();
162  sinb = mySUSY.getSinb();
163  double cdenc = sqrt(2.0)*MW*cosb;
164  double cdenn = MW*cosb;
165  double g2 = gW;
166  double g2t = g2/sqrt(2.0);
167 
169 
170  // Neutralino-Fermion-Sfermion couplings
171  for (int a=0;a<4;a++) {
172  for (int x=0;x<6;x++) {
173  // LL + RL TYPE MI
174  NRlE.assign(a, x, - (g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 0) + (mE/cdenn)*ON(a, 2)*myRl(x, 3)));
175  NRlMU.assign(a, x, -(g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 1) + (mMU/cdenn)*ON(a, 2)*myRl(x, 4)));
176  NRlTAU.assign(a, x, -(g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 2) + (mTAU/cdenn)*ON(a, 2)*myRl(x, 5)));
177  // RL + RR TYPE MI
178  NLlE.assign(a, x, -(g2t)*((mE/cdenn)*ON(a, 2)*myRl(x, 0) + 2.0*ON(a, 0)*ttw*myRl(x, 3)));
179  NLlMU.assign(a, x, -(g2t)*((mMU/cdenn)*ON(a, 2)*myRl(x, 1) + 2.0*ON(a, 0)*ttw*myRl(x, 4)));
180  NLlTAU.assign(a, x, -(g2t)*((mTAU/cdenn)*ON(a, 2)*myRl(x, 2) + 2.0*ON(a, 0)*ttw*myRl(x, 5)));
181 // Commented expressions might be useful for complex neutralino mixing matrices
182 // NLlE.assign(a, x, -(g2t)*((mE/cdenn)*ON(a, 2).conjugate()*myRl(x, 0) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 3)));
183 // NLlMU.assign(a, x, -(g2t)*((mMU/cdenn)*ON(a, 2).conjugate()*myRl(x, 1) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 4)));
184 // NLlTAU.assign(a, x, -(g2t)*((mTAU/cdenn)*ON(a, 2).conjugate()*myRl(x, 2) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 5)));
185  }
186  }
187 
188  // Chargino-Fermion-Sfermion couplings
189  for (int a=0;a<2;a++) {
190  for (int x=0;x<3;x++) {
191  // LL-TYPE
192  CRlE.assign(a, x, - (g2*myV(a, 0)*myRn(x, 0)));
193  CRlMU.assign(a, x, - (g2*myV(a, 0)*myRn(x, 1)));
194  CRlTAU.assign(a, x, - (g2*myV(a, 0)*myRn(x, 2)));
195  // LR-TYPE
196  CLlE.assign(a, x, g2*mE/cdenc*myU(a, 1).conjugate()*myRn(x, 0));
197  CLlMU.assign(a, x, g2*mMU/cdenc*myU(a, 1).conjugate()*myRn(x, 1));
198  CLlTAU.assign(a, x, g2*mTAU/cdenc*myU(a, 1).conjugate()*myRn(x, 2));
199  }
200  }
201 
202  // Definition of y and z - remember they are dimensionless quantities.
203  for (int a=0;a<4;a++) {
204  for (int x=0;x<6;x++) {
205  Lepty.assign(a, x, MNeig(a) * MNeig(a) / mym_se_sq(x) );
206  }
207  }
208 
209  for (int a=0;a<2;a++) {
210  for (int x=0;x<3;x++) {
211  Leptz.assign(a, x, MChi(a) * MChi(a) / mym_sn_sq(x) );
212  }
213  }
214 
215  for (int a=0;a<4;a++) {
216  for (int x=0;x<6;x++) {
217  if (fabs(1.0 - Lepty(a, x)) > 0.01) {
218  Leptfa1.assign(a, x, (1.0/mym_se_sq(x))*(1.0/pow((1.0 - Lepty(a, x)),4.0))*
219  (2.0 - 9.0*Lepty(a, x) + 18.0*pow(Lepty(a, x),2.0) - 11.0*pow(Lepty(a, x),3.0)
220  + 6.0*pow(Lepty(a, x),3.0)*log(Lepty(a, x))) );
221  Leptf1.assign(a, x, ((1.0 - 6.0*Lepty(a, x) + 3.0 * pow(Lepty(a, x),2.0) +
222  2.0*pow(Lepty(a, x),3.0) - 6.0*pow(Lepty(a,x),2.0)*log(Lepty(a, x))))/
223  (6.0 * pow((1.0 - Lepty(a,x)),4.0)) );
224  Leptf2.assign(a, x, (1.0 - pow(Lepty(a, x),2.0) + 2.0 * Lepty(a, x) * log(Lepty(a, x)))/
225  (pow((1.0-Lepty(a, x)),3.0)));
226  }
227  else {
228  Leptfa1.assign(a, x, (3.0/2.0)*(1.0/mym_se_sq(x)));
229  Leptf1.assign(a, x, 1.0/12.0 - (Lepty(a, x) - 1.0)/30.0);
230  Leptf2.assign(a, x, 1.0/3.0 - (Lepty(a, x) - 1.0)/6.0);
231  }
232  }
233  }
234 
235  for (int a=0;a<2;a++) {
236  for (int x=0;x<3;x++) {
237  if(fabs(1.0-Leptz(a, x)) > 0.01) {
238  Leptfa2.assign(a, x, (1.0/mym_sn_sq(x))*(1.0/pow(1.0 - Leptz(a,x),4.0))*
239  (16.0 - 45.0*Leptz(a,x) + 36.0*pow(Leptz(a,x),2.0) - 7.0*pow(Leptz(a,x),3.0)
240  + 6.0*(2.0 - 3.0*Leptz(a,x))*log(Leptz(a,x))) );
241  Leptf3.assign(a, x, ((2.0 + 3.0*Leptz(a, x) - 6.0*pow(Leptz(a, x),2.0)
242  + pow(Leptz(a, x),3.0) + 6.0*Leptz(a, x)*log(Leptz(a, x)))/
243  (6.0*pow((1.0 - Leptz(a, x)),4.0))) );
244  Leptf4.assign(a, x, ((-3.0 + 4.0*Leptz(a, x) - pow(Leptz(a, x),2.0)
245  - 2.0*log(Leptz(a, x)))/
246  pow((1.0 - Leptz(a, x)),3.0)) );
247  }
248  else {
249  Leptfa2.assign(a, x, (-9.0/2.0)*(1.0/mym_sn_sq(x)) );
250  Leptf3.assign(a, x, 1.0/12.0 - (Leptz(a, x) - 1.0)/20.0 );
251  Leptf4.assign(a, x, 2.0/3.0 - (Leptz(a, x) - 1.0)/2.0 );
252  }
253  }
254  }
255 
256  if (li_to_lj == 1) // mu -> (3)e
257  {
258  // Neutralino contributions
259  for (int a=0;a<4;a++) {
260  for (int x=0;x<6;x++) {
261  AmpA1RN.assign(a, x, (piconst/18.0)*NLlE(a,x)*NLlMU(a,x).conjugate()*Leptfa1(a,x) );
262  AmpA1LN.assign(a, x, (piconst/18.0)*NRlE(a,x)*NRlMU(a,x).conjugate()*Leptfa1(a,x) );
263  AmpARN.assign(a, x, piconst*(NRlE(a, x) * NRlMU(a, x).conjugate() * Leptf1(a, x)
264 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
265  + NLlE(a, x) * NLlMU(a, x).conjugate() * (mE/mMU) * Leptf1(a, x)
266 //...until here
267  + NRlE(a, x) * NLlMU(a, x).conjugate() * (MNeig(a)/mMU) * Leptf2(a, x)) /mym_se_sq(x) );
268  AmpALN.assign(a, x, piconst*(NLlE(a, x) * NLlMU(a, x).conjugate() * Leptf1(a, x)
269 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
270  + NRlE(a, x) * NRlMU(a, x).conjugate() * (mE/mMU) * Leptf1(a, x)
271 //...until here
272  + NLlE(a, x) * NRlMU(a, x).conjugate() * (MNeig(a)/mMU) * Leptf2(a, x)) /mym_se_sq(x) );
273  }
274  }
275  gslpp::complex A1RN = 0.0;
276  gslpp::complex A1LN = 0.0;
277  gslpp::complex ARN = 0.0;
278  gslpp::complex ALN = 0.0;
279  for (int a=0;a<4;a++) {
280  for (int x=0;x<6;x++) {
281  A1RN = A1RN + AmpA1RN(a,x);
282  A1LN = A1LN + AmpA1LN(a,x);
283  ARN = ARN + AmpARN(a,x);
284  ALN = ALN + AmpALN(a,x);
285  }
286  }
287  // Chargino contributions
288  for (int a=0;a<2;a++) {
289  for (int x=0;x<3;x++) {
290  AmpA1RC.assign(a, x, -(piconst/18.0)*CLlE(a,x)*CLlMU(a,x).conjugate()*Leptfa2(a,x) );
291  AmpA1LC.assign(a, x, -(piconst/18.0)*CRlE(a,x)*CRlMU(a,x).conjugate()*Leptfa2(a,x) );
292  AmpARC.assign(a, x, -(piconst/mym_sn_sq(x)) * (CRlE(a, x)*CRlMU(a, x).conjugate() * Leptf3(a, x)
293 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
294  + CLlE(a, x) * CLlMU(a, x).conjugate() * (mE/mMU) * Leptf3(a, x)
295 //...until here
296  + CRlE(a, x) * CLlMU(a, x).conjugate() * (MChi(a)/mMU) * Leptf4(a, x)) );
297  AmpALC.assign(a, x, -(piconst/mym_sn_sq(x)) * (CLlE(a, x) * CLlMU(a, x).conjugate() * Leptf3(a, x)
298 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
299  + CRlE(a, x) * CRlMU(a, x).conjugate() * (mE/mMU) * Leptf3(a, x)
300 //...until here
301  + CLlE(a, x) * CRlMU(a, x).conjugate() * (MChi(a)/mMU) * Leptf4(a, x)) );
302  }
303  }
304  gslpp::complex A1RC = 0.0;
305  gslpp::complex A1LC = 0.0;
306  gslpp::complex ARC = 0.0;
307  gslpp::complex ALC = 0.0;
308  for (int a=0;a<2;a++) {
309  for (int x=0;x<3;x++) {
310  A1RC = A1RC + AmpA1RC(a,x);
311  A1LC = A1LC + AmpA1LC(a,x);
312  ARC = ARC + AmpARC(a,x);
313  ALC = ALC + AmpALC(a,x);
314  }
315  }
316  // write AR and AL into a vector for mu->(3)e
317  AFunctions.assign(0, A1RN + A1RC );
318  AFunctions.assign(1, A1LN + A1LC );
319  AFunctions.assign(2, ARN + ARC );
320  AFunctions.assign(3, ALN + ALC );
321 
322  }
323 
324  if (li_to_lj == 2) // tau -> (3)mu
325  {
326  // Neutralino contributions
327  for (int a=0;a<4;a++) {
328  for (int x=0;x<6;x++) {
329  AmpTauA1RN.assign(a, x, (piconst/18.0)*NLlMU(a,x)*NLlTAU(a,x).conjugate()*Leptfa1(a,x) );
330  AmpTauA1LN.assign(a, x, (piconst/18.0)*NRlMU(a,x)*NRlTAU(a,x).conjugate()*Leptfa1(a,x) );
331  AmpTauARN.assign(a, x, piconst * (NRlMU(a, x) * NRlTAU(a, x).conjugate() * Leptf1(a, x)
332 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
333  + NLlMU(a, x) * NLlTAU(a, x).conjugate() * (mMU/mTAU) * Leptf1(a, x)
334 //...until here
335  + NRlMU(a, x) * NLlTAU(a, x).conjugate() * (MNeig(a)/mTAU) * Leptf2(a, x)) /mym_se_sq(x) );
336  AmpTauALN.assign(a, x, piconst * (NLlMU(a, x) * NLlTAU(a, x).conjugate() * Leptf1(a, x)
337 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
338  + NRlMU(a, x) * NRlTAU(a, x).conjugate() * (mMU/mTAU) * Leptf1(a, x)
339 //...until here
340  + NLlMU(a, x) * NRlTAU(a, x).conjugate() * (MNeig(a)/mTAU) * Leptf2(a, x)) /mym_se_sq(x) );
341  }
342  }
343  gslpp::complex TauA1RN = 0.0;
344  gslpp::complex TauA1LN = 0.0;
345  gslpp::complex TauARN = 0.0;
346  gslpp::complex TauALN = 0.0;
347  for (int a=0;a<4;a++) {
348  for (int x=0;x<6;x++) {
349  TauA1RN = TauA1RN + AmpTauA1RN(a,x);
350  TauA1LN = TauA1LN + AmpTauA1LN(a,x);
351  TauARN = TauARN + AmpTauARN(a,x);
352  TauALN = TauALN + AmpTauALN(a,x);
353  }
354  }
355  // Chargino contributions
356  for (int a=0;a<2;a++) {
357  for (int x=0;x<3;x++) {
358  AmpTauA1RC.assign(a, x, -(piconst/18.0)*CLlMU(a,x)*CLlTAU(a,x).conjugate()*Leptfa2(a,x) );
359  AmpTauA1LC.assign(a, x, -(piconst/18.0)*CRlMU(a,x)*CRlTAU(a,x).conjugate()*Leptfa2(a,x) );
360  AmpTauARC.assign(a, x, -piconst / mym_sn_sq(x) * (CRlMU(a, x) * CRlTAU(a, x).conjugate() * Leptf3(a, x)
361 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
362  + CLlMU(a, x) * CLlTAU(a, x).conjugate() * (mMU/mTAU) * Leptf3(a, x)
363 //...until here
364  + CRlMU(a, x) * CLlTAU(a, x).conjugate() * (MChi(a)/mTAU) * Leptf4(a, x)) );
365  AmpTauALC.assign(a, x, -piconst / mym_sn_sq(x) * (CLlMU(a, x) * CLlTAU(a, x).conjugate() * Leptf3(a, x)
366 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
367  + CRlMU(a, x) * CRlTAU(a, x).conjugate() * (mMU/mTAU) * Leptf3(a, x)
368 //...until here
369  + CLlMU(a, x) * CRlTAU(a, x).conjugate() * (MChi(a)/mTAU) * Leptf4(a, x)) );
370  }
371  }
372  gslpp::complex TauA1RC = 0.0;
373  gslpp::complex TauA1LC = 0.0;
374  gslpp::complex TauARC = 0.0;
375  gslpp::complex TauALC = 0.0;
376  for (int a=0;a<2;a++) {
377  for (int x=0;x<3;x++) {
378  TauA1RC = TauA1RC + AmpTauA1RC(a,x);
379  TauA1LC = TauA1LC + AmpTauA1LC(a,x);
380  TauARC = TauARC + AmpTauARC(a,x);
381  TauALC = TauALC + AmpTauALC(a,x);
382  }
383  }
384  // write AR and AL into a vector for tau->(3)mu
385  AFunctions.assign(0, TauA1RC + TauA1RN );
386  AFunctions.assign(1, TauA1LC + TauA1LN );
387  AFunctions.assign(2, TauARC + TauARN );
388  AFunctions.assign(3, TauALC + TauALN );
389  }
390 
391  if (li_to_lj == 3) // tau -> (3)e
392  {
393  // Neutralino contributions
394  for (int a=0;a<4;a++) {
395  for (int x=0;x<6;x++) {
396  AmpTEA1RN.assign(a, x, (piconst/18.0)*NLlE(a,x)*NLlTAU(a,x).conjugate()*Leptfa1(a,x) );
397  AmpTEA1LN.assign(a, x, (piconst/18.0)*NRlE(a,x)*NRlTAU(a,x).conjugate()*Leptfa1(a,x) );
398  AmpTEARN.assign(a, x, piconst * (NRlE(a, x) * NRlTAU(a, x).conjugate() * Leptf1(a,x)
399 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
400  + NLlE(a, x) * NLlTAU(a, x).conjugate() * (mE/mTAU) * Leptf1(a,x)
401 //...until here
402  + NRlE(a, x) * NLlTAU(a, x).conjugate() * (MNeig(a)/mTAU) * Leptf2(a, x)) / mym_se_sq(x) );
403  AmpTEALN.assign(a, x, piconst * (NLlE(a, x) * NLlTAU(a, x).conjugate() * Leptf1(a, x)
404 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
405  + NRlE(a, x) * NRlTAU(a, x).conjugate() * (mE/mTAU) * Leptf1(a, x)
406 //...until here
407  + NLlE(a, x) * NRlTAU(a, x).conjugate() * (MNeig(a)/mTAU) * Leptf2(a, x)) / mym_se_sq(x) );
408  }
409  }
410  gslpp::complex TEA1RN = 0.0;
411  gslpp::complex TEA1LN = 0.0;
412  gslpp::complex TEARN = 0.0;
413  gslpp::complex TEALN = 0.0;
414  for (int a=0;a<4;a++) {
415  for (int x=0;x<6;x++) {
416  TEA1RN = TEA1RN + AmpTEA1RN(a,x);
417  TEA1LN = TEA1LN + AmpTEA1LN(a,x);
418  TEARN = TEARN + AmpTEARN(a,x);
419  TEALN = TEALN + AmpTEALN(a,x);
420  }
421  }
422  // Chargino contributions
423  for (int a=0;a<2;a++) {
424  for (int x=0;x<3;x++) {
425  AmpTEA1RC.assign(a, x, -(piconst/18.0)*CLlE(a,x)*CLlTAU(a,x).conjugate()*Leptfa2(a,x) );
426  AmpTEA1LC.assign(a, x, -(piconst/18.0)*CRlE(a,x)*CRlTAU(a,x).conjugate()*Leptfa2(a,x) );
427  AmpTEARC.assign(a, x, -piconst / mym_sn_sq(x) * (CRlE(a, x) * CRlTAU(a, x).conjugate() * Leptf3(a, x)
428 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
429  + CLlE(a, x) * CLlTAU(a, x).conjugate() * (mE/mTAU) * Leptf3(a, x)
430 //...until here
431  + CRlE(a, x) * CLlTAU(a, x).conjugate() *(MChi(a)/mTAU) * Leptf4(a, x)) );
432  AmpTEALC.assign(a, x, -piconst / mym_sn_sq(x) * (CLlE(a, x) * CLlTAU(a, x).conjugate() * Leptf3(a, x)
433 // The following contribution is absent in PRD 53.2442 (Hisano et al.), but appears in PRD 73.055003 (Arganda & Herrero)...
434  + CRlE(a, x) * CRlTAU(a, x).conjugate() * (mE/mTAU) * Leptf3(a, x)
435 //...until here
436  + CLlE(a, x) * CRlTAU(a, x).conjugate() * (MChi(a)/mTAU) * Leptf4(a, x)) );
437  }
438  }
439  gslpp::complex TEA1RC = 0.0;
440  gslpp::complex TEA1LC = 0.0;
441  gslpp::complex TEARC = 0.0;
442  gslpp::complex TEALC = 0.0;
443  for (int a=0;a<2;a++) {
444  for (int x=0;x<3;x++) {
445  TEA1RC = TEA1RC + AmpTEA1RC(a, x);
446  TEA1LC = TEA1LC + AmpTEA1LC(a,x);
447  TEARC = TEARC + AmpTEARC(a,x);
448  TEALC = TEALC + AmpTEALC(a,x);
449  }
450  }
451  // write AR and AL into a vector for tau->(3)e
452  AFunctions.assign(0, TEA1RC + TEA1RN );
453  AFunctions.assign(1, TEA1LC + TEA1LN );
454  AFunctions.assign(2, TEARC + TEARN );
455  AFunctions.assign(3, TEALC + TEALN );
456 
457  }
458 
459  return(AFunctions);
460 }
gslpp::matrix< gslpp::complex > AmpTauA1LN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > myRl
Definition: SUSYMatching.h:190
gslpp::matrix< gslpp::complex > AmpTauALN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpTEARN
Definition: SUSYMatching.h:225
gslpp::vector< double > mym_se_sq
Definition: SUSYMatching.h:185
gslpp::matrix< gslpp::complex > NLlE
Definition: SUSYMatching.h:220
gslpp::vector< double > MChi
Chargino mass-eigenvalue.
Definition: SUSYMatching.h:195
gslpp::vector< double > mym_sn_sq
Definition: SUSYMatching.h:185
Particle getLeptons(const StandardModel::lepton p) const
A get method to retrieve the member object of a lepton.
virtual double Mw_tree() const
The tree-level mass of the boson, .
gslpp::matrix< double > Lepty
Functions needed to calculate various LFV observables.
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > ON
Definition: SUSYMatching.h:209
gslpp::matrix< double > Leptf4
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > CRlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< double > Leptfa1
Definition: SUSYMatching.h:214
complex pow(const complex &z1, const complex &z2)
gslpp::matrix< double > Leptf2
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > AmpTEA1LC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CLlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpA1LC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpTauALC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > NRlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpA1LN
Amplitudes of Chargino and Neutralino contribution to various LFV observables.
Definition: SUSYMatching.h:225
gslpp::vector< double > MNeig
Definition: SUSYMatching.h:205
gslpp::matrix< gslpp::complex > AmpTauA1RN
Definition: SUSYMatching.h:225
gslpp::matrix< double > Leptz
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > AmpTEALN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpARN
Definition: SUSYMatching.h:225
void assign(const size_t &i, const size_t &j, const double &a)
gslpp::matrix< gslpp::complex > AmpA1RN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CRlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > myRn
Definition: SUSYMatching.h:190
gslpp::matrix< gslpp::complex > AmpTEARC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > NRlE
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpTEA1RC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > myV
Chargino mixing matrices.
Definition: SUSYMatching.h:200
gslpp::matrix< double > Leptf3
Definition: SUSYMatching.h:214
gslpp::matrix< gslpp::complex > AmpTauARN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpARC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CRlE
Chargino and Neutralino couplings to sfermions.
Definition: SUSYMatching.h:220
void NeutralinoRemixing()
gslpp::matrix< gslpp::complex > AmpTEA1RN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CLlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NLlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NRlTAU
Definition: SUSYMatching.h:220
double getSinb() const
Gets .
Definition: SUSY.h:160
gslpp::matrix< double > Leptf1
Definition: SUSYMatching.h:214
gslpp::matrix< double > Leptfa2
Definition: SUSYMatching.h:214
virtual gslpp::vector< gslpp::complex > AFunctions(int n)
Calculates gamma penguin amplitudes for the process from .
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
gslpp::matrix< gslpp::complex > AmpTEALC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpA1RC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpTauA1RC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > CLlE
Definition: SUSYMatching.h:220
complex log(const complex &z)
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
gslpp::matrix< gslpp::complex > NLlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > AmpALN
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpALC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > AmpTauARC
Definition: SUSYMatching.h:225
gslpp::matrix< gslpp::complex > myU
Definition: SUSYMatching.h:200
gslpp::matrix< gslpp::complex > AmpTauA1LC
Definition: SUSYMatching.h:225
const SUSY & mySUSY
Definition: SUSYMatching.h:177
complex sqrt(const complex &z)
gslpp::matrix< gslpp::complex > AmpTEA1LN
Definition: SUSYMatching.h:225
gslpp::vector< gslpp::complex > SUSYMatching::BFunctions ( int  n)
virtual

Calculates box diagram amplitudes for the process \( \ell_j \to \ell_i \ell_i \ell_i \) from from [90] [14].

Calculates box amplitudes for m->3e (1), t->3m (2) and t->3e (3)

Parameters
[in]ndetermines the process, e.g., 1 = \( \mu \to eee \), 2 = \( \tau \to \mu \mu \mu \), 3 = \( \tau \to eee \)
Returns
returns the vector of gamma diagram amplitude

Definition at line 462 of file SUSYMatching.cpp.

462  {
463  //box diagram contributions
464 
466 
467  double MW = mySUSY.Mw_tree();
468  double pi = M_PI;
469  double sw2 = mySUSY.StandardModel::sW2(MW);
470  double stw = sqrt(sw2);
471  double ctw = sqrt(1.0 - sw2);
472  double ttw = stw/ctw;
474  double mMU = mySUSY.getLeptons(StandardModel::MU).getMass();
475  double mTAU = mySUSY.getLeptons(StandardModel::TAU).getMass();
476 
477  double cdenc = sqrt(2.0)*MW*cosb;
478  double cdenn = MW*cosb;
479  double g2 = gW;
480  double g2t = g2/sqrt(2.0);
481  double alph = mySUSY.getAle();
482 
484 
485  // Neutralino-Fermion-Sfermion couplings
486  for (int a=0;a<4;a++) {
487  for (int x=0;x<6;x++) {
488  // LL + RL TYPE MI
489  NRlE.assign(a, x, - (g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 0) + (mE/cdenn)*ON(a, 2)*myRl(x, 3)));
490  NRlMU.assign(a, x, -(g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 1) + (mMU/cdenn)*ON(a, 2)*myRl(x, 4)));
491  NRlTAU.assign(a, x, -(g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 2) + (mTAU/cdenn)*ON(a, 2)*myRl(x, 5)));
492  // RL + RR TYPE MI
493  NLlE.assign(a, x, -(g2t)*((mE/cdenn)*ON(a, 2)*myRl(x, 0) + 2.0*ON(a, 0)*ttw*myRl(x, 3)));
494  NLlMU.assign(a, x, -(g2t)*((mMU/cdenn)*ON(a, 2)*myRl(x, 1) + 2.0*ON(a, 0)*ttw*myRl(x, 4)));
495  NLlTAU.assign(a, x, -(g2t)*((mTAU/cdenn)*ON(a, 2)*myRl(x, 2) + 2.0*ON(a, 0)*ttw*myRl(x, 5)));
496 // Commented expressions might be useful for complex neutralino mixing matrices
497 // NLlE.assign(a, x, -(g2t)*((mE/cdenn)*ON(a, 2).conjugate()*myRl(x, 0) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 3)));
498 // NLlMU.assign(a, x, -(g2t)*((mMU/cdenn)*ON(a, 2).conjugate()*myRl(x, 1) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 4)));
499 // NLlTAU.assign(a, x, -(g2t)*((mTAU/cdenn)*ON(a, 2).conjugate()*myRl(x, 2) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 5)));
500  }
501  }
502 
503  // Chargino-Fermion-Sfermion couplings
504  for (int a=0;a<2;a++) {
505  for (int x=0;x<3;x++) {
506  // LL-TYPE
507  CRlE.assign(a, x, - (g2*myV(a, 0)*myRn(x, 0)));
508  CRlMU.assign(a, x, - (g2*myV(a, 0)*myRn(x, 1)));
509  CRlTAU.assign(a, x, - (g2*myV(a, 0)*myRn(x, 2)));
510  // LR-TYPE
511  CLlE.assign(a, x, g2*mE/cdenc*myU(a, 1).conjugate()*myRn(x, 0));
512  CLlMU.assign(a, x, g2*mMU/cdenc*myU(a, 1).conjugate()*myRn(x, 1));
513  CLlTAU.assign(a, x, g2*mTAU/cdenc*myU(a, 1).conjugate()*myRn(x, 2));
514  }
515  }
516 
517  if (li_to_lj == 1) // mu -> 3e
518  {
519  // Neutralino contributions
520 // J4n(a,b,x,t)=/*J4n(a,b,x,t)*/
521 // *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
522 // I4n(a,b,x,t)=/*I4n(a,b,x,t)*/
523 // *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi)
524  gslpp::complex B1nRMu3E = 0.0;
525  gslpp::complex B2nRMu3E = 0.0;
526  gslpp::complex B3nRMu3E = 0.0;
527  gslpp::complex B4nRMu3E = 0.0;
528  gslpp::complex B1nLMu3E = 0.0;
529  gslpp::complex B2nLMu3E = 0.0;
530  gslpp::complex B3nLMu3E = 0.0;
531  gslpp::complex B4nLMu3E = 0.0;
532  for (int a=0;a<4;a++) {
533  for (int b=0;b<4;b++) {
534  for (int x=0;x<6;x++) {
535  for (int t=0;t<6;t++) {
536  B1nRMu3E = B1nRMu3E + (1.0/(4.0*pi*alph))*(0.5*NLlMU(a,x)*NLlE(a,t)*NLlE(b,t)*NLlE(b,x)/*J4n(a,b,x,t)*/
537  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
538  +MNeig(a)*MNeig(b)*NLlMU(a,x)*NLlE(a,t)*NLlE(b,t)*NLlE(b,x)/*I4n(a,b,x,t)*/
539  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
540  B2nRMu3E = B2nRMu3E + (1.0/(4.0*pi*alph))*(0.25*(NLlMU(a,x)*NLlE(a,t)*NRlE(b,t)*NRlE(b,x)
541  +NLlMU(a,x)*NRlE(a,t)*NLlE(b,t)*NRlE(b,x)
542  -NLlMU(a,x)*NRlE(a,t)*NRlE(b,t)*NLlE(b,x))/*J4n(a,b,x,t)*/
543  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
544  -0.5*MNeig(a)*MNeig(b)*NLlMU(a,x)*NRlE(a,t)*NRlE(b,t)*NLlE(b,x)/*I4n(a,b,x,t)*/
545  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
546  B3nRMu3E = B3nRMu3E + (1.0/(4.0*pi*alph))*MNeig(a)*MNeig(b)*(NLlMU(a,x)*NRlE(a,t)*NLlE(b,t)*NRlE(b,x)
547  +0.5*NLlMU(a,x)*NLlE(a,t)*NRlE(b,t)*NRlE(b,x))/*I4n(a,b,x,t)*/
548  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
549  B4nRMu3E = B4nRMu3E + (1.0/(32.0*pi*alph))*MNeig(a)*MNeig(b)*NLlMU(a,x)*NLlE(a,t)*NRlE(b,t)*NRlE(b,x)/*I4n(a,b,x,t)*/
550  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
551  B1nLMu3E = B1nLMu3E + (1.0/(4.0*pi*alph))*(0.5*NRlMU(a,x)*NRlE(a,t)*NRlE(b,t)*NRlE(b,x)/*J4n(a,b,x,t)*/
552  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
553  +MNeig(a)*MNeig(b)*NRlMU(a,x)*NRlE(a,t)*NRlE(b,t)*NRlE(b,x)/*I4n(a,b,x,t)*/
554  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
555  B2nLMu3E = B2nLMu3E + (1.0/(4.0*pi*alph))*(0.25*(NRlMU(a,x)*NRlE(a,t)*NLlE(b,t)*NLlE(b,x)
556  +NRlMU(a,x)*NLlE(a,t)*NRlE(b,t)*NLlE(b,x)
557  -NRlMU(a,x)*NLlE(a,t)*NLlE(b,t)*NRlE(b,x))/*J4n(a,b,x,t)*/
558  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
559  -0.5*MNeig(a)*MNeig(b)*NRlMU(a,x)*NLlE(a,t)*NLlE(b,t)*NRlE(b,x)/*I4n(a,b,x,t)*/
560  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
561  B3nLMu3E = B3nLMu3E + (1.0/(4.0*pi*alph))*MNeig(a)*MNeig(b)*(NRlMU(a,x)*NLlE(a,t)*NRlE(b,t)*NLlE(b,x)
562  +0.5*NRlMU(a,x)*NRlE(a,t)*NLlE(b,t)*NLlE(b,x))/*I4n(a,b,x,t)*/
563  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
564  B4nLMu3E = B4nLMu3E + (1.0/(32.0*pi*alph))*MNeig(a)*MNeig(b)*NRlMU(a,x)*NRlE(a,t)*NLlE(b,t)*NLlE(b,x)/*I4n(a,b,x,t)*/
565  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
566  }
567  }
568  }
569  }
570  // Chargino contributions
571 // J4c(a,b,x,t)=/*J4c(a,b,x,t)*/
572 // *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi)
573 // I4c(a,b,x,t)=/*I4c(a,b,x,t)*/
574 // *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi)
575  gslpp::complex B1cRMu3E = 0.0;
576  gslpp::complex B2cRMu3E = 0.0;
577  gslpp::complex B3cRMu3E = 0.0;
578  gslpp::complex B1cLMu3E = 0.0;
579  gslpp::complex B2cLMu3E = 0.0;
580  gslpp::complex B3cLMu3E = 0.0;
581  for (int a=0;a<2;a++) {
582  for (int b=0;b<2;b++) {
583  for (int x=0;x<3;x++) {
584  for (int t=0;t<3;t++) {
585  B1cRMu3E = B1cRMu3E + (1.0/(8.0*pi*alph))*CLlMU(a,x)*CLlE(a,t)*CLlE(b,t)*CLlE(b,x)/*J4c(a,b,x,t)*/
586  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi);
587  B2cRMu3E = B2cRMu3E + (1.0/(4.0*pi*alph))*(0.25*CLlMU(a,x)*CLlE(a,t)*CRlE(b,t)*CRlE(b,x)/*J4c(a,b,x,t)*/
588  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi)
589  -0.5*MChi(a)*MChi(b)*CLlMU(a,x)*CRlE(a,t)*CRlE(b,t)*CLlE(b,x)/*I4c(a,b,x,t)*/
590  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi));
591  B3cRMu3E = B3cRMu3E + (1.0/(4.0*pi*alph))*MChi(a)*MChi(b)*CLlMU(a,x)*CRlE(a,t)*CLlE(b,t)*CRlE(b,x)/*I4c(a,b,x,t)*/
592  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi);
593  B1cLMu3E = B1cLMu3E + (1.0/(8.0*pi*alph))*(CRlMU(a,x)*CRlE(a,t)*CRlE(b,t)*CRlE(b,x)/*J4c(a,b,x,t)*/
594  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi));
595  B2cLMu3E = B2cLMu3E + (1.0/(4.0*pi*alph))*(0.25*CRlMU(a,x)*CRlE(a,t)*CLlE(b,t)*CLlE(b,x)/*J4c(a,b,x,t)*/
596  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi)
597  -0.5*MChi(a)*MChi(b)*CRlMU(a,x)*CLlE(a,t)*CLlE(b,t)*CRlE(b,x)/*I4c(a,b,x,t)*/
598  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi));
599  B3cLMu3E = B3cLMu3E + (1.0/(4.0*pi*alph))*MChi(a)*MChi(b)*CRlMU(a,x)*CLlE(a,t)*CRlE(b,t)*CLlE(b,x)/*I4c(a,b,x,t)*/
600  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi);
601  }
602  }
603  }
604  }
605 
606 // write BR and BL into a vector for mu->3e
607  BFunctions.assign(0, B1nRMu3E + B1cRMu3E );
608  BFunctions.assign(1, B1nLMu3E + B1cLMu3E );
609  BFunctions.assign(2, B2nRMu3E + B2cRMu3E );
610  BFunctions.assign(3, B2nLMu3E + B2cLMu3E );
611  BFunctions.assign(4, B3nRMu3E + B3cRMu3E );
612  BFunctions.assign(5, B3nLMu3E + B3cLMu3E );
613  BFunctions.assign(6, B4nRMu3E );
614  BFunctions.assign(7, B4nLMu3E );
615  }
616  if (li_to_lj == 2) // tau -> 3mu
617  {
618  // Neutralino contributions
619  gslpp::complex B1nRTau3Mu = 0.0;
620  gslpp::complex B2nRTau3Mu = 0.0;
621  gslpp::complex B3nRTau3Mu = 0.0;
622  gslpp::complex B4nRTau3Mu = 0.0;
623  gslpp::complex B1nLTau3Mu = 0.0;
624  gslpp::complex B2nLTau3Mu = 0.0;
625  gslpp::complex B3nLTau3Mu = 0.0;
626  gslpp::complex B4nLTau3Mu = 0.0;
627  for (int a=0;a<4;a++) {
628  for (int b=0;b<4;b++) {
629  for (int x=0;x<6;x++) {
630  for (int t=0;t<6;t++) {
631  B1nRTau3Mu = B1nRTau3Mu + (1.0/(4.0*pi*alph))*(0.5*NLlTAU(a,x)*NLlMU(a,t)*NLlMU(b,t)*NLlMU(b,x)/*J4n(a,b,x,t)*/
632  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
633  +MNeig(a)*MNeig(b)*NLlTAU(a,x)*NLlMU(a,t)*NLlMU(b,t)*NLlMU(b,x)/*I4n(a,b,x,t)*/
634  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
635  B2nRTau3Mu = B2nRTau3Mu + (1.0/(4.0*pi*alph))*(0.25*(NLlTAU(a,x)*NLlMU(a,t)*NRlMU(b,t)*NRlMU(b,x)
636  +NLlTAU(a,x)*NRlMU(a,t)*NLlMU(b,t)*NRlMU(b,x)
637  -NLlTAU(a,x)*NRlMU(a,t)*NRlMU(b,t)*NLlMU(b,x))/*J4n(a,b,x,t)*/
638  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
639  -0.5*MNeig(a)*MNeig(b)*NLlTAU(a,x)*NRlMU(a,t)*NRlMU(b,t)*NLlMU(b,x)/*I4n(a,b,x,t)*/
640  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
641  B3nRTau3Mu = B3nRTau3Mu + (1.0/(4.0*pi*alph))*MNeig(a)*MNeig(b)*(NLlTAU(a,x)*NRlMU(a,t)*NLlMU(b,t)*NRlMU(b,x)
642  +0.5*NLlTAU(a,x)*NLlMU(a,t)*NRlMU(b,t)*NRlMU(b,x))/*I4n(a,b,x,t)*/
643  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
644  B4nRTau3Mu = B4nRTau3Mu + (1.0/(32.0*pi*alph))*MNeig(a)*MNeig(b)*NLlTAU(a,x)*NLlMU(a,t)*NRlMU(b,t)*NRlMU(b,x)/*I4n(a,b,x,t)*/
645  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
646  B1nLTau3Mu = B1nLTau3Mu + (1.0/(4.0*pi*alph))*(0.5*NRlTAU(a,x)*NRlMU(a,t)*NRlMU(b,t)*NRlMU(b,x)/*J4n(a,b,x,t)*/
647  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
648  +MNeig(a)*MNeig(b)*NRlTAU(a,x)*NRlMU(a,t)*NRlMU(b,t)*NRlMU(b,x)/*I4n(a,b,x,t)*/
649  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
650  B2nLTau3Mu = B2nLTau3Mu + (1.0/(4.0*pi*alph))*(0.25*(NRlTAU(a,x)*NRlMU(a,t)*NLlMU(b,t)*NLlMU(b,x)
651  +NRlTAU(a,x)*NLlMU(a,t)*NRlMU(b,t)*NLlMU(b,x)
652  -NRlTAU(a,x)*NLlMU(a,t)*NLlMU(b,t)*NRlMU(b,x))/*J4n(a,b,x,t)*/
653  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
654  -0.5*MNeig(a)*MNeig(b)*NRlTAU(a,x)*NLlMU(a,t)*NLlMU(b,t)*NRlMU(b,x)/*I4n(a,b,x,t)*/
655  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
656  B3nLTau3Mu = B3nLTau3Mu + (1.0/(4.0*pi*alph))*MNeig(a)*MNeig(b)*(NRlTAU(a,x)*NLlMU(a,t)*NRlMU(b,t)*NLlMU(b,x)
657  +0.5*NRlTAU(a,x)*NRlMU(a,t)*NLlMU(b,t)*NLlMU(b,x))/*I4n(a,b,x,t)*/
658  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
659  B4nLTau3Mu = B4nLTau3Mu + (1.0/(32.0*pi*alph))*MNeig(a)*MNeig(b)*NRlTAU(a,x)*NRlMU(a,t)*NLlMU(b,t)*NLlMU(b,x)/*I4n(a,b,x,t)*/
660  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
661  }
662  }
663  }
664  }
665  // Chargino contributions
666  gslpp::complex B1cRTau3Mu = 0.0;
667  gslpp::complex B2cRTau3Mu = 0.0;
668  gslpp::complex B3cRTau3Mu = 0.0;
669  gslpp::complex B1cLTau3Mu = 0.0;
670  gslpp::complex B2cLTau3Mu = 0.0;
671  gslpp::complex B3cLTau3Mu = 0.0;
672  for (int a=0;a<2;a++) {
673  for (int b=0;b<2;b++) {
674  for (int x=0;x<3;x++) {
675  for (int t=0;t<3;t++) {
676  B1cRTau3Mu = B1cRTau3Mu + (1.0/(8.0*pi*alph))*CLlTAU(a,x)*CLlMU(a,t)*CLlMU(b,t)*CLlMU(b,x)/*J4c(a,b,x,t)*/
677  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi);
678  B2cRTau3Mu = B2cRTau3Mu + (1.0/(4.0*pi*alph))*(0.25*CLlTAU(a,x)*CLlMU(a,t)*CRlMU(b,t)*CRlMU(b,x)/*J4c(a,b,x,t)*/
679  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi)
680  -0.5*MChi(a)*MChi(b)*CLlTAU(a,x)*CRlMU(a,t)*CRlMU(b,t)*CLlMU(b,x)/*I4c(a,b,x,t)*/
681  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi));
682  B3cRTau3Mu = B3cRTau3Mu + (1.0/(4.0*pi*alph))*MChi(a)*MChi(b)*CLlTAU(a,x)*CRlMU(a,t)*CLlMU(b,t)*CRlMU(b,x)/*I4c(a,b,x,t)*/
683  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi);
684  B1cLTau3Mu = B1cLTau3Mu + (1.0/(8.0*pi*alph))*CRlTAU(a,x)*CRlMU(a,t)*CRlMU(b,t)*CRlMU(b,x)/*J4c(a,b,x,t)*/
685  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi);
686  B2cLTau3Mu = B2cLTau3Mu + (1.0/(4.0*pi*alph))*(0.25*CRlTAU(a,x)*CRlMU(a,t)*CLlMU(b,t)*CLlMU(b,x)/*J4c(a,b,x,t)*/
687  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi)
688  -0.5*MChi(a)*MChi(b)*CRlTAU(a,x)*CLlMU(a,t)*CLlMU(b,t)*CRlMU(b,x)/*I4c(a,b,x,t)*/
689  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi));
690  B3cLTau3Mu = B3cLTau3Mu + (1.0/(4.0*pi*alph))*MChi(a)*MChi(b)*CRlTAU(a,x)*CLlMU(a,t)*CRlMU(b,t)*CLlMU(b,x)/*I4c(a,b,x,t)*/
691  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi);
692  }
693  }
694  }
695  }
696  // write BR and BL into a vector for tau->3mu
697  BFunctions.assign(0, B1nRTau3Mu + B1cRTau3Mu );
698  BFunctions.assign(1, B1nLTau3Mu + B1cLTau3Mu );
699  BFunctions.assign(2, B2nRTau3Mu + B2cRTau3Mu );
700  BFunctions.assign(3, B2nLTau3Mu + B2cLTau3Mu );
701  BFunctions.assign(4, B3nRTau3Mu + B3cRTau3Mu );
702  BFunctions.assign(5, B3nLTau3Mu + B3cLTau3Mu );
703  BFunctions.assign(6, B4nRTau3Mu );
704  BFunctions.assign(7, B4nLTau3Mu );
705  }
706  if (li_to_lj == 3) // tau -> 3e
707  {
708  // Neutralino contributions
709  gslpp::complex B1nRTau3E = 0.0;
710  gslpp::complex B2nRTau3E = 0.0;
711  gslpp::complex B3nRTau3E = 0.0;
712  gslpp::complex B4nRTau3E = 0.0;
713  gslpp::complex B1nLTau3E = 0.0;
714  gslpp::complex B2nLTau3E = 0.0;
715  gslpp::complex B3nLTau3E = 0.0;
716  gslpp::complex B4nLTau3E = 0.0;
717  for (int a=0;a<4;a++) {
718  for (int b=0;b<4;b++) {
719  for (int x=0;x<6;x++) {
720  for (int t=0;t<6;t++) {
721  B1nRTau3E = B1nRTau3E + (1.0/(4.0*pi*alph))*(0.5*NLlTAU(a,x)*NLlE(a,t)*NLlE(b,t)*NLlE(b,x)/*J4n(a,b,x,t)*/
722  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
723  +MNeig(a)*MNeig(b)*NLlTAU(a,x)*NLlE(a,t)*NLlE(b,t)*NLlE(b,x)/*I4n(a,b,x,t)*/
724  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
725  B2nRTau3E = B2nRTau3E + (1.0/(4.0*pi*alph))*(0.25*(NLlTAU(a,x)*NLlE(a,t)*NRlE(b,t)*NRlE(b,x)
726  +NLlTAU(a,x)*NRlE(a,t)*NLlE(b,t)*NRlE(b,x)
727  -NLlTAU(a,x)*NRlE(a,t)*NRlE(b,t)*NLlE(b,x))/*J4n(a,b,x,t)*/
728  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
729  -0.5*MNeig(a)*MNeig(b)*NLlTAU(a,x)*NRlE(a,t)*NRlE(b,t)*NLlE(b,x)/*I4n(a,b,x,t)*/
730  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
731  B3nRTau3E = B3nRTau3E + (1.0/(4.0*pi*alph))*MNeig(a)*MNeig(b)*(NLlTAU(a,x)*NRlE(a,t)*NLlE(b,t)*NRlE(b,x)
732  +0.5*NLlTAU(a,x)*NLlE(a,t)*NRlE(b,t)*NRlE(b,x))/*I4n(a,b,x,t)*/
733  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
734  B4nRTau3E = B4nRTau3E + (1.0/(32.0*pi*alph))*MNeig(a)*MNeig(b)*NLlTAU(a,x)*NLlE(a,t)*NRlE(b,t)*NRlE(b,x)/*I4n(a,b,x,t)*/
735  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
736  B1nLTau3E = B1nLTau3E + (1.0/(4.0*pi*alph))*(0.5*NRlTAU(a,x)*NRlE(a,t)*NRlE(b,t)*NRlE(b,x)/*J4n(a,b,x,t)*/
737  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
738  +MNeig(a)*MNeig(b)*NRlTAU(a,x)*NRlE(a,t)*NRlE(b,t)*NRlE(b,x)/*I4n(a,b,x,t)*/
739  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
740  B2nLTau3E = B2nLTau3E + (1.0/(4.0*pi*alph))*(0.25*(NRlTAU(a,x)*NRlE(a,t)*NLlE(b,t)*NLlE(b,x)
741  +NRlTAU(a,x)*NLlE(a,t)*NRlE(b,t)*NLlE(b,x)
742  -NRlTAU(a,x)*NLlE(a,t)*NLlE(b,t)*NRlE(b,x))/*J4n(a,b,x,t)*/
743  *PV.D00(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(4.0*pi*pi)
744  -0.5*MNeig(a)*MNeig(b)*NRlTAU(a,x)*NLlE(a,t)*NLlE(b,t)*NRlE(b,x)/*I4n(a,b,x,t)*/
745  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi));
746  B3nLTau3E = B3nLTau3E + (1.0/(4.0*pi*alph))*MNeig(a)*MNeig(b)*(NRlTAU(a,x)*NLlE(a,t)*NRlE(b,t)*NLlE(b,x)
747  +0.5*NRlTAU(a,x)*NRlE(a,t)*NLlE(b,t)*NLlE(b,x))/*I4n(a,b,x,t)*/
748  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
749  B4nLTau3E = B4nLTau3E + (1.0/(32.0*pi*alph))*MNeig(a)*MNeig(b)*NRlTAU(a,x)*NRlE(a,t)*NLlE(b,t)*NLlE(b,x)/*I4n(a,b,x,t)*/
750  *PV.D0(0., 0., MNeig(a)*MNeig(a), MNeig(b)*MNeig(b), mym_se_sq(x), mym_se_sq(t))/(16.0*pi*pi);
751  }
752  }
753  }
754  }
755  // Chargino contributions
756  gslpp::complex B1cRTau3E = 0.0;
757  gslpp::complex B2cRTau3E = 0.0;
758  gslpp::complex B3cRTau3E = 0.0;
759  gslpp::complex B1cLTau3E = 0.0;
760  gslpp::complex B2cLTau3E = 0.0;
761  gslpp::complex B3cLTau3E = 0.0;
762  for (int a=0;a<2;a++) {
763  for (int b=0;b<2;b++) {
764  for (int x=0;x<3;x++) {
765  for (int t=0;t<3;t++) {
766  B1cRTau3E = B1cRTau3E + (1.0/(8.0*pi*alph))*CLlTAU(a,x)*CLlE(a,t)*CLlE(b,t)*CLlE(b,x)/*J4c(a,b,x,t)*/
767  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi);
768  B2cRTau3E = B2cRTau3E + (1.0/(4.0*pi*alph))*(0.25*CLlTAU(a,x)*CLlE(a,t)*CRlE(b,t)*CRlE(b,x)/*J4c(a,b,x,t)*/
769  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi)
770  -0.5*MChi(a)*MChi(b)*CLlTAU(a,x)*CRlE(a,t)*CRlE(b,t)*CLlE(b,x)/*I4c(a,b,x,t)*/
771  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi));
772  B3cRTau3E = B3cRTau3E + (1.0/(4.0*pi*alph))*MChi(a)*MChi(b)*CLlTAU(a,x)*CRlE(a,t)*CLlE(b,t)*CRlE(b,x)/*I4c(a,b,x,t)*/
773  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi);
774  B1cLTau3E = B1cLTau3E + (1.0/(8.0*pi*alph))*CRlTAU(a,x)*CRlE(a,t)*CRlE(b,t)*CRlE(b,x)/*J4c(a,b,x,t)*/
775  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi);
776  B2cLTau3E = B2cLTau3E + (1.0/(4.0*pi*alph))*(0.25*CRlTAU(a,x)*CRlE(a,t)*CLlE(b,t)*CLlE(b,x)/*J4c(a,b,x,t)*/
777  *PV.D00(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(4.0*pi*pi)
778  -0.5*MChi(a)*MChi(b)*CRlTAU(a,x)*CLlE(a,t)*CLlE(b,t)*CRlE(b,x)/*I4c(a,b,x,t)*/
779  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi));
780  B3cLTau3E = B3cLTau3E + (1.0/(4.0*pi*alph))*MChi(a)*MChi(b)*CRlTAU(a,x)*CLlE(a,t)*CRlE(b,t)*CLlE(b,x)/*I4c(a,b,x,t)*/
781  *PV.D0(0., 0., MChi(a)*MChi(a), MChi(b)*MChi(b), mym_sn_sq(x), mym_sn_sq(t))/(16.0*pi*pi);
782  }
783  }
784  }
785  }
786  // write BR and BL into a vector for tau->3e
787  BFunctions.assign(0, B1nRTau3E + B1cRTau3E );
788  BFunctions.assign(1, B1nLTau3E + B1cLTau3E );
789  BFunctions.assign(2, B2nRTau3E + B2cRTau3E );
790  BFunctions.assign(3, B2nLTau3E + B2cLTau3E );
791  BFunctions.assign(4, B3nRTau3E + B3cRTau3E );
792  BFunctions.assign(5, B3nLTau3E + B3cLTau3E );
793  BFunctions.assign(6, B4nRTau3E );
794  BFunctions.assign(7, B4nLTau3E );
795  }
796 // if (li_to_lj == 4) // tau -> muee
797 // {
798 // // write BR and BL into a vector for tau->muee
799 // BFunctions.assign(0, 0. );
800 // BFunctions.assign(1, 0. );
801 // BFunctions.assign(2, 0. );
802 // BFunctions.assign(3, 0. );
803 // BFunctions.assign(4, 0. );
804 // BFunctions.assign(5, 0. );
805 // BFunctions.assign(6, 0. );
806 // BFunctions.assign(7, 0. );
807 // }
808  std::cout<<"BFunctions("<<li_to_lj<<") = "<<BFunctions<<std::endl;
809 
810  return(BFunctions);
811 }
gslpp::matrix< gslpp::complex > myRl
Definition: SUSYMatching.h:190
gslpp::vector< double > mym_se_sq
Definition: SUSYMatching.h:185
gslpp::matrix< gslpp::complex > NLlE
Definition: SUSYMatching.h:220
gslpp::vector< double > MChi
Chargino mass-eigenvalue.
Definition: SUSYMatching.h:195
gslpp::vector< double > mym_sn_sq
Definition: SUSYMatching.h:185
Particle getLeptons(const StandardModel::lepton p) const
A get method to retrieve the member object of a lepton.
virtual double Mw_tree() const
The tree-level mass of the boson, .
gslpp::matrix< gslpp::complex > ON
Definition: SUSYMatching.h:209
gslpp::matrix< gslpp::complex > CRlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > CLlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NRlMU
Definition: SUSYMatching.h:220
gslpp::vector< double > MNeig
Definition: SUSYMatching.h:205
virtual gslpp::vector< gslpp::complex > BFunctions(int n)
Calculates box diagram amplitudes for the process from from .
gslpp::matrix< gslpp::complex > CRlMU
Definition: SUSYMatching.h:220
const PVfunctions PV
Definition: SUSYMatching.h:178
gslpp::matrix< gslpp::complex > myRn
Definition: SUSYMatching.h:190
gslpp::matrix< gslpp::complex > NRlE
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > myV
Chargino mixing matrices.
Definition: SUSYMatching.h:200
gslpp::matrix< gslpp::complex > CRlE
Chargino and Neutralino couplings to sfermions.
Definition: SUSYMatching.h:220
void NeutralinoRemixing()
gslpp::matrix< gslpp::complex > CLlMU
Definition: SUSYMatching.h:220
gslpp::complex D0(const double s, const double t, const double m02, const double m12, const double m22, const double m32) const
.
gslpp::matrix< gslpp::complex > NLlTAU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > NRlTAU
Definition: SUSYMatching.h:220
double getAle() const
A get method to retrieve the fine-structure constant .
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
gslpp::complex D00(const double s, const double t, const double m02, const double m12, const double m22, const double m32) const
.
gslpp::matrix< gslpp::complex > CLlE
Definition: SUSYMatching.h:220
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
gslpp::matrix< gslpp::complex > NLlMU
Definition: SUSYMatching.h:220
gslpp::matrix< gslpp::complex > myU
Definition: SUSYMatching.h:200
const SUSY & mySUSY
Definition: SUSYMatching.h:177
complex sqrt(const complex &z)
gslpp::vector< gslpp::complex > SUSYMatching::BHFunctions ( int  n)
virtual

Calculates Higgs penguin amplitudes for the process \( \ell_j \to \ell_i \ell_i \ell_i \) from [14].

Calculates Higgs penguin amplitudes for m->3e (1), t->3m (2) and t->3e (3)

Parameters
[in]ndetermines the process, e.g., 1 = \( \mu \to eee \), 2 = \( \tau \to \mu \mu \mu \), 3 = \( \tau \to eee \)
Returns
returns the vector of Higgs penguin amplitude

Definition at line 818 of file SUSYMatching.cpp.

818  {
819  //Higgs penguin contributions from PhysRevD.73.055003
820 
821 
822  // To do:
823  // Are the trilinear couplings from (A24) correct?
824  // Should we update the missing parameters in the update function (sinalpha, A_l, ...)?
825 
827 
828  double MZ = mySUSY.getMz();
829  double MW = mySUSY.Mw_tree();
830  double pi = M_PI;
831  double piconst = 1.0/(32.0 * pi * pi);
832  double sw2 = mySUSY.StandardModel::sW2(MW);
833  double stw = sqrt(sw2);
834  double ctw = sqrt(1.0 - sw2);
835  double ttw = stw/ctw;
837  double mMU = mySUSY.getLeptons(StandardModel::MU).getMass();
838  double mTAU = mySUSY.getLeptons(StandardModel::TAU).getMass();
839  double cos2b = 2.0*cosb*cosb-1.0;
840  gslpp::complex sina = mySUSY.getSaeff();
841  gslpp::complex cosa = sqrt(1.0-sina*sina);
842  gslpp::complex sinapb = sina*cosb+cosa*sinb;
843  gslpp::complex cosapb = cosa*cosb-sina*sinb;
844  gslpp::complex cosbma = cosb*cosa+sinb*sina;
845  gslpp::complex sinbma = sinb*cosa-cosb*sina;
846  double mh = mySUSY.getMHl();
847  double mH = mySUSY.getMHh();
848  double mA = mySUSY.getMHa();
849  gslpp::complex M1 = mySUSY.getM1();
850  gslpp::complex M2 = mySUSY.getM2();
851  gslpp::complex muH = mySUSY.getMuH();
852  TEhat = mySUSY.getTEhat();
853 
854  double cdenc = sqrt(2.0)*MW*cosb;
855  double cdenn = MW*cosb;
856  double g2 = gW;
857  double g2t = g2/sqrt(2.0);
858  double alph = mySUSY.getAle();
859 
860  //DEBUG
888 
890 
891 
892  // Neutralino-Fermion-Sfermion couplings
893  for (int a=0;a<4;a++) {
894  for (int x=0;x<6;x++) {
895  // LL + RL TYPE MI
896  NRlE.assign(a, x, - (g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 0) + (mE/cdenn)*ON(a, 2)*myRl(x, 3)));
897  NRlMU.assign(a, x, -(g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 1) + (mMU/cdenn)*ON(a, 2)*myRl(x, 4)));
898  NRlTAU.assign(a, x, -(g2t)*((-ON(a, 1) - ON(a, 0)*ttw)*myRl(x, 2) + (mTAU/cdenn)*ON(a, 2)*myRl(x, 5)));
899  // RL + RR TYPE MI
900  NLlE.assign(a, x, -(g2t)*((mE/cdenn)*ON(a, 2)*myRl(x, 0) + 2.0*ON(a, 0)*ttw*myRl(x, 3)));
901  NLlMU.assign(a, x, -(g2t)*((mMU/cdenn)*ON(a, 2)*myRl(x, 1) + 2.0*ON(a, 0)*ttw*myRl(x, 4)));
902  NLlTAU.assign(a, x, -(g2t)*((mTAU/cdenn)*ON(a, 2)*myRl(x, 2) + 2.0*ON(a, 0)*ttw*myRl(x, 5)));
903 // Commented expressions might be useful for complex neutralino mixing matrices
904 // NLlE.assign(a, x, -(g2t)*((mE/cdenn)*ON(a, 2).conjugate()*myRl(x, 0) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 3)));
905 // NLlMU.assign(a, x, -(g2t)*((mMU/cdenn)*ON(a, 2).conjugate()*myRl(x, 1) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 4)));
906 // NLlTAU.assign(a, x, -(g2t)*((mTAU/cdenn)*ON(a, 2).conjugate()*myRl(x, 2) + 2.0*ON(a, 0).conjugate()*ttw*myRl(x, 5)));
907  }
908  }
909 
910  // Chargino-Fermion-Sfermion couplings
911  for (int a=0;a<2;a++) {
912  for (int x=0;x<3;x++) {
913  // LL-TYPE
914  CRlE.assign(a, x, - (g2*myV(a, 0)*myRn(x, 0)));
915  CRlMU.assign(a, x, - (g2*myV(a, 0)*myRn(x, 1)));
916  CRlTAU.assign(a, x, - (g2*myV(a, 0)*myRn(x, 2)));
917  // LR-TYPE
918  CLlE.assign(a, x, g2*mE/cdenc*myU(a, 1).conjugate()*myRn(x, 0));
919  CLlMU.assign(a, x, g2*mMU/cdenc*myU(a, 1).conjugate()*myRn(x, 1));
920  CLlTAU.assign(a, x, g2*mTAU/cdenc*myU(a, 1).conjugate()*myRn(x, 2));
921  }
922  }
923 
924  gslpp::vector<gslpp::complex> sigma1(3, 0.);
925  gslpp::vector<gslpp::complex> sigma2(3, 0.);
926  gslpp::vector<gslpp::complex> sigma3(3, 0.);
927  gslpp::vector<gslpp::complex> sigma4(3, 0.);
928  gslpp::vector<gslpp::complex> sigma5(3, 0.);
929  sigma1.assign(0, sina);
930  sigma1.assign(1, -cosa);
931  sigma1.assign(2, sinb*gslpp::complex::i());
932  sigma2.assign(0, cosa);
933  sigma2.assign(1, sina);
934  sigma2.assign(2, -cosb*gslpp::complex::i());
935  sigma3.assign(0, sinapb);
936  sigma3.assign(1, -cosapb);
937  sigma3.assign(2, 0.);
938  sigma4.assign(0, -sina);
939  sigma4.assign(1, cosa);
940  sigma4.assign(2, 0.);
941  sigma5.assign(0, -cosbma);
942  sigma5.assign(1, sinbma);
943  sigma5.assign(2, cos2b*gslpp::complex::i());
944 
945  gslpp::matrix<gslpp::complex> Qpp(4, 4, 0.), Rpp(4, 4, 0.);
946  gslpp::matrix<gslpp::complex> DL0(4, 4, 0.), DR0(4, 4, 0.), DL1(4, 4, 0.), DR1(4, 4, 0.), DL2(4, 4, 0.), DR2(4, 4, 0.);
947  for (int a=0;a<4;a++) {
948  for (int b=0;b<4;b++) {
949  Qpp.assign(a, b, 0.5*(ON(a,2)*(ON(b,1)-ttw*ON(b,0))+ON(b,2)*(ON(a,1)-ttw*ON(a,0))) );
950  Rpp.assign(a, b, (M2.conjugate()*ON(a,1)*ON(b,1) +M1.conjugate()*ON(a,0)*ON(b,0) -muH.conjugate()*(ON(a,2)*ON(b,3)+ON(a,3)*ON(b,2)))/(2.0*MW) );
951  DL0.assign(b, a, -g2/sinb * (Qpp(a,b).conjugate()*sigma5(0) -Rpp(a,b).conjugate()*sigma2(0) +MNeig(a)/(2.0*MW)*sigma2(0)*delta_ab(a,b)) );
952  DR0.assign(b, a, DL0(b,a).conjugate() );
953  DL1.assign(b, a, -g2/sinb * (Qpp(a,b).conjugate()*sigma5(1) -Rpp(a,b).conjugate()*sigma2(1) +MNeig(a)/(2.0*MW)*sigma2(1)*delta_ab(a,b)) );
954  DR1.assign(b, a, DL1(b,a).conjugate() );
955  DL2.assign(b, a, -g2/sinb * (Qpp(a,b).conjugate()*sigma5(2) -Rpp(a,b).conjugate()*sigma2(2) +MNeig(a)/(2.0*MW)*sigma2(2)*delta_ab(a,b)) );
956  DR2.assign(b, a, DL2(b,a).conjugate() );
957  }
958  }
959 
960  gslpp::matrix<gslpp::complex> Qch(2, 2, 0.), Rch(2, 2, 0.);
961  gslpp::matrix<gslpp::complex> WL0(2, 2, 0.), WR0(2, 2, 0.), WL1(2, 2, 0.), WR1(2, 2, 0.), WL2(2, 2, 0.), WR2(2, 2, 0.);
962  for (int a=0;a<2;a++) {
963  for (int b=0;b<2;b++) {
964  Qch.assign(a, b, myU(a,1)*myV(b,0)/sqrt(2.0) );
965  Rch.assign(a, b, (M2.conjugate()*myU(a,0)*myV(b,0) +muH.conjugate()*myU(a,1)*myV(b,1))/(2.0*MW) );
966  WR0.assign(a, b, -g2/sinb * (Qch(a,b)*sigma5(0).conjugate() -Rch(a,b)*sigma2(0).conjugate() +MChi(a)/(2.0*MW)*sigma2(0).conjugate()*delta_ab(a,b)) );
967  WL0.assign(b, a, WR0(a,b).conjugate() );
968  WR1.assign(a, b, -g2/sinb * (Qch(a,b)*sigma5(1).conjugate() -Rch(a,b)*sigma2(1).conjugate() +MChi(a)/(2.0*MW)*sigma2(1).conjugate()*delta_ab(a,b)) );
969  WL1.assign(b, a, WR1(a,b).conjugate() );
970  WR2.assign(a, b, -g2/sinb * (Qch(a,b)*sigma5(2).conjugate() -Rch(a,b)*sigma2(2).conjugate() +MChi(a)/(2.0*MW)*sigma2(2).conjugate()*delta_ab(a,b)) );
971  WL2.assign(b, a, WR2(a,b).conjugate() );
972  }
973  }
974 
975  gslpp::vector<gslpp::complex> gLLE(3, 0.), gRRE(3, 0.), gLRE(3, 0.), gRLE(3, 0.);
976  gslpp::vector<gslpp::complex> gLLMU(3, 0.), gRRMU(3, 0.), gLRMU(3, 0.), gRLMU(3, 0.);
977  gslpp::vector<gslpp::complex> gLLTAU(3, 0.), gRRTAU(3, 0.), gLRTAU(3, 0.), gRLTAU(3, 0.);
978  gslpp::vector<gslpp::complex> gLLNU(3, 0.);
979  for (int p=0;p<3;p++) {
980  gLLE.assign(p, MZ/ctw*sigma3(p)*(0.5-sw2) + mE*mE/(MW*cosb)*sigma4(p));
981  gLLMU.assign(p, MZ/ctw*sigma3(p)*(0.5-sw2) + mMU*mMU/(MW*cosb)*sigma4(p));
982  gLLTAU.assign(p, MZ/ctw*sigma3(p)*(0.5-sw2) + mTAU*mTAU/(MW*cosb)*sigma4(p));
983  gRRE.assign(p, MZ/ctw*sigma3(p)*sw2 + mE*mE/(MW*cosb)*sigma4(p));
984  gRRMU.assign(p, MZ/ctw*sigma3(p)*sw2 + mMU*mMU/(MW*cosb)*sigma4(p));
985  gRRTAU.assign(p, MZ/ctw*sigma3(p)*sw2 + mTAU*mTAU/(MW*cosb)*sigma4(p));
986  gLRE.assign(p, (-sigma1(p)*TEhat(0,0)/mE*v1/sqrt(2.0)-sigma2(p).conjugate()*muH)*mE/(2.0*MW*cosb));
987  gLRMU.assign(p, (-sigma1(p)*TEhat(1,1)/mMU*v1/sqrt(2.0)-sigma2(p).conjugate()*muH)*mMU/(2.0*MW*cosb));
988  gLRTAU.assign(p, (-sigma1(p)*TEhat(2,2)/mTAU*v1/sqrt(2.0)-sigma2(p).conjugate()*muH)*mTAU/(2.0*MW*cosb));
989  gRLE.assign(p, gLRE(p).conjugate());
990  gRLMU.assign(p, gLRMU(p).conjugate());
991  gRLTAU.assign(p, gLRTAU(p).conjugate());
992  gLLNU.assign(p, -0.5*MZ/ctw*sigma3(p));
993  }
994 
995  //Note that the dependence on the Rl and Rn elements here is different than the one in PhysRevD.73.055003
996  //due to the different choice of flavour basis.
997  //The replacement rules are:
998  //AH basis (PhysRevD.73.055003) -> SLHA basis
999  // (from 1 to 6) -> (from 0 to 5)
1000  // Rl(x,1) -> Rl(0,x)
1001  // Rl(x,2) -> Rl(3,x)
1002  // Rl(x,3) -> Rl(1,x)
1003  // Rl(x,4) -> Rl(4,x)
1004  // Rl(x,5) -> Rl(2,x)
1005  // Rl(x,6) -> Rl(5,x)
1006  //For Rn, the indices have to be swapped.
1007 
1008  gslpp::matrix<gslpp::complex> Gl0(6, 6, 0.), Gl1(6, 6, 0.), Gl2(6, 6, 0.);
1009  for (int x=0;x<6;x++) {
1010  for (int y=0;y<6;y++) {
1011  Gl0.assign(x, y, -g2*( gLLE(0)*myRl(x,0).conjugate()*myRl(y,0) +gRRE(0)*myRl(x,3).conjugate()*myRl(y,3) +gLRE(0)*myRl(x,0).conjugate()*myRl(y,3) +gRLE(0)*myRl(x,3).conjugate()*myRl(y,0)
1012  +gLLMU(0)*myRl(x,1).conjugate()*myRl(y,1) +gRRMU(0)*myRl(x,4).conjugate()*myRl(y,4) +gLRMU(0)*myRl(x,1).conjugate()*myRl(y,4) +gRLMU(0)*myRl(x,4).conjugate()*myRl(y,1)
1013  +gLLTAU(0)*myRl(x,2).conjugate()*myRl(y,2) +gRRTAU(0)*myRl(x,5).conjugate()*myRl(y,5) +gLRTAU(0)*myRl(x,2).conjugate()*myRl(y,5) +gRLTAU(0)*myRl(x,5).conjugate()*myRl(y,2)));
1014  Gl1.assign(x, y, -g2*( gLLE(1)*myRl(x,0).conjugate()*myRl(y,0) +gRRE(1)*myRl(x,3).conjugate()*myRl(y,3) +gLRE(1)*myRl(x,0).conjugate()*myRl(y,3) +gRLE(1)*myRl(x,3).conjugate()*myRl(y,0)
1015  +gLLMU(1)*myRl(x,1).conjugate()*myRl(y,1) +gRRMU(1)*myRl(x,4).conjugate()*myRl(y,4) +gLRMU(1)*myRl(x,1).conjugate()*myRl(y,4) +gRLMU(1)*myRl(x,4).conjugate()*myRl(y,1)
1016  +gLLTAU(1)*myRl(x,2).conjugate()*myRl(y,2) +gRRTAU(1)*myRl(x,5).conjugate()*myRl(y,5) +gLRTAU(1)*myRl(x,2).conjugate()*myRl(y,5) +gRLTAU(1)*myRl(x,5).conjugate()*myRl(y,2)));
1017  Gl2.assign(x, y, -g2*( gLLE(2)*myRl(x,0).conjugate()*myRl(y,0) +gRRE(2)*myRl(x,3).conjugate()*myRl(y,3) +gLRE(2)*myRl(x,0).conjugate()*myRl(y,3) +gRLE(2)*myRl(x,3).conjugate()*myRl(y,0)
1018  +gLLMU(2)*myRl(x,1).conjugate()*myRl(y,1) +gRRMU(2)*myRl(x,4).conjugate()*myRl(y,4) +gLRMU(2)*myRl(x,1).conjugate()*myRl(y,4) +gRLMU(2)*myRl(x,4).conjugate()*myRl(y,1)
1019  +gLLTAU(2)*myRl(x,2).conjugate()*myRl(y,2) +gRRTAU(2)*myRl(x,5).conjugate()*myRl(y,5) +gLRTAU(2)*myRl(x,2).conjugate()*myRl(y,5) +gRLTAU(2)*myRl(x,5).conjugate()*myRl(y,2)));
1020  }
1021  }
1022 
1023  gslpp::matrix<gslpp::complex> Gnu0(3, 3, 0.), Gnu1(3, 3, 0.), Gnu2(3, 3, 0.);
1024  for (int x=0;x<3;x++) {
1025  Gnu0.assign(x, x, -g2*gLLNU(0) );
1026  Gnu1.assign(x, x, -g2*gLLNU(1) );
1027  Gnu2.assign(x, x, -g2*gLLNU(2) );
1028  }
1029 
1030  gslpp::vector<gslpp::complex> SRE(3, 0.), SLE(3, 0.), SRMU(3, 0.), SLMU(3, 0.), SRTAU(3, 0.), SLTAU(3, 0.);
1031  for (int p=0;p<3;p++) {
1032  SRE.assign(p, g2*mE/(2.0*MW*cosb) * sigma1(p));
1033  SLE.assign(p, g2*mE/(2.0*MW*cosb) * sigma1(p).conjugate());
1034  SRMU.assign(p, g2*mMU/(2.0*MW*cosb) * sigma1(p));
1035  SLMU.assign(p, g2*mMU/(2.0*MW*cosb) * sigma1(p).conjugate());
1036  SRTAU.assign(p, g2*mTAU/(2.0*MW*cosb) * sigma1(p));
1037  SLTAU.assign(p, g2*mTAU/(2.0*MW*cosb) * sigma1(p).conjugate());
1038  }
1039 
1040  if (li_to_lj == 1) // mu -> 3e
1041  {
1042  // Neutralino contributions
1043  gslpp::complex HpengMuEEENR0 = 0.0;
1044  gslpp::complex HpengMuEEENL0 = 0.0;
1045  gslpp::complex HpengMuEEENR1 = 0.0;
1046  gslpp::complex HpengMuEEENL1 = 0.0;
1047  gslpp::complex HpengMuEEENR2 = 0.0;
1048  gslpp::complex HpengMuEEENL2 = 0.0;
1049  for (int x=0;x<6;x++) {
1050  for (int a=0;a<4;a++) {
1051  for (int b=0;b<4;b++) {
1052  // h R contribution
1053  HpengMuEEENR0 = HpengMuEEENR0 - 2.0*piconst*(NRlE(a,x)*DL0(a,b)*NLlMU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1054  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1055  +mMU*mMU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1056  +mE*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1057  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1058  +NLlE(a,x)*DR0(a,b)*NRlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1059  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1060  +NLlE(a,x)*DR0(a,b)*NLlMU(b,x).conjugate()*mE*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1061  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1062  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1063  +NRlE(a,x)*DL0(a,b)*NRlMU(b,x).conjugate()*mMU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1064  +NLlE(a,x)*DL0(a,b)*NLlMU(b,x).conjugate()*mE*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1065  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1066  +NRlE(a,x)*DR0(a,b)*NRlMU(b,x).conjugate()*mMU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1067  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1068  +NRlE(a,x)*DR0(a,b)*NLlMU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1069  // h L contribution
1070  HpengMuEEENL0 = HpengMuEEENL0 - 2.0*piconst*(NLlE(a,x)*DR0(a,b)*NRlMU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1071  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1072  +mMU*mMU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1073  +mE*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1074  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1075  +NRlE(a,x)*DL0(a,b)*NLlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1076  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1077  +NRlE(a,x)*DL0(a,b)*NRlMU(b,x).conjugate()*mE*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1078  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1079  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1080  +NLlE(a,x)*DR0(a,b)*NLlMU(b,x).conjugate()*mMU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1081  +NRlE(a,x)*DR0(a,b)*NRlMU(b,x).conjugate()*mE*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1082  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1083  +NLlE(a,x)*DL0(a,b)*NLlMU(b,x).conjugate()*mMU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1084  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1085  +NLlE(a,x)*DL0(a,b)*NRlMU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1086  // H R contribution
1087  HpengMuEEENR1 = HpengMuEEENR1 - 2.0*piconst*(NRlE(a,x)*DL1(a,b)*NLlMU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1088  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1089  +mMU*mMU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1090  +mE*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1091  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1092  +NLlE(a,x)*DR1(a,b)*NRlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1093  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1094  +NLlE(a,x)*DR1(a,b)*NLlMU(b,x).conjugate()*mE*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1095  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1096  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1097  +NRlE(a,x)*DL1(a,b)*NRlMU(b,x).conjugate()*mMU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1098  +NLlE(a,x)*DL1(a,b)*NLlMU(b,x).conjugate()*mE*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1099  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1100  +NRlE(a,x)*DR1(a,b)*NRlMU(b,x).conjugate()*mMU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1101  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1102  +NRlE(a,x)*DR1(a,b)*NLlMU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1103  // H L contribution
1104  HpengMuEEENL1 = HpengMuEEENL1 - 2.0*piconst*(NLlE(a,x)*DR1(a,b)*NRlMU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1105  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1106  +mMU*mMU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1107  +mE*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1108  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1109  +NRlE(a,x)*DL1(a,b)*NLlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1110  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1111  +NRlE(a,x)*DL1(a,b)*NRlMU(b,x).conjugate()*mE*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1112  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1113  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1114  +NLlE(a,x)*DR1(a,b)*NLlMU(b,x).conjugate()*mMU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1115  +NRlE(a,x)*DR1(a,b)*NRlMU(b,x).conjugate()*mE*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1116  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1117  +NLlE(a,x)*DL1(a,b)*NLlMU(b,x).conjugate()*mMU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1118  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1119  +NLlE(a,x)*DL1(a,b)*NRlMU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1120  // A R contribution
1121  HpengMuEEENR2 = HpengMuEEENR2 - 2.0*piconst*(NRlE(a,x)*DL2(a,b)*NLlMU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1122  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1123  +mMU*mMU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1124  +mE*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1125  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1126  +NLlE(a,x)*DR2(a,b)*NRlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1127  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1128  +NLlE(a,x)*DR2(a,b)*NLlMU(b,x).conjugate()*mE*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1129  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1130  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1131  +NRlE(a,x)*DL2(a,b)*NRlMU(b,x).conjugate()*mMU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1132  +NLlE(a,x)*DL2(a,b)*NLlMU(b,x).conjugate()*mE*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1133  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1134  +NRlE(a,x)*DR2(a,b)*NRlMU(b,x).conjugate()*mMU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1135  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1136  +NRlE(a,x)*DR2(a,b)*NLlMU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1137  // A L contribution
1138  HpengMuEEENL2 = HpengMuEEENL2 - 2.0*piconst*(NLlE(a,x)*DR2(a,b)*NRlMU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1139  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1140  +mMU*mMU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1141  +mE*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1142  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1143  +NRlE(a,x)*DL2(a,b)*NLlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1144  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1145  +NRlE(a,x)*DL2(a,b)*NRlMU(b,x).conjugate()*mE*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1146  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1147  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1148  +NLlE(a,x)*DR2(a,b)*NLlMU(b,x).conjugate()*mMU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1149  +NRlE(a,x)*DR2(a,b)*NRlMU(b,x).conjugate()*mE*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1150  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1151  +NLlE(a,x)*DL2(a,b)*NLlMU(b,x).conjugate()*mMU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1152  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1153  +NLlE(a,x)*DL2(a,b)*NRlMU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1154  }
1155  for (int y=0;y<6;y++) {
1156  // h R contribution
1157  HpengMuEEENR0 = HpengMuEEENR0 - 2.0*piconst*Gl0(x,y)*(-NLlE(a,x)*NLlMU(a,y).conjugate()*mE*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1158  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1159  -NRlE(a,x)*NRlMU(a,y).conjugate()*mMU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1160  -NRlE(a,x)*NLlMU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1161  // h L contribution
1162  HpengMuEEENL0 = HpengMuEEENL0 - 2.0*piconst*Gl0(x,y)*(-NRlE(a,x)*NRlMU(a,y).conjugate()*mE*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1163  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1164  -NLlE(a,x)*NLlMU(a,y).conjugate()*mMU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1165  -NLlE(a,x)*NRlMU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1166  // H R contribution
1167  HpengMuEEENR1 = HpengMuEEENR1 - 2.0*piconst*Gl1(x,y)*(-NLlE(a,x)*NLlMU(a,y).conjugate()*mE*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1168  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1169  -NRlE(a,x)*NRlMU(a,y).conjugate()*mMU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1170  -NRlE(a,x)*NLlMU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1171  // H L contribution
1172  HpengMuEEENL1 = HpengMuEEENL1 - 2.0*piconst*Gl1(x,y)*(-NRlE(a,x)*NRlMU(a,y).conjugate()*mE*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1173  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1174  -NLlE(a,x)*NLlMU(a,y).conjugate()*mMU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1175  -NLlE(a,x)*NRlMU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1176  // A R contribution
1177  HpengMuEEENR2 = HpengMuEEENR2 - 2.0*piconst*Gl2(x,y)*(-NLlE(a,x)*NLlMU(a,y).conjugate()*mE*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1178  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1179  -NRlE(a,x)*NRlMU(a,y).conjugate()*mMU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1180  -NRlE(a,x)*NLlMU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1181  // A L contribution
1182  HpengMuEEENL2 = HpengMuEEENL2 - 2.0*piconst*Gl2(x,y)*(-NRlE(a,x)*NRlMU(a,y).conjugate()*mE*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1183  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1184  -NLlE(a,x)*NLlMU(a,y).conjugate()*mMU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1185  -NLlE(a,x)*NRlMU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1186  }
1187  // h R contribution
1188  HpengMuEEENR0 = HpengMuEEENR0 - 2.0*piconst*(SRMU(0)/(mE*mE-mMU*mMU)*(-NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1189  +NLlE(a,x)*NRlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1190  -NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1191  +NRlE(a,x)*NLlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x)))
1192  +SRE(0)/(mMU*mMU-mE*mE)*(-NLlE(a,x)*NLlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1193  +NLlE(a,x)*NRlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1194  -NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1195  +NRlE(a,x)*NLlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))));
1196  // h L contribution
1197  HpengMuEEENL0 = HpengMuEEENL0 - 2.0*piconst*(SLMU(0)/(mE*mE-mMU*mMU)*(-NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1198  +NRlE(a,x)*NLlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1199  -NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1200  +NLlE(a,x)*NRlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x)))
1201  +SLE(0)/(mMU*mMU-mE*mE)*(-NRlE(a,x)*NRlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1202  +NRlE(a,x)*NLlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1203  -NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1204  +NLlE(a,x)*NRlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))));
1205  // H R contribution
1206  HpengMuEEENR1 = HpengMuEEENR1 - 2.0*piconst*(SRMU(1)/(mE*mE-mMU*mMU)*(-NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1207  +NLlE(a,x)*NRlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1208  -NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1209  +NRlE(a,x)*NLlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x)))
1210  +SRE(1)/(mMU*mMU-mE*mE)*(-NLlE(a,x)*NLlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1211  +NLlE(a,x)*NRlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1212  -NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1213  +NRlE(a,x)*NLlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))));
1214  // H L contribution
1215  HpengMuEEENL1 = HpengMuEEENL1 - 2.0*piconst*(SLMU(1)/(mE*mE-mMU*mMU)*(-NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1216  +NRlE(a,x)*NLlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1217  -NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1218  +NLlE(a,x)*NRlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x)))
1219  +SLE(1)/(mMU*mMU-mE*mE)*(-NRlE(a,x)*NRlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1220  +NRlE(a,x)*NLlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1221  -NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1222  +NLlE(a,x)*NRlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))));
1223  // A R contribution
1224  HpengMuEEENR2 = HpengMuEEENR2 - 2.0*piconst*(SRMU(2)/(mE*mE-mMU*mMU)*(-NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1225  +NLlE(a,x)*NRlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1226  -NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1227  +NRlE(a,x)*NLlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x)))
1228  +SRE(2)/(mMU*mMU-mE*mE)*(-NLlE(a,x)*NLlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1229  +NLlE(a,x)*NRlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1230  -NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1231  +NRlE(a,x)*NLlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))));
1232  // A L contribution
1233  HpengMuEEENL2 = HpengMuEEENL2 - 2.0*piconst*(SLMU(2)/(mE*mE-mMU*mMU)*(-NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1234  +NRlE(a,x)*NLlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1235  -NRlE(a,x)*NRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1236  +NLlE(a,x)*NRlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x)))
1237  +SLE(2)/(mMU*mMU-mE*mE)*(-NRlE(a,x)*NRlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1238  +NRlE(a,x)*NLlMU(a,x).conjugate()*mMU*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1239  -NLlE(a,x)*NLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))
1240  +NLlE(a,x)*NRlMU(a,x).conjugate()*mE*MNeig(a)*PV.B0(1.,0.,MNeig(a)*MNeig(a),mym_se_sq(x))));
1241  }
1242  }
1243 
1244  // summing up the h (0), H (1) and A (2) parts
1245  gslpp::complex B2HiggsnR = (-0.5*HpengMuEEENR0*SLE(0)/(mh*mh)-0.5*HpengMuEEENR1*SLE(1)/(mH*mH)-0.5*HpengMuEEENR2*SLE(2)/(mA*mA))/(4.0*pi*alph);
1246  gslpp::complex B2HiggsnL = (-0.5*HpengMuEEENL0*SRE(0)/(mh*mh)-0.5*HpengMuEEENL1*SRE(1)/(mH*mH)-0.5*HpengMuEEENL2*SRE(2)/(mA*mA))/(4.0*pi*alph);
1247  gslpp::complex B3HiggsnR = (HpengMuEEENR0*SRE(0)/(mh*mh)+HpengMuEEENR1*SRE(1)/(mH*mH)+HpengMuEEENR2*SRE(2)/(mA*mA))/(4.0*pi*alph);
1248  gslpp::complex B3HiggsnL = (HpengMuEEENL0*SLE(0)/(mh*mh)+HpengMuEEENL1*SLE(1)/(mH*mH)+HpengMuEEENL2*SLE(2)/(mA*mA))/(4.0*pi*alph);
1249 
1250  // Chargino contributions
1251  gslpp::complex HpengMuEEECR0 = 0.0;
1252  gslpp::complex HpengMuEEECL0 = 0.0;
1253  gslpp::complex HpengMuEEECR1 = 0.0;
1254  gslpp::complex HpengMuEEECL1 = 0.0;
1255  gslpp::complex HpengMuEEECR2 = 0.0;
1256  gslpp::complex HpengMuEEECL2 = 0.0;
1257  for (int x=0;x<3;x++) {
1258  for (int a=0;a<2;a++) {
1259  for (int b=0;b<2;b++) {
1260  // h R contribution
1261  HpengMuEEECR0 = HpengMuEEECR0 - 2.0*piconst*(CRlE(a,x)*WL0(a,b)*CLlMU(b,x).conjugate()*(PV.B0(1.,0.,MChi(a)*MChi(a),MChi(b)*MChi(b))
1262  -mym_sn_sq(x)*PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1263  +mMU*mMU*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1264  +mE*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1265  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))))
1266  +CLlE(a,x)*WR0(a,b)*CRlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1267  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1268  +CLlE(a,x)*WR0(a,b)*CLlMU(b,x).conjugate()*mE*MChi(b)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1269  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1270  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1271  +CRlE(a,x)*WL0(a,b)*CRlMU(b,x).conjugate()*mMU*MChi(b)*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1272  +CLlE(a,x)*WL0(a,b)*CLlMU(b,x).conjugate()*mE*MChi(a)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1273  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1274  +CRlE(a,x)*WR0(a,b)*CRlMU(b,x).conjugate()*mMU*MChi(a)*(PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1275  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1276  +CRlE(a,x)*WR0(a,b)*CLlMU(b,x).conjugate()*MChi(a)*MChi(b)*(-PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))));
1277  // h L contribution
1278  HpengMuEEECL0 = HpengMuEEECL0 - 2.0*piconst*(CLlE(a,x)*WR0(a,b)*CRlMU(b,x).conjugate()*(PV.B0(1.,0.,MChi(a)*MChi(a),MChi(b)*MChi(b))
1279  -mym_sn_sq(x)*PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1280  +mMU*mMU*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1281  +mE*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1282  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))))
1283  +CRlE(a,x)*WL0(a,b)*CLlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1284  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1285  +CRlE(a,x)*WL0(a,b)*CRlMU(b,x).conjugate()*mE*MChi(b)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1286  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1287  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1288  +CLlE(a,x)*WR0(a,b)*CLlMU(b,x).conjugate()*mMU*MChi(b)*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1289  +CRlE(a,x)*WR0(a,b)*CRlMU(b,x).conjugate()*mE*MChi(a)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1290  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1291  +CLlE(a,x)*WL0(a,b)*CLlMU(b,x).conjugate()*mMU*MChi(a)*(PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1292  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1293  +CLlE(a,x)*WL0(a,b)*CRlMU(b,x).conjugate()*MChi(a)*MChi(b)*(-PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))));
1294  // H R contribution
1295  HpengMuEEECR1 = HpengMuEEECR1 - 2.0*piconst*(CRlE(a,x)*WL1(a,b)*CLlMU(b,x).conjugate()*(PV.B0(1.,0.,MChi(a)*MChi(a),MChi(b)*MChi(b))
1296  -mym_sn_sq(x)*PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1297  +mMU*mMU*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1298  +mE*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1299  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))))
1300  +CLlE(a,x)*WR1(a,b)*CRlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1301  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1302  +CLlE(a,x)*WR1(a,b)*CLlMU(b,x).conjugate()*mE*MChi(b)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1303  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1304  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1305  +CRlE(a,x)*WL1(a,b)*CRlMU(b,x).conjugate()*mMU*MChi(b)*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1306  +CLlE(a,x)*WL1(a,b)*CLlMU(b,x).conjugate()*mE*MChi(a)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1307  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1308  +CRlE(a,x)*WR1(a,b)*CRlMU(b,x).conjugate()*mMU*MChi(a)*(PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1309  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1310  +CRlE(a,x)*WR1(a,b)*CLlMU(b,x).conjugate()*MChi(a)*MChi(b)*(-PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))));
1311  // H L contribution
1312  HpengMuEEECL1 = HpengMuEEECL1 - 2.0*piconst*(CLlE(a,x)*WR1(a,b)*CRlMU(b,x).conjugate()*(PV.B0(1.,0.,MChi(a)*MChi(a),MChi(b)*MChi(b))
1313  -mym_sn_sq(x)*PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1314  +mMU*mMU*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1315  +mE*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1316  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))))
1317  +CRlE(a,x)*WL1(a,b)*CLlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1318  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1319  +CRlE(a,x)*WL1(a,b)*CRlMU(b,x).conjugate()*mE*MChi(b)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1320  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1321  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1322  +CLlE(a,x)*WR1(a,b)*CLlMU(b,x).conjugate()*mMU*MChi(b)*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1323  +CRlE(a,x)*WR1(a,b)*CRlMU(b,x).conjugate()*mE*MChi(a)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1324  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1325  +CLlE(a,x)*WL1(a,b)*CLlMU(b,x).conjugate()*mMU*MChi(a)*(PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1326  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1327  +CLlE(a,x)*WL1(a,b)*CRlMU(b,x).conjugate()*MChi(a)*MChi(b)*(-PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))));
1328  // A R contribution
1329  HpengMuEEECR2 = HpengMuEEECR2 - 2.0*piconst*(CRlE(a,x)*WL2(a,b)*CLlMU(b,x).conjugate()*(PV.B0(1.,0.,MChi(a)*MChi(a),MChi(b)*MChi(b))
1330  -mym_sn_sq(x)*PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1331  +mMU*mMU*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1332  +mE*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1333  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))))
1334  +CLlE(a,x)*WR2(a,b)*CRlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1335  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1336  +CLlE(a,x)*WR2(a,b)*CLlMU(b,x).conjugate()*mE*MChi(b)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1337  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1338  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1339  +CRlE(a,x)*WL2(a,b)*CRlMU(b,x).conjugate()*mMU*MChi(b)*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1340  +CLlE(a,x)*WL2(a,b)*CLlMU(b,x).conjugate()*mE*MChi(a)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1341  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1342  +CRlE(a,x)*WR2(a,b)*CRlMU(b,x).conjugate()*mMU*MChi(a)*(PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1343  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1344  +CRlE(a,x)*WR2(a,b)*CLlMU(b,x).conjugate()*MChi(a)*MChi(b)*(-PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))));
1345  // A L contribution
1346  HpengMuEEECL2 = HpengMuEEECL2 - 2.0*piconst*(CLlE(a,x)*WR2(a,b)*CRlMU(b,x).conjugate()*(PV.B0(1.,0.,MChi(a)*MChi(a),MChi(b)*MChi(b))
1347  -mym_sn_sq(x)*PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1348  +mMU*mMU*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1349  +mE*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1350  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))))
1351  +CRlE(a,x)*WL2(a,b)*CLlMU(b,x).conjugate()*mMU*mE*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1352  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1353  +CRlE(a,x)*WL2(a,b)*CRlMU(b,x).conjugate()*mE*MChi(b)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1354  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1355  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1356  +CLlE(a,x)*WR2(a,b)*CLlMU(b,x).conjugate()*mMU*MChi(b)*PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1357  +CRlE(a,x)*WR2(a,b)*CRlMU(b,x).conjugate()*mE*MChi(a)*(PV.C11(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1358  -PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1359  +CLlE(a,x)*WL2(a,b)*CLlMU(b,x).conjugate()*mMU*MChi(a)*(PV.C12(mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))
1360  -PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b)))
1361  +CLlE(a,x)*WL2(a,b)*CRlMU(b,x).conjugate()*MChi(a)*MChi(b)*(-PV.C0(0.,mym_sn_sq(x),MChi(a)*MChi(a),MChi(b)*MChi(b))));
1362  }
1363  for (int y=0;y<3;y++) {
1364  // h R contribution
1365  HpengMuEEECR0 = HpengMuEEECR0 - 2.0*piconst*Gnu0(x,y)*(-CLlE(a,x)*CLlMU(a,y).conjugate()*mE*(PV.C11(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1366  -PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)))
1367  -CRlE(a,x)*CRlMU(a,y).conjugate()*mMU*PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1368  -CRlE(a,x)*CLlMU(a,y).conjugate()*MChi(a)*PV.C0(0.,MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)));
1369  // h L contribution
1370  HpengMuEEECL0 = HpengMuEEECL0 - 2.0*piconst*Gnu0(x,y)*(-CRlE(a,x)*CRlMU(a,y).conjugate()*mE*(PV.C11(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1371  -PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)))
1372  -CLlE(a,x)*CLlMU(a,y).conjugate()*mMU*PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1373  -CLlE(a,x)*CRlMU(a,y).conjugate()*MChi(a)*PV.C0(0.,MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)));
1374  // H R contribution
1375  HpengMuEEECR1 = HpengMuEEECR1 - 2.0*piconst*Gnu1(x,y)*(-CLlE(a,x)*CLlMU(a,y).conjugate()*mE*(PV.C11(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1376  -PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)))
1377  -CRlE(a,x)*CRlMU(a,y).conjugate()*mMU*PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1378  -CRlE(a,x)*CLlMU(a,y).conjugate()*MChi(a)*PV.C0(0.,MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)));
1379  // H L contribution
1380  HpengMuEEECL1 = HpengMuEEECL1 - 2.0*piconst*Gnu1(x,y)*(-CRlE(a,x)*CRlMU(a,y).conjugate()*mE*(PV.C11(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1381  -PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)))
1382  -CLlE(a,x)*CLlMU(a,y).conjugate()*mMU*PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1383  -CLlE(a,x)*CRlMU(a,y).conjugate()*MChi(a)*PV.C0(0.,MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)));
1384  // A R contribution
1385  HpengMuEEECR2 = HpengMuEEECR2 - 2.0*piconst*Gnu2(x,y)*(-CLlE(a,x)*CLlMU(a,y).conjugate()*mE*(PV.C11(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1386  -PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)))
1387  -CRlE(a,x)*CRlMU(a,y).conjugate()*mMU*PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1388  -CRlE(a,x)*CLlMU(a,y).conjugate()*MChi(a)*PV.C0(0.,MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)));
1389  // A L contribution
1390  HpengMuEEECL2 = HpengMuEEECL2 - 2.0*piconst*Gnu2(x,y)*(-CRlE(a,x)*CRlMU(a,y).conjugate()*mE*(PV.C11(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1391  -PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)))
1392  -CLlE(a,x)*CLlMU(a,y).conjugate()*mMU*PV.C12(MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y))
1393  -CLlE(a,x)*CRlMU(a,y).conjugate()*MChi(a)*PV.C0(0.,MChi(a)*MChi(a),mym_sn_sq(x),mym_sn_sq(y)));
1394  }
1395  // h R contribution
1396  HpengMuEEECR0 = HpengMuEEECR0 - 2.0*piconst*(SRMU(0)/(mE*mE-mMU*mMU)*(-CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1397  +CLlE(a,x)*CRlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1398  -CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1399  +CRlE(a,x)*CLlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x)))
1400  +SRE(0)/(mMU*mMU-mE*mE)*(-CLlE(a,x)*CLlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1401  +CLlE(a,x)*CRlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1402  -CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1403  +CRlE(a,x)*CLlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))));
1404  // h L contribution
1405  HpengMuEEECL0 = HpengMuEEECL0 - 2.0*piconst*(SLMU(0)/(mE*mE-mMU*mMU)*(-CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1406  +CRlE(a,x)*CLlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1407  -CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1408  +CLlE(a,x)*CRlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x)))
1409  +SLE(0)/(mMU*mMU-mE*mE)*(-CRlE(a,x)*CRlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1410  +CRlE(a,x)*CLlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1411  -CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1412  +CLlE(a,x)*CRlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))));
1413  // H R contribution
1414  HpengMuEEECR1 = HpengMuEEECR1 - 2.0*piconst*(SRMU(1)/(mE*mE-mMU*mMU)*(-CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1415  +CLlE(a,x)*CRlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1416  -CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1417  +CRlE(a,x)*CLlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x)))
1418  +SRE(1)/(mMU*mMU-mE*mE)*(-CLlE(a,x)*CLlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1419  +CLlE(a,x)*CRlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1420  -CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1421  +CRlE(a,x)*CLlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))));
1422  // H L contribution
1423  HpengMuEEECL1 = HpengMuEEECL1 - 2.0*piconst*(SLMU(1)/(mE*mE-mMU*mMU)*(-CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1424  +CRlE(a,x)*CLlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1425  -CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1426  +CLlE(a,x)*CRlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x)))
1427  +SLE(1)/(mMU*mMU-mE*mE)*(-CRlE(a,x)*CRlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1428  +CRlE(a,x)*CLlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1429  -CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1430  +CLlE(a,x)*CRlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))));
1431  // A R contribution
1432  HpengMuEEECR2 = HpengMuEEECR2 - 2.0*piconst*(SRMU(2)/(mE*mE-mMU*mMU)*(-CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1433  +CLlE(a,x)*CRlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1434  -CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1435  +CRlE(a,x)*CLlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x)))
1436  +SRE(2)/(mMU*mMU-mE*mE)*(-CLlE(a,x)*CLlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1437  +CLlE(a,x)*CRlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1438  -CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1439  +CRlE(a,x)*CLlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))));
1440  // A L contribution
1441  HpengMuEEECL2 = HpengMuEEECL2 - 2.0*piconst*(SLMU(2)/(mE*mE-mMU*mMU)*(-CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mE*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1442  +CRlE(a,x)*CLlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1443  -CRlE(a,x)*CRlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1444  +CLlE(a,x)*CRlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x)))
1445  +SLE(2)/(mMU*mMU-mE*mE)*(-CRlE(a,x)*CRlMU(a,x).conjugate()*mMU*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1446  +CRlE(a,x)*CLlMU(a,x).conjugate()*mMU*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1447  -CLlE(a,x)*CLlMU(a,x).conjugate()*mE*mMU*PV.B1(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))
1448  +CLlE(a,x)*CRlMU(a,x).conjugate()*mE*MChi(a)*PV.B0(1.,0.,MChi(a)*MChi(a),mym_sn_sq(x))));
1449  }
1450  }
1451 
1452  // summing up the h (0), H (1) and A (2) parts
1453  gslpp::complex B2HiggscR = (-0.5*HpengMuEEECR0*SLE(0)/(mh*mh)-0.5*HpengMuEEECR1*SLE(1)/(mH*mH)-0.5*HpengMuEEECR2*SLE(2)/(mA*mA))/(4.0*pi*alph);
1454  gslpp::complex B2HiggscL = (-0.5*HpengMuEEECL0*SRE(0)/(mh*mh)-0.5*HpengMuEEECL1*SRE(1)/(mH*mH)-0.5*HpengMuEEECL2*SRE(2)/(mA*mA))/(4.0*pi*alph);
1455  gslpp::complex B3HiggscR = (HpengMuEEECR0*SRE(0)/(mh*mh)+HpengMuEEECR1*SRE(1)/(mH*mH)+HpengMuEEECR2*SRE(2)/(mA*mA))/(4.0*pi*alph);
1456  gslpp::complex B3HiggscL = (HpengMuEEECL0*SLE(0)/(mh*mh)+HpengMuEEECL1*SLE(1)/(mH*mH)+HpengMuEEECL2*SLE(2)/(mA*mA))/(4.0*pi*alph);
1457 
1458  // write B2H and B3H into a vector for mu->3e
1459  BHFunctions.assign(0, B2HiggsnR+B2HiggscR );
1460  BHFunctions.assign(1, B2HiggsnL+B2HiggscL );
1461  BHFunctions.assign(2, B3HiggsnR+B3HiggscR );
1462  BHFunctions.assign(3, B3HiggsnL+B3HiggscL );
1463  }
1464 
1465  if (li_to_lj == 2) // tau -> 3mu
1466  {
1467  // Neutralino contributions
1468  gslpp::complex HpengTauMUMUMUNR0 = 0.0;
1469  gslpp::complex HpengTauMUMUMUNL0 = 0.0;
1470  gslpp::complex HpengTauMUMUMUNR1 = 0.0;
1471  gslpp::complex HpengTauMUMUMUNL1 = 0.0;
1472  gslpp::complex HpengTauMUMUMUNR2 = 0.0;
1473  gslpp::complex HpengTauMUMUMUNL2 = 0.0;
1474  for (int x=0;x<6;x++) {
1475  for (int a=0;a<4;a++) {
1476  for (int b=0;b<4;b++) {
1477  // h R contribution
1478  HpengTauMUMUMUNR0 = HpengTauMUMUMUNR0 - 2.0*piconst*(NRlMU(a,x)*DL0(a,b)*NLlTAU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1479  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1480  +mTAU*mTAU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1481  +mMU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1482  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1483  +NLlMU(a,x)*DR0(a,b)*NRlTAU(b,x).conjugate()*mTAU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1484  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1485  +NLlMU(a,x)*DR0(a,b)*NLlTAU(b,x).conjugate()*mMU*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1486  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1487  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1488  +NRlMU(a,x)*DL0(a,b)*NRlTAU(b,x).conjugate()*mTAU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1489  +NLlMU(a,x)*DL0(a,b)*NLlTAU(b,x).conjugate()*mMU*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1490  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1491  +NRlMU(a,x)*DR0(a,b)*NRlTAU(b,x).conjugate()*mTAU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1492  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1493  +NRlMU(a,x)*DR0(a,b)*NLlTAU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1494  // h L contribution
1495  HpengTauMUMUMUNL0 = HpengTauMUMUMUNL0 - 2.0*piconst*(NLlMU(a,x)*DR0(a,b)*NRlTAU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1496  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1497  +mTAU*mTAU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1498  +mMU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1499  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1500  +NRlMU(a,x)*DL0(a,b)*NLlTAU(b,x).conjugate()*mTAU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1501  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1502  +NRlMU(a,x)*DL0(a,b)*NRlTAU(b,x).conjugate()*mMU*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1503  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1504  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1505  +NLlMU(a,x)*DR0(a,b)*NLlTAU(b,x).conjugate()*mTAU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1506  +NRlMU(a,x)*DR0(a,b)*NRlTAU(b,x).conjugate()*mMU*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1507  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1508  +NLlMU(a,x)*DL0(a,b)*NLlTAU(b,x).conjugate()*mTAU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1509  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1510  +NLlMU(a,x)*DL0(a,b)*NRlTAU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1511  // H R contribution
1512  HpengTauMUMUMUNR1 = HpengTauMUMUMUNR1 - 2.0*piconst*(NRlMU(a,x)*DL1(a,b)*NLlTAU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1513  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1514  +mTAU*mTAU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1515  +mMU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1516  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1517  +NLlMU(a,x)*DR1(a,b)*NRlTAU(b,x).conjugate()*mTAU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1518  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1519  +NLlMU(a,x)*DR1(a,b)*NLlTAU(b,x).conjugate()*mMU*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1520  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1521  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1522  +NRlMU(a,x)*DL1(a,b)*NRlTAU(b,x).conjugate()*mTAU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1523  +NLlMU(a,x)*DL1(a,b)*NLlTAU(b,x).conjugate()*mMU*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1524  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1525  +NRlMU(a,x)*DR1(a,b)*NRlTAU(b,x).conjugate()*mTAU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1526  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1527  +NRlMU(a,x)*DR1(a,b)*NLlTAU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1528  // H L contribution
1529  HpengTauMUMUMUNL1 = HpengTauMUMUMUNL1 - 2.0*piconst*(NLlMU(a,x)*DR1(a,b)*NRlTAU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1530  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1531  +mTAU*mTAU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1532  +mMU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1533  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1534  +NRlMU(a,x)*DL1(a,b)*NLlTAU(b,x).conjugate()*mTAU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1535  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1536  +NRlMU(a,x)*DL1(a,b)*NRlTAU(b,x).conjugate()*mMU*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1537  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1538  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1539  +NLlMU(a,x)*DR1(a,b)*NLlTAU(b,x).conjugate()*mTAU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1540  +NRlMU(a,x)*DR1(a,b)*NRlTAU(b,x).conjugate()*mMU*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1541  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1542  +NLlMU(a,x)*DL1(a,b)*NLlTAU(b,x).conjugate()*mTAU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1543  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1544  +NLlMU(a,x)*DL1(a,b)*NRlTAU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1545  // A R contribution
1546  HpengTauMUMUMUNR2 = HpengTauMUMUMUNR2 - 2.0*piconst*(NRlMU(a,x)*DL2(a,b)*NLlTAU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1547  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1548  +mTAU*mTAU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1549  +mMU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1550  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1551  +NLlMU(a,x)*DR2(a,b)*NRlTAU(b,x).conjugate()*mTAU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1552  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1553  +NLlMU(a,x)*DR2(a,b)*NLlTAU(b,x).conjugate()*mMU*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1554  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1555  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1556  +NRlMU(a,x)*DL2(a,b)*NRlTAU(b,x).conjugate()*mTAU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1557  +NLlMU(a,x)*DL2(a,b)*NLlTAU(b,x).conjugate()*mMU*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1558  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1559  +NRlMU(a,x)*DR2(a,b)*NRlTAU(b,x).conjugate()*mTAU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1560  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1561  +NRlMU(a,x)*DR2(a,b)*NLlTAU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1562  // A L contribution
1563  HpengTauMUMUMUNL2 = HpengTauMUMUMUNL2 - 2.0*piconst*(NLlMU(a,x)*DR2(a,b)*NRlTAU(b,x).conjugate()*(PV.B0(1.,0.,MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1564  -mym_se_sq(x)*PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1565  +mTAU*mTAU*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1566  +mMU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1567  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))))
1568  +NRlMU(a,x)*DL2(a,b)*NLlTAU(b,x).conjugate()*mTAU*mMU*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1569  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1570  +NRlMU(a,x)*DL2(a,b)*NRlTAU(b,x).conjugate()*mMU*MNeig(b)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1571  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1572  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1573  +NLlMU(a,x)*DR2(a,b)*NLlTAU(b,x).conjugate()*mTAU*MNeig(b)*PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1574  +NRlMU(a,x)*DR2(a,b)*NRlTAU(b,x).conjugate()*mMU*MNeig(a)*(PV.C11(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1575  -PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1576  +NLlMU(a,x)*DL2(a,b)*NLlTAU(b,x).conjugate()*mTAU*MNeig(a)*(PV.C12(mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))
1577  -PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b)))
1578  +NLlMU(a,x)*DL2(a,b)*NRlTAU(b,x).conjugate()*MNeig(a)*MNeig(b)*(-PV.C0(0.,mym_se_sq(x),MNeig(a)*MNeig(a),MNeig(b)*MNeig(b))));
1579  }
1580  for (int y=0;y<6;y++) {
1581  // h R contribution
1582  HpengTauMUMUMUNR0 = HpengTauMUMUMUNR0 - 2.0*piconst*Gl0(x,y)*(-NLlMU(a,x)*NLlTAU(a,y).conjugate()*mMU*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1583  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1584  -NRlMU(a,x)*NRlTAU(a,y).conjugate()*mTAU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1585  -NRlMU(a,x)*NLlTAU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1586  // h L contribution
1587  HpengTauMUMUMUNL0 = HpengTauMUMUMUNL0 - 2.0*piconst*Gl0(x,y)*(-NRlMU(a,x)*NRlTAU(a,y).conjugate()*mMU*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1588  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1589  -NLlMU(a,x)*NLlTAU(a,y).conjugate()*mTAU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1590  -NLlMU(a,x)*NRlTAU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1591  // H R contribution
1592  HpengTauMUMUMUNR1 = HpengTauMUMUMUNR1 - 2.0*piconst*Gl1(x,y)*(-NLlMU(a,x)*NLlTAU(a,y).conjugate()*mMU*(PV.C11(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1593  -PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)))
1594  -NRlMU(a,x)*NRlTAU(a,y).conjugate()*mTAU*PV.C12(MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y))
1595  -NRlMU(a,x)*NLlTAU(a,y).conjugate()*MNeig(a)*PV.C0(0.,MNeig(a)*MNeig(a),mym_se_sq(x),mym_se_sq(y)));
1596  // H L contribution
1597  HpengTauMUMUMUNL1 = HpengTauMUMUMUNL1 - 2.0*piconst*Gl1(x,y)*(-NRlMU(a,x)*NRlTAU(a,y).conjugate()*mMU*(PV.C11(MNeig(a)*MNeig(a),