Page MenuHomeHEPForge
Files F7879331

No OneTemporary
Actions

View Options

diff --git a/Models/Susy/RPV/Makefile.am b/Models/Susy/RPV/Makefile.am
--- a/Models/Susy/RPV/Makefile.am
+++ b/Models/Susy/RPV/Makefile.am
@@ -1,9 +1,14 @@
if WANT_RPV
pkglib_LTLIBRARIES = HwRPV.la
endif
HwRPV_la_SOURCES = \
RPV.cc RPV.h RPV.fh \
LLEVertex.h LLEVertex.cc\
LQDVertex.h LQDVertex.cc\
-UDDVertex.h UDDVertex.cc
+UDDVertex.h UDDVertex.cc\
+RPVFFZVertex.h RPVFFZVertex.cc\
+RPVFFWVertex.h RPVFFWVertex.cc\
+RPVWSSVertex.h RPVWSSVertex.cc\
+RPVFFSVertex.h RPVFFSVertex.cc\
+RPVWWHVertex.h RPVWWHVertex.cc
HwRPV_la_LDFLAGS = -module
diff --git a/Models/Susy/RPV/RPV.cc b/Models/Susy/RPV/RPV.cc
--- a/Models/Susy/RPV/RPV.cc
+++ b/Models/Susy/RPV/RPV.cc
@@ -1,384 +1,314 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the RPV class.
//
#include "RPV.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
void RPV::persistentOutput(PersistentOStream & os ) const {
- os << lambdaLLE_ << lambdaLQD_ << lambdaUDD_
- << LLEVertex_ << LQDVertex_ << UDDVertex_;
+ os << lambdaLLE_ << lambdaLQD_ << lambdaUDD_ << ounit(vnu_,GeV)
+ << LLEVertex_ << LQDVertex_ << UDDVertex_
+ << HiggsAMix_ << HiggsPMix_;
}
void RPV::persistentInput(PersistentIStream & is, int) {
- is >> lambdaLLE_ >> lambdaLQD_ >> lambdaUDD_
- >> LLEVertex_ >> LQDVertex_ >> UDDVertex_;
+ is >> lambdaLLE_ >> lambdaLQD_ >> lambdaUDD_ >> iunit(vnu_,GeV)
+ >> LLEVertex_ >> LQDVertex_ >> UDDVertex_
+ >> HiggsAMix_ >> HiggsPMix_;
}
ClassDescription<RPV> RPV::initRPV;
// Definition of the static class description member.
void RPV::Init() {
static ClassDocumentation<RPV> documentation
("The RPV class is the base class for the implementation of the"
" R-parity violating MSSM.");
static Reference<RPV,AbstractFFSVertex> interfaceLLEVertex
("Vertex/LLE",
"The vertex for the trillinear LLE interaction",
&RPV::LLEVertex_, false, false, true, false, false);
static Reference<RPV,AbstractFFSVertex> interfaceLQDVertex
("Vertex/LQD",
"The vertex for the trillinear LQD interaction",
&RPV::LQDVertex_, false, false, true, false, false);
static Reference<RPV,AbstractFFSVertex> interfaceUDDVertex
("Vertex/UDD",
"The vertex for the trillinear UDD interaction",
&RPV::UDDVertex_, false, false, true, false, false);
}
void RPV::extractParameters(bool checkmodel) {
MSSM::extractParameters(false);
if(checkmodel) {
map<string,ParamMap>::const_iterator pit;
pit = parameters().find("modsel");
if(pit == parameters().end()) return;
ParamMap::const_iterator it;
// nmssm or mssm
it = pit->second.find(3);
int inmssm = (it != pit->second.end()) ? int(it->second) : 0;
if(inmssm != 0)
throw Exception() << "R-parity violating MSSM model"
<< " used but NMSSM read in." << Exception::runerror;
// RPV
it = pit->second.find(4);
int irpv = (it != pit->second.end()) ? int(it->second) : 0;
if(irpv != 1) throw Exception() << "RPV model used but no RPV in input file"
<< Exception::runerror;
// CPV
it = pit->second.find(5);
int icpv = (it != pit->second.end()) ? int(it->second) : 0;
if(icpv != 0) throw Exception() << "RPV model does not support CPV"
<< Exception::runerror;
// flavour violation
it = pit->second.find(6);
int ifv = (it != pit->second.end()) ? int(it->second) : 0;
if(ifv != 0) throw Exception() << "RPV model does not support "
<< "flavour violation"
<< Exception::runerror;
}
// get the RPV parameters
+ // lambda
map<string,ParamMap>::const_iterator pit;
pit=parameters().find("rvlamlle");
if( pit != parameters().end() ) {
for(ParamMap::const_iterator it = pit->second.begin();
it!=pit->second.end();++it) {
if(it->first==-1) continue;
int i = it->first/100-1;
int k = it->first%10-1;
int j = (it->first%100)/10-1;
lambdaLLE_[i][j][k] = it->second;
}
}
+ // lambda'
pit=parameters().find("rvlamlqd");
if( pit != parameters().end() ) {
for(ParamMap::const_iterator it = pit->second.begin();
it!=pit->second.end();++it) {
if(it->first==-1) continue;
int i = it->first/100-1;
int k = it->first%10-1;
int j = (it->first%100)/10-1;
lambdaLQD_[i][j][k] = it->second;
}
}
+ // lambda''
pit=parameters().find("rvlamudd");
if( pit != parameters().end() ) {
for(ParamMap::const_iterator it = pit->second.begin();
it!=pit->second.end();++it) {
if(it->first==-1) continue;
int i = it->first/100-1;
int k = it->first%10-1;
int j = (it->first%100)/10-1;
lambdaUDD_[i][j][k] = it->second;
}
}
+ // sneutrino vevs
+ pit=parameters().find("rvsnvev");
+ vnu_.resize(3);
+ if( pit != parameters().end() ) {
+ for(ParamMap::const_iterator it = pit->second.begin();
+ it!=pit->second.end();++it) {
+ if(it->first>0) {
+ assert(it->first>=1&&it->first<=3);
+ vnu_[it->first-1] = it->second*GeV;
+ }
+ }
+ }
}
void RPV::createMixingMatrices() {
-// map<string,pair<MatrixSize, MixingVector> >::const_iterator it;
-// for(it=mixings().begin();it!=mixings().end();++it) {
-// string name=it->first;
-// // pseudo-scalar higgs mixing
-// if (name == "nmamix") {
-// createMixingMatrix(theHiggsAMix,name,it->second.second,it->second.first);
-// }
-// }
+ map<string,pair<MatrixSize, MixingVector> >::const_iterator it;
+ for(it=mixings().begin();it!=mixings().end();++it) {
+ string name=it->first;
+ cerr << "testing in mixings loop " << name << "\n";
+ // pseudo-scalar higgs mixing
+ if (name == "rvamix") {
+ cerr << "testing in mixing create A\n";
+ createMixingMatrix(HiggsAMix_,name,it->second.second,it->second.first);
+ }
+ else if (name == "rvlmix") {
+ cerr << "testing in mixing create C\n";
+ createMixingMatrix(HiggsPMix_,name,it->second.second,it->second.first);
+ }
+ }
// base class for neutralinos and charginos
MSSM::createMixingMatrices();
}
-// Block MINPAR # SUSY breaking input parameters
-// 3 1.000000000000000e+01 # tanb
-// 4 1.000000000000000e+00 # sign(mu)
-// 1 1.000000000000000e+02 # m0
-// 2 2.500000000000000e+02 # m12
-// 5 -1.000000000000000e+02 # A0
+
+
+
+
+
+
+
+
+
+
+
// # Higgs mixing
// Block alpha # Effective Higgs mixing parameter
// -1.149203839391468e-01 # alpha
-// Block stopmix # stop mixing matrix
-// 1 1 5.567136252574828e-01 # O_{11}
-// 1 2 8.307044838284376e-01 # O_{12}
-// 2 1 8.307044838284376e-01 # O_{21}
-// 2 2 -5.567136252574828e-01 # O_{22}
-
-// Block sbotmix # sbottom mixing matrix
-// 1 1 9.494595966791360e-01 # O_{11}
-// 1 2 3.138892707212089e-01 # O_{12}
-// 2 1 -3.138892707212089e-01 # O_{21}
-// 2 2 9.494595966791360e-01 # O_{22}
-
-// Block staumix # stau mixing matrix
-// 1 1 2.666587767723425e-01 # O_{11}
-// 1 2 9.637910026402394e-01 # O_{12}
-// 2 1 9.637910026402394e-01 # O_{21}
-// 2 2 -2.666587767723425e-01 # O_{22}
-
-// Block nmix # neutralino mixing matrix
-// 1 1 9.851357669950006e-01 # N_{1,1}
-// 1 2 -5.771661350693690e-02 # N_{1,2}
-// 1 3 1.517355210244914e-01 # N_{1,3}
-// 1 4 -5.614841735871662e-02 # N_{1,4}
-// 2 1 1.073787916649285e-01 # N_{2,1}
-// 2 2 9.402548967081881e-01 # N_{2,2}
-// 2 3 -2.795671827746073e-01 # N_{2,3}
-// 2 4 1.619651648729555e-01 # N_{2,4}
-// 3 1 -6.112790189200206e-02 # N_{3,1}
-// 3 2 9.103171953480492e-02 # N_{3,2}
-// 3 3 6.945942965173417e-01 # N_{3,3}
-// 3 4 7.109960399990973e-01 # N_{3,4}
-// 4 1 -1.193343843912277e-01 # N_{4,1}
-// 4 2 3.229593593319475e-01 # N_{4,2}
-// 4 3 6.452575340284490e-01 # N_{4,3}
-// 4 4 -6.819818705078531e-01 # N_{4,4}
-
-// Block Umix # chargino U mixing matrix
-// 1 1 9.140759450766424e-01 # U_{1,1}
-// 1 2 -4.055430515151790e-01 # U_{1,2}
-// 2 1 4.055430515151790e-01 # U_{2,1}
-// 2 2 9.140759450766424e-01 # U_{2,2}
-
-// Block Vmix # chargino V mixing matrix
-// 1 1 9.711123714587470e-01 # V_{1,1}
-// 1 2 -2.386226351370898e-01 # V_{1,2}
-// 2 1 2.386226351370898e-01 # V_{2,1}
-// 2 2 9.711123714587470e-01 # V_{2,2}
-
// Block RVT Q= 9.118760000000000e+01 # R-Parity violating LLE soft terms
// 1 1 1 0.000000000000000e+00 # T_{111}
// 1 1 2 0.000000000000000e+00 # T_{112}
// 1 1 3 0.000000000000000e+00 # T_{113}
// 1 2 1 0.000000000000000e+00 # T_{121}
// 1 2 2 0.000000000000000e+00 # T_{122}
// 1 2 3 4.243543810455972e+01 # T_{123}
// 1 3 1 0.000000000000000e+00 # T_{131}
// 1 3 2 0.000000000000000e+00 # T_{132}
// 1 3 3 0.000000000000000e+00 # T_{133}
// 2 1 1 0.000000000000000e+00 # T_{211}
// 2 1 2 0.000000000000000e+00 # T_{212}
// 2 1 3 -4.243543810455972e+01 # T_{213}
// 2 2 1 0.000000000000000e+00 # T_{221}
// 2 2 2 0.000000000000000e+00 # T_{222}
// 2 2 3 0.000000000000000e+00 # T_{223}
// 2 3 1 0.000000000000000e+00 # T_{231}
// 2 3 2 0.000000000000000e+00 # T_{232}
// 2 3 3 0.000000000000000e+00 # T_{233}
// 3 1 1 0.000000000000000e+00 # T_{311}
// 3 1 2 0.000000000000000e+00 # T_{312}
// 3 1 3 0.000000000000000e+00 # T_{313}
// 3 2 1 0.000000000000000e+00 # T_{321}
// 3 2 2 0.000000000000000e+00 # T_{322}
// 3 2 3 0.000000000000000e+00 # T_{323}
// 3 3 1 0.000000000000000e+00 # T_{331}
// 3 3 2 0.000000000000000e+00 # T_{332}
// 3 3 3 0.000000000000000e+00 # T_{333}
// Block RVTP Q= 9.118760000000000e+01 # R-Parity violating LQD soft terms
// 1 1 1 0.000000000000000e+00 # T'_{111}
// 1 1 2 0.000000000000000e+00 # T'_{112}
// 1 1 3 0.000000000000000e+00 # T'_{113}
// 1 2 1 0.000000000000000e+00 # T'_{121}
// 1 2 2 0.000000000000000e+00 # T'_{122}
// 1 2 3 0.000000000000000e+00 # T'_{123}
// 1 3 1 0.000000000000000e+00 # T'_{131}
// 1 3 2 0.000000000000000e+00 # T'_{132}
// 1 3 3 0.000000000000000e+00 # T'_{133}
// 2 1 1 0.000000000000000e+00 # T'_{211}
// 2 1 2 0.000000000000000e+00 # T'_{212}
// 2 1 3 0.000000000000000e+00 # T'_{213}
// 2 2 1 0.000000000000000e+00 # T'_{221}
// 2 2 2 0.000000000000000e+00 # T'_{222}
// 2 2 3 0.000000000000000e+00 # T'_{223}
// 2 3 1 0.000000000000000e+00 # T'_{231}
// 2 3 2 0.000000000000000e+00 # T'_{232}
// 2 3 3 0.000000000000000e+00 # T'_{233}
// 3 1 1 0.000000000000000e+00 # T'_{311}
// 3 1 2 0.000000000000000e+00 # T'_{312}
// 3 1 3 0.000000000000000e+00 # T'_{313}
// 3 2 1 0.000000000000000e+00 # T'_{321}
// 3 2 2 0.000000000000000e+00 # T'_{322}
// 3 2 3 0.000000000000000e+00 # T'_{323}
// 3 3 1 0.000000000000000e+00 # T'_{331}
// 3 3 2 0.000000000000000e+00 # T'_{332}
// 3 3 3 0.000000000000000e+00 # T'_{333}
// Block RVTPP Q= 9.118760000000000e+01 # R-Parity violating UDD soft terms
// 1 1 1 0.000000000000000e+00 # T''_{111}
// 1 1 2 0.000000000000000e+00 # T''_{112}
// 1 1 3 0.000000000000000e+00 # T''_{113}
// 1 2 1 0.000000000000000e+00 # T''_{121}
// 1 2 2 0.000000000000000e+00 # T''_{122}
// 1 2 3 0.000000000000000e+00 # T''_{123}
// 1 3 1 0.000000000000000e+00 # T''_{131}
// 1 3 2 0.000000000000000e+00 # T''_{132}
// 1 3 3 0.000000000000000e+00 # T''_{133}
// 2 1 1 0.000000000000000e+00 # T''_{211}
// 2 1 2 0.000000000000000e+00 # T''_{212}
// 2 1 3 0.000000000000000e+00 # T''_{213}
// 2 2 1 0.000000000000000e+00 # T''_{221}
// 2 2 2 0.000000000000000e+00 # T''_{222}
// 2 2 3 0.000000000000000e+00 # T''_{223}
// 2 3 1 0.000000000000000e+00 # T''_{231}
// 2 3 2 0.000000000000000e+00 # T''_{232}
// 2 3 3 0.000000000000000e+00 # T''_{233}
// 3 1 1 0.000000000000000e+00 # T''_{311}
// 3 1 2 0.000000000000000e+00 # T''_{312}
// 3 1 3 0.000000000000000e+00 # T''_{313}
// 3 2 1 0.000000000000000e+00 # T''_{321}
// 3 2 2 0.000000000000000e+00 # T''_{322}
// 3 2 3 0.000000000000000e+00 # T''_{323}
// 3 3 1 0.000000000000000e+00 # T''_{331}
// 3 3 2 0.000000000000000e+00 # T''_{332}
// 3 3 3 0.000000000000000e+00 # T''_{333}
// Block RVKAPPA Q= 9.118760000000000e+01 # R-Parity violating kappa
// 1 0.000000000000000e+00 # kappa_{1}
// 2 0.000000000000000e+00 # kappa_{2}
// 3 0.000000000000000e+00 # kappa_{3}
// Block RVD Q= 9.118760000000000e+01 # R-Parity violating D
// 1 0.000000000000000e+00 # D_{1}
// 2 0.000000000000000e+00 # D_{2}
// 3 0.000000000000000e+00 # D_{3}
// Block RVSNVEV Q= 9.118760000000000e+01 # sneutrino VEVs D
// 1 0.000000000000000e+00 # SneutrinoVev_{1}
// 2 0.000000000000000e+00 # SneutrinoVev_{2}
// 3 0.000000000000000e+00 # SneutrinoVev_{3}
// Block RVM2LH1 Q= 9.118760000000000e+01 # M2LH1
// 1 0.000000000000000e+00 # M2LH1_{1}
// 2 0.000000000000000e+00 # M2LH1_{2}
// 3 0.000000000000000e+00 # M2LH1_{3}
-// Block RVNMIX Q= 9.118760000000000e+01 # neutrino-neutralino mixing matrix
-// 1 1 1.000000000000000e+00 # N_{11}
-// 1 2 0.000000000000000e+00 # N_{12}
-// 1 3 0.000000000000000e+00 # N_{13}
-// 1 4 0.000000000000000e+00 # N_{14}
-// 1 5 0.000000000000000e+00 # N_{15}
-// 1 6 0.000000000000000e+00 # N_{16}
-// 1 7 0.000000000000000e+00 # N_{17}
-// 2 1 0.000000000000000e+00 # N_{21}
-// 2 2 1.000000000000000e+00 # N_{22}
-// 2 3 0.000000000000000e+00 # N_{23}
-// 2 4 0.000000000000000e+00 # N_{24}
-// 2 5 0.000000000000000e+00 # N_{25}
-// 2 6 0.000000000000000e+00 # N_{26}
-// 2 7 0.000000000000000e+00 # N_{27}
-// 3 1 0.000000000000000e+00 # N_{31}
-// 3 2 0.000000000000000e+00 # N_{32}
-// 3 3 1.000000000000000e+00 # N_{33}
-// 3 4 0.000000000000000e+00 # N_{34}
-// 3 5 0.000000000000000e+00 # N_{35}
-// 3 6 0.000000000000000e+00 # N_{36}
-// 3 7 0.000000000000000e+00 # N_{37}
-// 4 1 0.000000000000000e+00 # N_{41}
-// 4 2 0.000000000000000e+00 # N_{42}
-// 4 3 0.000000000000000e+00 # N_{43}
-// 4 4 9.847565328754548e-01 # N_{44}
-// 4 5 1.103472257924572e-01 # N_{45}
-// 4 6 -6.109421671307692e-02 # N_{46}
-// 4 7 -1.197729410310909e-01 # N_{47}
-// 5 1 0.000000000000000e+00 # N_{51}
-// 5 2 0.000000000000000e+00 # N_{52}
-// 5 3 0.000000000000000e+00 # N_{53}
-// 5 4 -6.117909738186333e-02 # N_{54}
-// 5 5 9.417013337872188e-01 # N_{55}
-// 5 6 9.139986887133542e-02 # N_{56}
-// 5 7 3.179650609064091e-01 # N_{57}
-// 6 1 0.000000000000000e+00 # N_{61}
-// 6 2 0.000000000000000e+00 # N_{62}
-// 6 3 0.000000000000000e+00 # N_{63}
-// 6 4 1.526096816059464e-01 # N_{64}
-// 6 5 -2.757226329309796e-01 # N_{65}
-// 6 6 6.951143481105823e-01 # N_{66}
-// 6 7 6.461449975203246e-01 # N_{67}
-// 7 1 0.000000000000000e+00 # N_{71}
-// 7 2 0.000000000000000e+00 # N_{72}
-// 7 3 0.000000000000000e+00 # N_{73}
-// 7 4 -5.676243549025390e-02 # N_{74}
-// 7 5 1.581110919350379e-01 # N_{75}
-// 7 6 7.104432445349301e-01 # N_{76}
-// 7 7 -6.834100561295584e-01 # N_{77}
// Block hmix Q= 4.658779932434547e+02 # Higgs mixing parameters
// 1 3.523836232180990e+02 # mu(Q)MSSM DRbar
// 2 9.750797435881633e+00 # tan beta(Q)MSSM DRbar
// 3 2.450549911158079e+02 # higgs vev(Q)MSSM DRbar
// 4 1.540022571859735e+05 # mA^2(Q)MSSM DRbar
// Block msoft Q= 4.658779932434547e+02 # MSSM DRbar SUSY breaking parameters
// 1 1.019307676894905e+02 # M_1(Q)
// 2 1.918845353043126e+02 # M_2(Q)
// 3 5.864985652616149e+02 # M_3(Q)
// 21 3.224468700882639e+04 # mH1^2(Q)
// 22 -1.265998776222425e+05 # mH2^2(Q)
// 31 1.929790251783819e+02 # meL(Q)
// 32 1.929762221011005e+02 # mmuL(Q)
// 33 1.943010234240223e+02 # mtauL(Q)
// 34 1.359061860029394e+02 # meR(Q)
// 35 1.358981170295141e+02 # mmuR(Q)
// 36 1.270729528972383e+02 # mtauR(Q)
// 41 5.456782453568811e+02 # mqL1(Q)
// 42 5.456765542407577e+02 # mqL2(Q)
// 43 4.973670810960051e+02 # mqL3(Q)
// 44 5.277735920289771e+02 # muR(Q)
// 45 5.277718568028397e+02 # mcR(Q)
// 46 4.228101648527092e+02 # mtR(Q)
// 47 5.256849074418041e+02 # mdR(Q)
// 48 5.256830975517062e+02 # msR(Q)
// 49 5.224091548560939e+02 # mbR(Q)
// Block au Q= 4.658779932434547e+02
// 1 1 -6.800325328454909e+02 # Au(Q)MSSM DRbar
// 2 2 -6.800290625507030e+02 # Ac(Q)MSSM DRbar
// 3 3 -5.004599927360533e+02 # At(Q)MSSM DRbar
// Block ad Q= 4.658779932434547e+02
// 1 1 -8.573073511584474e+02 # Ad(Q)MSSM DRbar
// 2 2 -8.573041034412643e+02 # As(Q)MSSM DRbar
// 3 3 -7.942004406414295e+02 # Ab(Q)MSSM DRbar
// Block ae Q= 4.658779932434547e+02
// 1 1 -2.510763383610091e+02 # Ae(Q)MSSM DRbar
// 2 2 -2.510705792217225e+02 # Amu(Q)MSSM DRbar
// 3 3 -2.478683315008348e+02 # Atau(Q)MSSM DRbar
void RPV::doinit() {
MSSM::doinit();
addVertex(LLEVertex_);
addVertex(LQDVertex_);
addVertex(UDDVertex_);
}
diff --git a/Models/Susy/RPV/RPV.h b/Models/Susy/RPV/RPV.h
--- a/Models/Susy/RPV/RPV.h
+++ b/Models/Susy/RPV/RPV.h
@@ -1,211 +1,250 @@
// -*- C++ -*-
#ifndef HERWIG_RPV_H
#define HERWIG_RPV_H
//
// This is the declaration of the RPV class.
//
#include "Herwig++/Models/Susy/MSSM.h"
#include "RPV.fh"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the RPV class.
*
* @see \ref RPVInterfaces "The interfaces"
* defined for RPV.
*/
class RPV: public MSSM {
public:
/**
* The default constructor.
*/
RPV() : lambdaLLE_(3,vector<vector<double> >(3,vector<double>(3,0.))),
lambdaLQD_(3,vector<vector<double> >(3,vector<double>(3,0.))),
lambdaUDD_(3,vector<vector<double> >(3,vector<double>(3,0.)))
{}
/**
* LLE couplings
*/
const vector<vector<vector<double> > > & lambdaLLE() {return lambdaLLE_;}
/**
* LQD couplings
*/
const vector<vector<vector<double> > > & lambdaLQD() {return lambdaLQD_;}
/**
* UDD couplings
*/
const vector<vector<vector<double> > > & lambdaUDD() {return lambdaUDD_;}
+ /**
+ * Sneutrino vevs
+ */
+ const vector<Energy> & sneutrinoVEVs() {return vnu_;}
+
+ /**
+ * Mixing matrix for the neutral CP-odd Higgs bosons
+ */
+ const MixingMatrixPtr & CPoddHiggsMix() const {
+ return HiggsAMix_;
+ }
+
+ /**
+ * Mixing matrix for the charged Higgs bosons
+ */
+ const MixingMatrixPtr & ChargedHiggsMix() const {
+ return HiggsPMix_;
+ }
+
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* Extract the parameters from the input blocks
*/
virtual void extractParameters(bool checkModel=true);
/**
* Create the mixing matrices for the model
*/
virtual void createMixingMatrices();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<RPV> initRPV;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
RPV & operator=(const RPV &);
private:
/**
* Trillinear couplings
*/
//@
/**
* LLE couplings
*/
vector<vector<vector<double> > > lambdaLLE_;
/**
* LQD couplings
*/
vector<vector<vector<double> > > lambdaLQD_;
/**
* UDD couplings
*/
vector<vector<vector<double> > > lambdaUDD_;
//@
/**
+ * Sneutrino vevs
+ */
+ vector<Energy> vnu_;
+
+ /**
* New vertices
*/
//@{
/**
* LLE vertex
*/
AbstractFFSVertexPtr LLEVertex_;
/**
* LQD vertex
*/
AbstractFFSVertexPtr LQDVertex_;
/**
* UDD vertex
*/
AbstractFFSVertexPtr UDDVertex_;
//@}
+ /**
+ * Mixing matrices
+ */
+ //@{
+ /**
+ * Pseudoscalar Higgs mixing
+ */
+ MixingMatrixPtr HiggsAMix_;
+
+ /**
+ * Charged Higgs mixing
+ */
+ MixingMatrixPtr HiggsPMix_;
+ //@}
+
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of RPV. */
template <>
struct BaseClassTrait<Herwig::RPV,1> {
/** Typedef of the first base class of RPV. */
typedef Herwig::MSSM NthBase;
};
/** This template specialization informs ThePEG about the name of
* the RPV class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::RPV>
: public ClassTraitsBase<Herwig::RPV> {
/** Return a platform-independent class name */
static string className() { return "Herwig::RPV"; }
/**
* The name of a file containing the dynamic library where the class
* RPV is implemented. It may also include several, space-separated,
* libraries if the class RPV depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwRPV.so"; }
};
/** @endcond */
}
#endif /* HERWIG_RPV_H */
diff --git a/Models/Susy/RPV/RPVFFSVertex.cc b/Models/Susy/RPV/RPVFFSVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVFFSVertex.cc
@@ -0,0 +1,787 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the RPVFFSVertex class.
+//
+
+#include "RPVFFSVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+
+RPVFFSVertex::RPVFFSVertex() : interactions_(0), mw_(ZERO),
+ _q2last(ZERO), _couplast(0.),
+ _leftlast(0.),_rightlast(0.),
+ _id1last(0), _id2last(0), _id3last(0),
+ yukawa_(true),
+ _sw(0.), _cw(0.),_sb(0.), _cb(0.),
+ vd_(ZERO), vu_(ZERO),
+ _massLast(make_pair(ZERO,ZERO)) {
+ orderInGem(1);
+ orderInGs(0);
+}
+
+
+// GOGOH
+
+// :theSij(2, vector<Complex>(2,0.0)),
+// theQij(2, vector<Complex>(2,0.0)),
+// theQijLp(4, vector<Complex>(2,0.0)),
+// theQijRp(4, vector<Complex>(2,0.0)),
+// theSijdp(4, vector<Complex>(4,0.0)),
+// theQijdp(4, vector<Complex>(4,0.0)),
+// theSa(0.0), theSb(0.0),
+// theCa(0.0), theCb(0.0)
+// theLLast(0.0), theRLast(0.0), theHLast(0),
+// theID1Last(0), theID2Last(0)
+
+// FFH
+
+// thetanb(0.0)
+// theSa(0.0), theSb(0.0),
+// theCa(0.0), theCb(0.0),
+// theFLast(make_pair(0,0)), theGlast(0.),
+//
+
+IBPtr RPVFFSVertex::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr RPVFFSVertex::fullclone() const {
+ return new_ptr(*this);
+}
+
+void RPVFFSVertex::doinit() {
+ // cast the model to the RPV one
+ tRPVPtr model = dynamic_ptr_cast<tRPVPtr>(generator()->standardModel());
+ if( !model ) throw InitException() << "RPVFFSVertex::doinit() - "
+ << "The pointer to the MSSM object is null!"
+ << Exception::abortnow;
+ // get the various mixing matrices
+ _stop = model->stopMix();
+ _sbot = model->sbottomMix();
+ _stau = model->stauMix();
+ _nmix = model->neutralinoMix();
+ _umix = model->charginoUMix();
+ _vmix = model->charginoVMix();
+ _mixH = model->CPevenHiggsMix();
+ _mixP = model->CPoddHiggsMix();
+ _mixC = model->ChargedHiggsMix();
+ if(!_stop || !_sbot ) throw InitException() << "RPVFFSVertex::doinit() - "
+ << "A squark mixing matrix pointer is null."
+ << " stop: " << _stop << " sbottom: "
+ << _sbot << Exception::abortnow;
+ if(!_nmix) throw InitException() << "RPVFFSVertex::doinit() - "
+ << "A gaugino mixing matrix pointer is null."
+ << " N: " << _nmix << " U: " << _umix
+ << " V: " << _vmix << Exception::abortnow;
+ if(! ( _stau || (!_stau&& _mixC->size().first>1)))
+ throw InitException() << "RPVFFSVertex::doinit() - "
+ << "Must have either the stau mixing matrix or it"
+ << " must be part of the charged Higgs boson mixing."
+ << Exception::abortnow;
+ // various interactions
+ // scalar Higgs bosons
+ vector<long> h0(3);
+ h0[0] = 25; h0[1] = 35;
+ if(_mixH&&_mixH->size().first>2) {
+ h0.push_back(1000012); h0.push_back(1000014); h0.push_back(1000016);
+ }
+ // pseudoscalar Higgs bosons
+ vector<long> a0(1,36);
+ if(_mixP&&_mixP->size().first>1) {
+ a0.push_back(1000017); a0.push_back(1000018); a0.push_back(1000019);
+ }
+ // charged Higgs bosons
+ vector<long> hp(1,37);
+ if(_mixC->size().first>1) {
+ hp.push_back(-1000011); hp.push_back(-1000013); hp.push_back(-1000015);
+ hp.push_back(-2000011); hp.push_back(-2000013); hp.push_back(-2000015);
+ }
+ // neutralinos
+ vector<long> neu(4);
+ neu[0] = 1000022; neu[1] = 1000023;
+ neu[2] = 1000025; neu[3] = 1000035;
+ if(_nmix->size().first>4) {
+ neu.push_back(12); neu.push_back(14); neu.push_back(16);
+ }
+ // charginos
+ vector<long> chg(2);
+ chg[0] = 1000024; chg[1] = 1000037;
+ if(_umix->size().first>2) {
+ chg.push_back(-11); chg.push_back(-13); chg.push_back(-15);
+ }
+ // FFH
+ if(interactions_==0||interactions_==1) {
+ // quarks neutral scalar
+ for ( unsigned int h = 0; h < h0.size(); ++h ) {
+ for(long ix=1;ix<7;++ix) addToList(-ix,ix,h0[h]);
+ }
+ // quarks neutral pseudoscalar
+ for ( unsigned int h = 0; h < a0.size(); ++h ) {
+ for(long ix=1;ix<7;++ix) addToList(-ix,ix,a0[h]);
+ }
+ // quarks charged higgs
+ for(unsigned int h=0;h<hp.size();++h) {
+ for(long ix=1;ix<6;ix+=2) {
+ //outgoing H+
+ addToList(-ix-1, ix, hp[h]);
+ //outgoing H-
+ addToList(-ix ,ix+1,-hp[h]);
+ }
+ }
+ // charged leptons neutral scalar
+ for ( unsigned int h = 0; h < h0.size(); ++h ) {
+ for(long ix=11;ix<16;ix+=2) addToList(-ix,ix,h0[h]);
+ }
+ // charged leptons neutral pseudoscalar
+ for ( unsigned int h = 0; h < a0.size(); ++h ) {
+ for(long ix=12;ix<16;ix+=2) addToList(-ix,ix,a0[h]);
+ }
+ // charged higgs to leptons, no mixing
+ if(_nmix->size().first<=4) {
+ for(unsigned int h=0;h<hp.size();++h) {
+ for(long ix=11;ix<16;ix+=2) {
+ //outgoing H+
+ addToList(-ix-1, ix, hp[h]);
+ //outgoing H-
+ addToList(-ix ,ix+1,-hp[h]);
+ }
+ }
+ }
+ }
+ // GOGOH
+ if(interactions_==0||interactions_==2) {
+ // neutral scalar Higgs neutralino
+ for(unsigned int i=0;i<h0.size();++i) {
+ for(unsigned int j = 0; j < neu.size(); ++j) {
+ for(unsigned int k = j; k < neu.size(); ++k) {
+ addToList(neu[j], neu[k], h0[i]);
+ }
+ }
+ }
+ // neutral pseudoscalar Higgs neutralino
+ for(unsigned int i=0;i<a0.size();++i) {
+ for(unsigned int j = 0; j < neu.size(); ++j) {
+ for(unsigned int k = j; k < neu.size(); ++k) {
+ addToList(neu[j], neu[k], a0[i]);
+ }
+ }
+ }
+ // neutral scalar Higgs chargino
+ for(unsigned int i=0;i<h0.size();++i) {
+ for(unsigned int j = 0; j < chg.size(); ++j) {
+ for(unsigned int k = 0; k < chg.size(); ++k) {
+ if(j==k&&abs(chg[j])<16) continue;
+ addToList(-chg[j], chg[k], h0[i]);
+ }
+ }
+ }
+ // neutral scalar Higgs chargino
+ for(unsigned int i=0;i<a0.size();++i) {
+ for(unsigned int j = 0; j < chg.size(); ++j) {
+ for(unsigned int k = 0; k < chg.size(); ++k) {
+ if(j==k&&abs(chg[j])<16) continue;
+ addToList(-chg[j], chg[k], a0[i]);
+ }
+ }
+ }
+ // charged higgs
+ for(unsigned int i=0;i<hp.size();++i) {
+ for(unsigned int j = 0; j < neu.size(); ++j) {
+ for(unsigned int k = 0; k < chg.size(); ++k) {
+ addToList(-chg[k], neu[j], hp[i]);
+ addToList( chg[k], neu[j],-hp[i]);
+ }
+ }
+ }
+ }
+ // neutralino sfermion
+ if(interactions_==0||interactions_==3) {
+ // quarks
+ for(unsigned int nl = 0; nl < neu.size(); ++nl) {
+ for(long ix=1;ix<7;++ix){
+ addToList( neu[nl], ix, -(1000000+ix) );
+ addToList( neu[nl], ix, -(2000000+ix) );
+ addToList( neu[nl], -ix, (1000000+ix) );
+ addToList( neu[nl], -ix, (2000000+ix) );
+ }
+ }
+ // neutrino
+ if(_nmix->size().first<=4) {
+ for(unsigned int nl = 0; nl < neu.size(); ++nl) {
+ for(long ix=11;ix<17;ix+=2) {
+ addToList( neu[nl], ix, -(1000000+ix) );
+ addToList( neu[nl], -ix, (1000000+ix) );
+ }
+ }
+ }
+ // charged leptons no mixing
+ if(!_mixC || _mixC->size().first==1) {
+ for(unsigned int nl = 0; nl < neu.size(); ++nl) {
+ for(long ix=11;ix<17;ix+=2) {
+ addToList( neu[nl], ix, -(1000000+ix) );
+ addToList( neu[nl], -ix, (1000000+ix) );
+ }
+ }
+ }
+ }
+ // chargino sfermion
+ if(interactions_==0||interactions_==4) {
+ //quarks
+ for(unsigned int ic = 0; ic < chg.size(); ++ic) {
+ for(long ix = 1; ix < 7; ++ix) {
+ if( ix % 2 == 0 ) {
+ addToList(-chg[ic], ix,-( 999999+ix));
+ addToList(-chg[ic], ix,-(1999999+ix));
+ addToList( chg[ic],-ix, 999999+ix );
+ addToList( chg[ic],-ix, 1999999+ix );
+ }
+ else {
+ addToList(-chg[ic],-ix, 1000001+ix );
+ addToList(-chg[ic],-ix, 2000001+ix );
+ addToList( chg[ic], ix,-(1000001+ix));
+ addToList( chg[ic], ix,-(2000001+ix));
+ }
+ }
+ }
+ // sneutrinos
+ if(!_mixH || _mixH->size().first<=2) {
+ for(unsigned int ic = 0; ic < chg.size(); ++ic) {
+ for(long ix = 11; ix < 17; ix+=2) {
+ addToList(-chg[ic],-ix,1000001+ix);
+ addToList(chg[ic],ix,-(1000001+ix));
+ }
+ }
+ }
+ // charged leptons
+ if(!_mixC || _mixC->size().first==1) {
+ for(unsigned int ic = 0; ic < chg.size(); ++ic) {
+ for(long ix = 12; ix < 17; ix+=2) {
+ addToList(-chg[ic], ix,-( 999999+ix));
+ addToList(-chg[ic], ix,-(1999999+ix));
+ addToList( chg[ic],-ix, ( 999999+ix));
+ addToList( chg[ic],-ix, (1999999+ix));
+ }
+ }
+ }
+ }
+ FFSVertex::doinit();
+ // various couplings and parameters
+ mw_ = getParticleData(ParticleID::Wplus)->mass();
+ _sw = sqrt(sin2ThetaW());
+ double tb = model_->tanBeta();
+ _cw = sqrt(1. - sqr(_sw));
+ _sb = tb/sqrt(1 + sqr(tb));
+ _cb = sqrt(1 - sqr(_sb));
+ vector<Energy> vnu = model->sneutrinoVEVs();
+ double g = electroMagneticCoupling(sqr(mw_))/_sw;
+ Energy v = 2.*mw_/g;
+ vd_ = sqrt((sqr(v)-sqr(vnu[0])-sqr(vnu[1])-sqr(vnu[2]))/
+ (1.+sqr(tb)))*g;
+ vu_ = vd_*tb;
+ // couplings of the neutral scalar Higgs to charginos
+ Energy me = getParticleData(ParticleID::eminus )->mass();
+ Energy mmu = getParticleData(ParticleID::muminus )->mass();
+ Energy mtau = getParticleData(ParticleID::tauminus)->mass();
+ for(unsigned int ih=0;ih<_mixH->size().first;++ih ) {
+ OCCHL_.push_back(vector<vector<Complex> >
+ (_umix->size().first,vector<Complex>(_umix->size().first,0.)));
+ for(unsigned int i=0;i<_umix->size().first;++i) {
+ for(unsigned int j=0;j<_umix->size().first;++j) {
+ OCCHL_[ih][i][j] = -sqrt(0.5)*
+ ( (*_mixH)(ih,0)*(*_vmix)(i,0)*(*_umix)(j,1) +
+ (*_mixH)(ih,1)*(*_vmix)(i,1)*(*_umix)(j,0));
+ if(_umix->size().first>2) {
+ OCCHL_[ih][i][j] -= sqrt(0.5)*
+ ( (*_mixH)(ih,2)*(*_vmix)(i,0)*(*_umix)(j,2) +
+ (*_mixH)(ih,3)*(*_vmix)(i,0)*(*_umix)(j,3) +
+ (*_mixH)(ih,4)*(*_vmix)(i,0)*(*_umix)(j,4) ) +
+ ( me *(*_umix)(j,2)*(*_vmix)(i,2) +
+ mmu *(*_umix)(j,3)*(*_vmix)(i,3) +
+ mtau*(*_umix)(j,4)*(*_vmix)(i,4))/vd_*(*_mixH)(ih,0)
+ -me /vd_*(*_umix)(j,1)*(*_vmix)(i,2)
+ -mmu /vd_*(*_umix)(j,1)*(*_vmix)(i,3)
+ -mtau/vd_*(*_umix)(j,1)*(*_vmix)(i,4);
+ }
+ }
+ }
+ }
+ // couplings of the neutral scalar Higgs to neutralinos
+ double tw = _sw/_cw;
+ for(unsigned int ih=0;ih<_mixH->size().first;++ih) {
+ ONNHL_.push_back(vector<vector<Complex> >
+ (_nmix->size().first,vector<Complex>(_nmix->size().first,0.)));
+ for(unsigned int i=0;i<_nmix->size().first;++i) {
+ for(unsigned int j=0;j<_nmix->size().first;++j) {
+ ONNHL_[ih][i][j] =0.;
+ for(unsigned int in=2;in<_nmix->size().first;++in) {
+ double sign = in!=3 ? 1. : -1.;
+ ONNHL_[ih][i][j] += 0.5*sign*(*_mixH)(ih,in-2)*
+ ( - (*_nmix)(i,1)*(*_nmix)(j,in)
+ - (*_nmix)(j,1)*(*_nmix)(i,in)
+ + tw*(*_nmix)(i,0)*(*_nmix)(j,in)
+ + tw*(*_nmix)(j,0)*(*_nmix)(i,in));
+ }
+ }
+ }
+ }
+ // couplings of the neutral pseudoscalar Higgs to neutralinos
+ for(unsigned int ih=0;ih<_mixP->size().first;++ih) {
+ ONNAL_.push_back(vector<vector<Complex> >
+ (_nmix->size().first,vector<Complex>(_nmix->size().first,0.)));
+ for(unsigned int i=0;i<_nmix->size().first;++i) {
+ for(unsigned int j=0;j<_nmix->size().first;++j) {
+ ONNAL_[ih][i][j] =0.;
+ for(unsigned int in=2;in<_nmix->size().first;++in) {
+ double sign = in!=3 ? 1. : -1.;
+ ONNAL_[ih][i][j] += -0.5*sign*(*_mixP)(ih,in-2)*
+ ( - (*_nmix)(i,1)*(*_nmix)(j,in)
+ - (*_nmix)(j,1)*(*_nmix)(i,in)
+ + tw*(*_nmix)(i,0)*(*_nmix)(j,in)
+ + tw*(*_nmix)(j,0)*(*_nmix)(i,in));
+ }
+ }
+ }
+ }
+
+// for(unsigned int i = 0; i < 4; ++i) {
+// for(unsigned int j = 0; j < 4; ++j) {
+// if( j < 2 ) {
+// theQijLp[i][j] = conj(nmix(i, 3)*vmix(j,0)
+// + (nmix(i,1) + nmix(i,0)*tw)*vmix(j,1)/sqrt(2));
+// theQijRp[i][j] = nmix(i, 2)*umix(j,0)
+// - (nmix(i,1) + nmix(i,0)*tw)*umix(j,1)/sqrt(2);
+// }
+// theQijdp[i][j] = 0.5*( nmix(i,2)*( nmix(j,1) - tw*nmix(j,0) )
+// + nmix(j,2)*( nmix(i,1) - tw*nmix(i,0) ) );
+// theSijdp[i][j] = 0.5*( nmix(i,3)*( nmix(j,1) - tw*nmix(j,0) )
+// + nmix(j,3)*( nmix(i,1) - tw*nmix(i,0) ) );
+// }
+// }
+}
+
+void RPVFFSVertex::persistentOutput(PersistentOStream & os) const {
+ os << interactions_ << _stop << _sbot << _stau << _umix << _vmix
+ << _nmix << _mixH << _mixP << _mixC << ounit(mw_,GeV) << yukawa_
+ << model_ << _sw << _cw << _sb << _cb << ounit(vd_,GeV) << ounit(vu_,GeV)
+ << OCCHL_ << ONNHL_ << ONNAL_;
+}
+
+void RPVFFSVertex::persistentInput(PersistentIStream & is, int) {
+ is >> interactions_ >> _stop >> _sbot >> _stau >> _umix >> _vmix
+ >> _nmix >> _mixH >> _mixP >> _mixC >> iunit(mw_,GeV) >> yukawa_
+ >> model_ >> _sw >> _cw >> _sb >> _cb >> iunit(vd_,GeV) >> iunit(vu_,GeV)
+ >> OCCHL_ >> ONNHL_ >> ONNAL_;
+}
+
+// *** Attention *** The following static variable is needed for the type
+// description system in ThePEG. Please check that the template arguments
+// are correct (the class and its base class), and that the constructor
+// arguments are correct (the class name and the name of the dynamically
+// loadable library where the class implementation can be found).
+DescribeClass<RPVFFSVertex,Helicity::FFSVertex>
+ describeHerwigRPVFFSVertex("Herwig::RPVFFSVertex", "RPVFFSVertex.so");
+
+void RPVFFSVertex::Init() {
+
+ static ClassDocumentation<RPVFFSVertex> documentation
+ ("There is no documentation for the RPVFFSVertex class");
+
+ static Switch<RPVFFSVertex,unsigned int> interfaceInteractions
+ ("Interactions",
+ "Whice interactions to include",
+ &RPVFFSVertex::interactions_, 0, false, false);
+ static SwitchOption interfaceInteractionsAll
+ (interfaceInteractions,
+ "All",
+ "Include all the interactions",
+ 0);
+ static SwitchOption interfaceInteractionsHiggsSMFermions
+ (interfaceInteractions,
+ "HiggsSMFermions",
+ "Interactions of Higgs with SM fermions",
+ 1);
+ static SwitchOption interfaceInteractionsHiggsGaugino
+ (interfaceInteractions,
+ "HiggsGaugino",
+ "Interactions of the Higgs with the gauginos",
+ 2);
+ static SwitchOption interfaceInteractionsNeutralinoSfermion
+ (interfaceInteractions,
+ "NeutralinoSfermion",
+ "Include the neutralino sfermion interactions",
+ 3);
+ static SwitchOption interfaceInteractionsCharginoSfermions
+ (interfaceInteractions,
+ "CharginoSfermions",
+ "Include the chargino sfermion interactions",
+ 4);
+
+ static Switch<RPVFFSVertex,bool> interfaceYukawa
+ ("Yukawa",
+ "Whether or not to include the Yukawa type couplings in neutralino/chargino interactions",
+ &RPVFFSVertex::yukawa_, true, false, false);
+ static SwitchOption interfaceYukawaYes
+ (interfaceYukawa,
+ "Yes",
+ "Include the terms",
+ true);
+ static SwitchOption interfaceYukawaNo
+ (interfaceYukawa,
+ "No",
+ "Don't include them",
+ false);
+}
+
+
+void RPVFFSVertex::setCoupling(Energy2 q2, tcPDPtr part1,
+ tcPDPtr part2,tcPDPtr part3) {
+ long f1ID(part1->id()), f2ID(part2->id()), isc(part3->id());
+ // gaugino sfermion
+ long ism(abs(part1->id())),ig(abs(part2->id()));
+ tcPDPtr smfermion = part1, gaugino = part2;
+ if( ism / 1000000 == 1 ) {
+ swap( ism, ig);
+ swap(smfermion,gaugino);
+ }
+ // squarks and sleptons + neutralino/chargino
+ if(((abs(isc)>=1000000&&abs(isc)<=1000006) ||
+ (abs(isc)>=2000000&&abs(isc)<=2000006)) ||
+ (((abs(isc)>=1000011&&abs(isc)<=1000016) ||
+ (abs(isc)>=2000011&&abs(isc)<=2000016)) &&
+ (!_mixC || _mixC->size().first<=1 ||
+ !_mixP || _mixP->size().first<=1 ))) {
+ if( ig != _id1last || ism != _id2last || isc != _id3last ) {
+ _id1last = ig;
+ _id2last = ism;
+ _id3last = isc;
+ // sfermion mass eigenstate
+ unsigned int alpha(isc/1000000 - 1);
+ // chargino
+ if(gaugino->charged()) {
+ // chargino
+ unsigned int ch = ig>1000000 ? (ig-1000024)/13 : (ig-9)/2;
+ Complex ul1 = (*_umix)(ch,0);
+ Complex ul2 = (*_umix)(ch,1);
+ Complex vl1 = (*_vmix)(ch,0);
+ Complex vl2 = (*_vmix)(ch,1);
+ if( ism >= 11 && ism <= 16 ) {
+ long lept = ( ism % 2 == 0 ) ? ism - 1 : ism;
+ double y = yukawa_ ?
+ double(model_->mass(q2, getParticleData(lept))/mw_/sqrt(2)/_cb) : 0.;
+ if( ism == 12 || ism == 14 ) {
+ _leftlast = 0.;
+ _rightlast = alpha == 0 ? - ul1 : y*ul2;
+ }
+ else if( ism == 16 ) {
+ _leftlast = 0.;
+ _rightlast = -ul1*(*_stau)(alpha, 0) + y*(*_stau)(alpha,1)*ul2;
+ }
+ else if( ism == 11 || ism == 13 || ism == 15 ) {
+ _leftlast = y*conj(ul2);
+ _rightlast = -vl1;
+ }
+ }
+ else {
+ double yd(0.), yu(0.);
+ if(yukawa_) {
+ if( ism % 2 == 0) {
+ yu = model_->mass(q2, getParticleData(ism))/mw_/sqrt(2)/_sb;
+ yd = model_->mass(q2, getParticleData(ism - 1))/mw_/sqrt(2)/_cb;
+ }
+ else {
+ yu = model_->mass(q2, getParticleData(ism + 1))/mw_/sqrt(2)/_sb;
+ yd = model_->mass(q2, getParticleData(ism))/mw_/sqrt(2)/_cb;
+ }
+ }
+ //heavy quarks
+ if( ism == 5 ) {
+ _leftlast = yd*conj(ul2)*(*_stop)(alpha,0);
+ _rightlast = -vl1*(*_stop)(alpha, 0) + yu*vl2*(*_stop)(alpha,1);
+ }
+ else if( ism == 6 ) {
+ _leftlast = yu*conj(vl2)*(*_sbot)(alpha,0);
+ _rightlast = -ul1*(*_sbot)(alpha, 0) + yd*ul2*(*_sbot)(alpha,1);
+ }
+ else {
+ if( alpha == 0 ) {
+ _leftlast = (ism % 2 == 0) ? yu*conj(vl2) : yd*conj(ul2);
+ _rightlast = (ism % 2 == 0) ? -ul1 : -vl1;
+ }
+ else {
+ _leftlast = 0.;
+ _rightlast = (ism % 2 == 0) ? yd*ul2 : yu*vl2;
+ }
+ }
+ }
+ }
+ // neutralino
+ else {
+ unsigned int nl = ig> 1000000 ?
+ (ig<1000025 ? ig-1000022 : (ig-1000005)/10) : (ig-4)/2;
+ // common primed neutralino matrices
+ Complex n2prime = (*_nmix)(nl,1)*_cw - (*_nmix)(nl,0)*_sw;
+ //handle neutrinos first
+ if( ism == 12 || ism == 14 || ism == 16 ) {
+ _leftlast = Complex(0., 0.);
+ _rightlast = -sqrt(0.5)*n2prime/_cw;
+ }
+ else {
+ Complex n1prime = (*_nmix)(nl,0)*_cw + (*_nmix)(nl,1)*_sw;
+ tcPDPtr smf = getParticleData(ism);
+ double qf = smf->charge()/eplus;
+ Complex bracketl = qf*_sw*( conj(n1prime) - _sw*conj(n2prime)/_cw );
+ double y = yukawa_ ? double(model_->mass(q2, smf)/2./mw_) : 0.;
+ double lambda(0.);
+ //neutralino mixing element
+ Complex nlf(0.);
+ if( ism % 2 == 0 ) {
+ y /= _sb;
+ lambda = -0.5 + qf*sqr(_sw);
+ nlf = (*_nmix)(nl,3);
+ }
+ else {
+ y /= _cb;
+ lambda = 0.5 + qf*sqr(_sw);
+ nlf = (*_nmix)(nl,2);
+ }
+ Complex bracketr = _sw*qf*n1prime - n2prime*lambda/_cw;
+ //heavy quarks/sleptons
+ if( ism == 5 || ism == 6 || ism == 15 ) {
+ Complex ma1(0.), ma2(0.);
+ if( ism == 5 ) {
+ ma1 = (*_sbot)(alpha,0);
+ ma2 = (*_sbot)(alpha,1);
+ }
+ else if( ism == 6 ) {
+ ma1 = (*_stop)(alpha,0);
+ ma2 = (*_stop)(alpha,1);
+ }
+ else {
+ ma1 = (*_stau)(alpha,0);
+ ma2 = (*_stau)(alpha,1);
+ }
+ _leftlast = y*conj(nlf)*ma1 - ma2*bracketl;
+ _rightlast = y*nlf*ma2 + ma1*bracketr;
+ }
+ else {
+ if( alpha == 0 ) {
+ _leftlast = y*conj(nlf);
+ _rightlast = bracketr;
+ }
+ else {
+ _leftlast = -bracketl;
+ _rightlast = y*nlf;
+ }
+ }
+ _leftlast *= -sqrt(2.);
+ _rightlast *= -sqrt(2.);
+ }
+ }
+ }
+ //determine the helicity order of the vertex
+ if( smfermion->id() < 0 ) {
+ left (conj(_rightlast));
+ right(conj( _leftlast));
+ }
+ else {
+ left ( _leftlast);
+ right(_rightlast);
+ }
+ norm(1.);
+ }
+ // SM quarks and Higgs
+ else if((abs(f1ID) <= 6 && abs(f2ID) <= 6) ||
+ ((abs(f1ID) >= 11 && abs(f1ID) <= 16) &&
+ (abs(f2ID) >= 11 && abs(f2ID) <= 16) &&
+ _umix->size().first==2) ) {
+ if( q2 != _q2last || _id1last != f1ID) {
+ _massLast.first = model_->mass(q2,part1);
+ _id1last = f1ID;
+ }
+ if( q2 != _q2last || _id2last != f2ID) {
+ _massLast.second = model_->mass(q2,part2);
+ _id2last = f2ID;
+ }
+ if( q2 != _q2last) _id3last = isc;
+ Complex fact(0.);
+ // scalar neutral Higgs
+ if(isc == ParticleID::h0 || isc == ParticleID::H0 ||
+ isc == ParticleID::SUSY_nu_eL || isc == ParticleID::SUSY_nu_muL ||
+ isc == ParticleID::SUSY_nu_tauL ) {
+ unsigned int ih = isc < 1000000 ? (isc-25)/10 : (isc-1000008)/2;
+ unsigned int id = abs(f1ID);
+ fact = -_massLast.first*
+ ((id%2==0) ? (*_mixH)(ih,1)/vu_ : (*_mixH)(ih,0)/vd_);
+ left (1.);
+ right(1.);
+ }
+ // pseudoscalar neutral Higgs
+ else if(isc == ParticleID::A0 || isc == 1000017 || isc == 1000018 ||
+ isc == 1000019 ) {
+ unsigned int ih = isc < 1000000 ? 0 : (isc-1000016);
+ unsigned int id = abs(f1ID);
+ if(_mixP) {
+ fact = -Complex(0., 1.)*_massLast.first*
+ ( (id%2==0) ? (*_mixP)(ih,1)/vu_ : (*_mixP)(ih,0)/vd_);
+ }
+ else {
+ fact = -Complex(0., 1.)*_massLast.first*
+ ( (id%2==0) ? _cb/vu_ : _sb/vd_);
+ }
+ left ( 1.);
+ right(-1.);
+ }
+ // charged higgs
+ else {
+ if(!_mixC) {
+ if( abs(f1ID) % 2 == 0 ) {
+ _leftlast = _massLast.first /vu_*_cb;
+ _rightlast = _massLast.second/vd_*_sb;
+ }
+ else {
+ _leftlast = _massLast.second/vu_*_cb;
+ _rightlast = _massLast.first /vd_*_sb;
+ }
+ }
+ else {
+ unsigned int ih;
+ if(abs(isc)==ParticleID::Hplus) {
+ ih = 0;
+ }
+ else {
+ isc *= -1;
+ ih = abs(isc)<2000000 ? (abs(isc)-1000009)/2 : (abs(isc)-2000003)/2;
+ }
+ if( abs(f1ID) % 2 == 0 ) {
+ _leftlast = _massLast.first /vu_*(*_mixC)(ih,1);
+ _rightlast = -_massLast.second/vd_*(*_mixC)(ih,0);
+ }
+ else {
+ _leftlast = _massLast.second/vu_*(*_mixC)(ih,1);
+ _rightlast = -_massLast.first /vd_*(*_mixC)(ih,0);
+ }
+ }
+ if( isc > 0 ) swap(_leftlast,_rightlast);
+ fact = sqrt(2.);
+ left ( _leftlast);
+ right(_rightlast);
+ }
+ norm(fact);
+ }
+ // gauginos and the Higgs
+ else {
+ if( isc == _id3last && f1ID == _id1last && f2ID == _id2last ) {
+ norm(_couplast);
+ left ( _leftlast);
+ right(_rightlast);
+ }
+ else {
+ _id1last = f1ID;
+ _id2last = f2ID;
+ _id3last = isc;
+ // scalar neutral Higgs
+ if(isc == ParticleID::h0 || isc == ParticleID::H0 ||
+ isc == ParticleID::SUSY_nu_eL || isc == ParticleID::SUSY_nu_muL ||
+ isc == ParticleID::SUSY_nu_tauL ) {
+ unsigned int ih = isc < 1000000 ? (isc-25)/10 : (isc-1000008)/2;
+ // charginos
+ if(part1->charged()) {
+ unsigned int ei = abs(f1ID)>1000000 ? (abs(f1ID)-1000024)/13 : (abs(f1ID)-9)/2;
+ unsigned int ej = abs(f2ID)>1000000 ? (abs(f2ID)-1000024)/13 : (abs(f2ID)-9)/2;
+ _leftlast = conj(OCCHL_[ih][ej][ei]);
+ _rightlast = OCCHL_[ih][ei][ej] ;
+ }
+ // neutralinos
+ else {
+ unsigned int ei = f1ID> 1000000 ?
+ (f1ID<1000025 ? f1ID-1000022 : (f1ID-1000005)/10) : (f1ID-4)/2;
+ unsigned int ej = f2ID> 1000000 ?
+ (f2ID<1000025 ? f2ID-1000022 : (f2ID-1000005)/10) : (f2ID-4)/2;
+ _leftlast = conj(ONNHL_[ih][ej][ei]);
+ _rightlast = ONNHL_[ih][ei][ej] ;
+ }
+ }
+ // pseudoscalar neutral Higgs
+ else if(isc == ParticleID::A0 || isc == 1000017 || isc == 1000018 ||
+ isc == 1000019 ) {
+ unsigned int ih = isc < 1000000 ? 0 : (isc-1000016);
+ // charginos
+ if(part1->charged()) {
+ assert(false);
+
+// else if( isc == ParticleID::A0 ) {
+// if( abs(f2ID) == ParticleID::SUSY_chi_1plus ||
+// abs(f2ID) == ParticleID::SUSY_chi_2plus ) {
+// unsigned int ei = (abs(f1ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
+// unsigned int ej = (abs(f2ID) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
+
+// theLLast = Complex(0.,1.)*( conj(theQij[ej][ei])*theSb
+// + conj(theSij[ej][ei])*theCb );
+// theRLast = -Complex(0.,1.)*(theQij[ei][ej]*theSb + theSij[ei][ej]*theCb);
+// }
+// //neutralinos
+// else {
+ }
+ // neutralinos
+ else {
+ unsigned int ei = f1ID> 1000000 ?
+ (f1ID<1000025 ? f1ID-1000022 : (f1ID-1000005)/10) : (f1ID-4)/2;
+ unsigned int ej = f2ID> 1000000 ?
+ (f2ID<1000025 ? f2ID-1000022 : (f2ID-1000005)/10) : (f2ID-4)/2;
+
+// theLLast = Complex(0.,1.)*( conj(theQijdp[ej][ei])*theSb
+// - conj(theSijdp[ej][ei])*theCb );
+// theRLast = -Complex(0.,1.)*(theQijdp[ei][ej]*theSb - theSijdp[ei][ej]*theCb);
+// }
+// }
+ }
+ }
+ // charged higgs
+ else {
+
+// unsigned int ei(0),ej(0);
+// long chg(f2ID), neu(f1ID);
+// if( abs(neu) == ParticleID::SUSY_chi_1plus ||
+// abs(neu) == ParticleID::SUSY_chi_2plus ) swap(chg, neu);
+// ej = ( abs(chg) == ParticleID::SUSY_chi_1plus) ? 0 : 1;
+// ei = neu - ParticleID::SUSY_chi_10;
+// if( ei > 1 ) ei = ( ei == 13 ) ? 3 : 2;
+// theLLast = -theQijLp[ei][ej]*theCb;
+// theRLast = -theQijRp[ei][ej]*theSb;
+// if( isc < 0 ) {
+// Complex tmp = theLLast;
+// theLLast = conj(theRLast);
+// theRLast = conj(tmp);
+// }
+// }
+
+ }
+
+
+
+
+
+
+
+
+
+ left ( _leftlast);
+ right(_rightlast);
+ }
+ }
+ // overall normalisation
+ if(q2!=_q2last || _couplast==0.) {
+ _couplast = weakCoupling(q2);
+ _q2last=q2;
+ }
+ norm(_couplast*norm());
+}
diff --git a/Models/Susy/RPV/RPVFFSVertex.h b/Models/Susy/RPV/RPVFFSVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVFFSVertex.h
@@ -0,0 +1,264 @@
+// -*- C++ -*-
+#ifndef Herwig_RPVFFSVertex_H
+#define Herwig_RPVFFSVertex_H
+//
+// This is the declaration of the RPVFFSVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/FFSVertex.h"
+#include "RPV.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the RPVFFSVertex class.
+ *
+ * @see \ref RPVFFSVertexInterfaces "The interfaces"
+ * defined for RPVFFSVertex.
+ */
+class RPVFFSVertex: public Helicity::FFSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ RPVFFSVertex();
+
+ /**
+ * Calculate the coupling for the vertex
+ * @param q2 The scale to at which evaluate the coupling.
+ * @param particle1 The first particle in the vertex.
+ * @param particle2 The second particle in the vertex.
+ * @param particle3 The third particle in the vertex.
+ */
+ virtual void setCoupling(Energy2 q2, tcPDPtr particle1, tcPDPtr particle2,
+ tcPDPtr particle3);
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr fullclone() const;
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ * @throws InitException if object could not be initialized properly.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ RPVFFSVertex & operator=(const RPVFFSVertex &);
+
+private:
+
+ /**
+ * Which interactions to include
+ */
+ unsigned int interactions_;
+
+ /**
+ * Mixing Matrices
+ */
+ //@{
+ /**
+ * Pointer to stop mixing matrix
+ */
+ tMixingMatrixPtr _stop;
+
+ /**
+ * Pointer to sbottom mixing matrix
+ */
+ tMixingMatrixPtr _sbot;
+
+ /**
+ * Pointer to stau mixing matrix
+ */
+ tMixingMatrixPtr _stau;
+
+ /**
+ * Pointer to U chargino mixing matrix
+ */
+ tMixingMatrixPtr _umix;
+
+ /**
+ * Pointer to V chargino mixing matrix
+ */
+ tMixingMatrixPtr _vmix;
+
+ /**
+ * Pointer to the neutralino mixing matrix
+ */
+ tMixingMatrixPtr _nmix;
+
+ /**
+ * Neutral scalar Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixH;
+
+ /**
+ * Neutral pseudoscalar Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * Charged Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixC;
+ //@}
+
+ /**
+ * The mass of the \f$W\f$.
+ */
+ Energy mw_;
+
+ /**
+ * The energy scale at which the coupling
+ * was last evaluated
+ */
+ Energy2 _q2last;
+
+ /**
+ * The value of the coupling at the scale last evaluated
+ */
+ Complex _couplast;
+
+ /**
+ * Store the value of the left coupling when it was last evaluated
+ */
+ Complex _leftlast;
+
+ /**
+ * Store the value of the right coupling when it was last evaluated
+ */
+ Complex _rightlast;
+
+ /**
+ * Store the id of the last gaugino to be evaluated
+ */
+ long _id1last;
+
+ /**
+ * Store the id of the last SM fermion to be evaluated
+ */
+ long _id2last;
+
+ /**
+ * Store the id of the last scalar to be evaluated
+ */
+ long _id3last;
+
+ /**
+ * Include Yukawa's ? in neutralino/chargino interactions
+ */
+ bool yukawa_;
+
+ /**
+ * Pointer to the Susy Model object
+ */
+ tRPVPtr model_;
+
+ /**
+ * \f$\sin(\theta_w)\f$
+ */
+ double _sw;
+
+ /**
+ * \f$\cos(\theta_w)\f$
+ */
+ double _cw;
+
+ /**
+ * \f$\sin(\beta)\f$
+ */
+ double _sb;
+
+ /**
+ * \f$\cos(\beta)\f$
+ */
+ double _cb;
+
+ /**
+ * VEV for down type Higgs
+ */
+ Energy vd_;
+
+ /**
+ * VEV for up type Higgs
+ */
+ Energy vu_;
+
+ /**
+ * Values of the masses
+ */
+ pair<Energy,Energy> _massLast;
+
+ /**
+ * OCCHL
+ */
+ vector<vector<vector<Complex> > > OCCHL_;
+
+ /**
+ * ONNHL
+ */
+ vector<vector<vector<Complex> > > ONNHL_;
+
+ /**
+ * ONNAL
+ */
+ vector<vector<vector<Complex> > > ONNAL_;
+};
+
+}
+
+#endif /* Herwig_RPVFFSVertex_H */
diff --git a/Models/Susy/RPV/RPVFFWVertex.cc b/Models/Susy/RPV/RPVFFWVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVFFWVertex.cc
@@ -0,0 +1,247 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the RPVFFWVertex class.
+//
+
+#include "RPVFFWVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "Herwig++/Models/StandardModel/StandardCKM.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+
+RPVFFWVertex::RPVFFWVertex() : _ckm(3,vector<Complex>(3,0.0)),
+ _sw(0.), _couplast(0.), _q2last(ZERO),
+ _id1last(0), _id2last(0), _leftlast(0.),
+ _rightlast(0.), _interactions(0) {
+ orderInGs(0);
+ orderInGem(1);
+}
+
+IBPtr RPVFFWVertex::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr RPVFFWVertex::fullclone() const {
+ return new_ptr(*this);
+}
+
+void RPVFFWVertex::doinit() {
+ // SUSY mixing matrices
+ tSusyBasePtr model = dynamic_ptr_cast<SusyBasePtr>(generator()->standardModel());
+ if(!model)
+ throw InitException() << "RPVFFWVertex::doinit() - The model pointer is null!"
+ << Exception::abortnow;
+ _theN = model->neutralinoMix();
+ _theU = model->charginoUMix();
+ _theV = model->charginoVMix();
+ if(!_theN || !_theU || ! _theV)
+ throw InitException() << "RPVFFWVertex::doinit() - "
+ << "A mixing matrix pointer is null."
+ << " N: " << _theN << " U: " << _theU << " V: "
+ << _theV << Exception::abortnow;
+ // SM interactions
+ if(_interactions==0 || _interactions==1) {
+ // particles for outgoing W-
+ // quarks
+ for(int ix=1;ix<6;ix+=2) {
+ for(int iy=2;iy<7;iy+=2) {
+ addToList(-ix, iy, -24);
+ }
+ }
+ // leptons
+ for(int ix=11;ix<17;ix+=2) {
+ addToList(-ix, ix+1, -24);
+ }
+ // particles for outgoing W+
+ // quarks
+ for(int ix=2;ix<7;ix+=2) {
+ for(int iy=1;iy<6;iy+=2) {
+ addToList(-ix, iy, 24);
+ }
+ }
+ // leptons
+ for(int ix=11;ix<17;ix+=2) {
+ addToList(-ix-1, ix, 24);
+ }
+ }
+ // neutralino and chargino
+ if(_interactions==0 || _interactions==2) {
+ vector<long> neu(4);
+ neu[0] = 1000022; neu[1] = 1000023;
+ neu[2] = 1000025; neu[3] = 1000035;
+ if(_theN->size().first==7) {
+ neu.push_back(12);
+ neu.push_back(14);
+ neu.push_back(16);
+ }
+ vector<long> cha(2);
+ cha[0] = 1000024; cha[1] = 1000037;
+ if(_theV->size().first==5) {
+ cha.push_back(11);
+ cha.push_back(13);
+ cha.push_back(15);
+ }
+ // sign == -1 outgoing W-, sign == +1 outgoing W+
+ for(int sign = -1; sign < 2; sign += 2) {
+ for(unsigned int ine = 0; ine < neu.size(); ++ine) {
+ for(unsigned int ic = 0; ic < cha.size(); ++ic ) {
+ if(ic>1&&ine>3&&ic==ine-2) continue;
+ addToList(-sign*cha[ic], neu[ine], sign*24);
+ }
+ }
+ }
+ }
+ Helicity::FFVVertex::doinit();
+ // CKM matric
+ Ptr<CKMBase>::transient_pointer CKM = model->CKM();
+ // cast the CKM object to the HERWIG one
+ ThePEG::Ptr<Herwig::StandardCKM>::transient_const_pointer
+ hwCKM = ThePEG::dynamic_ptr_cast< ThePEG::Ptr<Herwig::StandardCKM>::
+ transient_const_pointer>(CKM);
+ if(hwCKM) {
+ vector< vector<Complex > > CKM;
+ CKM = hwCKM->getUnsquaredMatrix(generator()->standardModel()->families());
+ for(unsigned int ix=0;ix<3;++ix) {
+ for(unsigned int iy=0;iy<3;++iy) {
+ _ckm[ix][iy]=CKM[ix][iy];
+ }
+ }
+ }
+ else {
+ throw Exception() << "Must have access to the Herwig::StandardCKM object"
+ << "for the CKM matrix in SMFFWVertex::doinit()"
+ << Exception::runerror;
+ }
+ _sw = sqrt(sin2ThetaW());
+}
+
+void RPVFFWVertex::persistentOutput(PersistentOStream & os) const {
+ os << _sw << _theN << _theU << _theV << _ckm << _interactions;
+}
+
+void RPVFFWVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _sw >> _theN >> _theU >> _theV >> _ckm >> _interactions;
+}
+
+
+// *** Attention *** The following static variable is needed for the type
+// description system in ThePEG. Please check that the template arguments
+// are correct (the class and its base class), and that the constructor
+// arguments are correct (the class name and the name of the dynamically
+// loadable library where the class implementation can be found).
+DescribeClass<RPVFFWVertex,Helicity::FFVVertex>
+ describeHerwigRPVFFWVertex("Herwig::RPVFFWVertex", "RPVFFWVertex.so");
+
+void RPVFFWVertex::Init() {
+
+ static ClassDocumentation<RPVFFWVertex> documentation
+ ("The couplings of the fermions to the W boson in the RPV model"
+ " with bilinear R-parity violation");
+
+
+ static Switch<RPVFFWVertex,unsigned int> interfaceInteractions
+ ("Interactions",
+ "Which interactions to include",
+ &RPVFFWVertex::_interactions, 0, false, false);
+ static SwitchOption interfaceInteractionsAll
+ (interfaceInteractions,
+ "All",
+ "Include all the interactions",
+ 0);
+ static SwitchOption interfaceInteractionsSM
+ (interfaceInteractions,
+ "SM",
+ "Only include the MS terms",
+ 1);
+ static SwitchOption interfaceInteractionsSUSY
+ (interfaceInteractions,
+ "SUSY",
+ "Include the neutralino/chargino terms",
+ 2);
+
+}
+
+void RPVFFWVertex::setCoupling(Energy2 q2,tcPDPtr part1,
+ tcPDPtr part2,tcPDPtr part3) {
+ assert(abs(part3->id()) == ParticleID::Wplus);
+ // normalization
+ // first the overall normalisation
+ if(q2 != _q2last||_couplast==0.) {
+ _couplast = weakCoupling(q2);
+ _q2last=q2;
+ }
+ norm(_couplast);
+ // left and right couplings for quarks
+ if(abs(part1->id()) <= 6) {
+ int iferm=abs(part1->id());
+ int ianti=abs(part2->id());
+ if(iferm%2!=0) swap(iferm,ianti);
+ iferm = iferm/2;
+ ianti = (ianti+1)/2;
+ assert( iferm>=1 && iferm<=3 && ianti>=1 && ianti<=3);
+ left(-sqrt(0.5)*_ckm[iferm-1][ianti-1]);
+ right(0.);
+ }
+ else {
+ long neu, cha;
+ if(part1->charged()) {
+ cha = part1->id();
+ neu = part2->id();
+ }
+ else {
+ cha = part2->id();
+ neu = part1->id();
+ }
+ if(_theV->size().first==2&&abs(neu)<=16) {
+ left(-sqrt(0.5));
+ right(0.);
+ }
+ else {
+ // assert((abs(cha) == 1000024 || abs(cha) == 1000037) &&
+ // (neu == 1000022 || neu == 1000023 ||
+ // neu == 1000025 || neu == 1000035 ||
+ // neu == 1000045) );
+
+ //
+ if(cha != _id1last || neu != _id2last) {
+ _id1last = cha;
+ _id2last = neu;
+ unsigned int eigc = abs(cha) < 16 ? (abs(cha)-11)/2+2 :
+ (abs(cha)-1000024)/13;
+ unsigned int eign = neu<=16 ? 4+(abs(neu)-12)/2 :
+ (neu<1000025 ? neu - 1000022 : (neu-1000025)/10+2);
+ _leftlast = (*_theN)(eign, 1)*conj((*_theV)(eigc, 0)) -
+ ( (*_theN)(eign, 3)*conj((*_theV)(eigc, 1))/sqrt(2));
+ _rightlast = conj((*_theN)(eign, 1))*(*_theU)(eigc, 0) +
+ ( conj((*_theN)(eign, 2))*(*_theU)(eigc, 1)/sqrt(2));
+ if(_theV->size().first==5) {
+ for(unsigned int k=0;k<3;++k)
+ _rightlast += ( conj((*_theN)(eign, 4+k))*(*_theU)(eigc, 2+k)/sqrt(2));
+
+ }
+ }
+ Complex ltemp = _leftlast;
+ Complex rtemp = _rightlast;
+ // conjugate if +ve chargino
+ if(cha>0) {
+ ltemp = conj(ltemp);
+ rtemp = conj(rtemp);
+ }
+ if((part1->id()==cha&&cha>0)||(part2->id()==cha&&cha<0)) {
+ Complex temp = ltemp;
+ ltemp = -rtemp;
+ rtemp = -temp;
+ }
+ left (ltemp);
+ right(rtemp);
+ }
+ }
+}
diff --git a/Models/Susy/RPV/RPVFFWVertex.h b/Models/Susy/RPV/RPVFFWVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVFFWVertex.h
@@ -0,0 +1,169 @@
+// -*- C++ -*-
+#ifndef Herwig_RPVFFWVertex_H
+#define Herwig_RPVFFWVertex_H
+//
+// This is the declaration of the RPVFFWVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
+#include "Herwig++/Models/Susy/SusyBase.h"
+#include "Herwig++/Models/Susy/MixingMatrix.fh"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the RPVFFWVertex class.
+ *
+ * @see \ref RPVFFWVertexInterfaces "The interfaces"
+ * defined for RPVFFWVertex.
+ */
+class RPVFFWVertex: public Helicity::FFVVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ RPVFFWVertex();
+
+ /**
+ * Calculate the couplings.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ */
+ virtual void setCoupling(Energy2 q2, tcPDPtr part1,
+ tcPDPtr part2, tcPDPtr part3);
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr fullclone() const;
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ * @throws InitException if object could not be initialized properly.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ RPVFFWVertex & operator=(const RPVFFWVertex &);
+
+private:
+
+ /**
+ * The elements of the CKM matrix.
+ */
+ vector<vector<Complex> > _ckm;
+
+ /**
+ * Store \f$sin(\theta_W)\f$
+ */
+ double _sw;
+
+ /**
+ * Store the neutralino mixing matrix
+ */
+ tMixingMatrixPtr _theN;
+
+ /**
+ * Store the U-type chargino mixing matrix
+ */
+ tMixingMatrixPtr _theU;
+
+ /**
+ * Store the V-type chargino mixing matrix
+ */
+ tMixingMatrixPtr _theV;
+
+ /**
+ * The value of the coupling when it was last evaluated
+ */
+ Complex _couplast;
+
+ /**
+ * The scale at which the coupling was last evaluated
+ */
+ Energy2 _q2last;
+
+ /**
+ * The id of the first chargino the last time the vertex was evaluated
+ */
+ long _id1last;
+
+ /**
+ * The id of the second chargino the last time the vertex was evaluated
+ */
+ long _id2last;
+
+ /**
+ * The value of the left coupling when it was last evaluated
+ */
+ Complex _leftlast;
+
+ /**
+ * The value of the right coupling when it was last evaluated
+ */
+ Complex _rightlast;
+
+ /**
+ * Which interactions to include
+ */
+ unsigned int _interactions;
+};
+
+}
+
+#endif /* Herwig_RPVFFWVertex_H */
diff --git a/Models/Susy/RPV/RPVFFZVertex.cc b/Models/Susy/RPV/RPVFFZVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVFFZVertex.cc
@@ -0,0 +1,277 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the RPVFFZVertex class.
+//
+
+#include "RPVFFZVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+
+
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+
+RPVFFZVertex::RPVFFZVertex() : _sw(0.), _cw(0.), _id1last(0),
+ _id2last(0), _q2last(), _couplast(0.),
+ _leftlast(0.), _rightlast(0.),
+ _gl(17,0.0), _gr(17,0.0), _gblast(0),
+ _interactions(0) {
+ orderInGem(1);
+ orderInGs(0);
+}
+
+IBPtr RPVFFZVertex::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr RPVFFZVertex::fullclone() const {
+ return new_ptr(*this);
+}
+
+void RPVFFZVertex::persistentOutput(PersistentOStream & os) const {
+ os << _sw << _cw << _theN << _theU << _theV
+ << _gl << _gr << _interactions;
+}
+
+void RPVFFZVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _sw >> _cw >> _theN >> _theU >> _theV
+ >> _gl >> _gr >> _interactions;
+}
+
+
+// *** Attention *** The following static variable is needed for the type
+// description system in ThePEG. Please check that the template arguments
+// are correct (the class and its base class), and that the constructor
+// arguments are correct (the class name and the name of the dynamically
+// loadable library where the class implementation can be found).
+DescribeClass<RPVFFZVertex,Helicity::FFVVertex>
+describeHerwigRPVFFZVertex("Herwig::RPVFFZVertex", "RPVFFZVertex.so");
+
+void RPVFFZVertex::Init() {
+
+ static ClassDocumentation<RPVFFZVertex> documentation
+ ("The RPVFFZVertex class implements trhe coupling of the Z to all"
+ " fermion-antifermion pairs in models with bilinear RPV.");
+
+ static Switch<RPVFFZVertex,unsigned int> interfaceInteractions
+ ("Interactions",
+ "Which interactions to include",
+ &RPVFFZVertex::_interactions, 0, false, false);
+ static SwitchOption interfaceInteractionsAll
+ (interfaceInteractions,
+ "All",
+ "Include all the interactions",
+ 0);
+ static SwitchOption interfaceInteractionsSM
+ (interfaceInteractions,
+ "SM",
+ "Only include what would have been the interactions with the SM"
+ " fermions in the absence of mixing",
+ 1);
+ static SwitchOption interfaceInteractionsNeutralino
+ (interfaceInteractions,
+ "Neutralino",
+ "Only include what would have been the interactions with the "
+ "neutralinos in the absence of mixing",
+ 2);
+ static SwitchOption interfaceInteractionsChargino
+ (interfaceInteractions,
+ "Chargino",
+ "Only include what would have been the interactions with the "
+ "charginos in the absence of mixing",
+ 3);
+}
+
+void RPVFFZVertex::doinit() {
+ // extract the mixing matrices
+ tSusyBasePtr model = dynamic_ptr_cast<SusyBasePtr>(generator()->standardModel());
+ if(!model) throw InitException() << "RPVFFZVertex::doinit() - "
+ << "The model pointer is null."
+ << Exception::abortnow;
+ _theN = model->neutralinoMix();
+ _theU = model->charginoUMix();
+ _theV = model->charginoVMix();
+ if( !_theN || !_theU || !_theV )
+ throw InitException() << "RPVFFZVertex::doinit - "
+ << "A mixing matrix pointer is null. U: "
+ << _theU << " V: " << _theV << " N: " << _theN
+ << Exception::abortnow;
+ // Standard Model fermions
+ if(_interactions==0||_interactions==1) {
+ // PDG codes for the particles
+ // the quarks
+ for(int ix=1;ix<7;++ix) {
+ addToList(-ix, ix, 23);
+ }
+ // the leptons
+ for(int ix=11;ix<17;ix+=2) {
+ addToList(-ix, ix, 23);
+ }
+ for(int ix=12;ix<17;ix+=2) {
+ if(_theN->size().first==7)
+ addToList( ix, ix, 23);
+ else
+ addToList(-ix, ix, 23);
+ }
+ }
+ // neutralinos
+ if(_interactions==0||_interactions==2) {
+ vector<long> neu(4);
+ neu[0] = 1000022; neu[1] = 1000023;
+ neu[2] = 1000025; neu[3] = 1000035;
+ if(_theN->size().first==7) {
+ neu.push_back(12);
+ neu.push_back(14);
+ neu.push_back(16);
+ }
+ for(unsigned int i = 0; i < neu.size(); ++i) {
+ for(unsigned int j = 0; j < neu.size(); ++j) {
+ if(i>3&&i==j) addToList(neu[i], neu[j], 23);
+ }
+ }
+ }
+ // charginos
+ if(_interactions==0||_interactions==3) {
+ addToList(-1000024, 1000024, 22);
+ addToList(-1000037, 1000037, 22);
+ vector<long> cha(2);
+ cha[0] = 1000024; cha[1] = 1000037;
+ if(_theV->size().first==5) {
+ cha.push_back(11);
+ cha.push_back(13);
+ cha.push_back(15);
+ }
+ for(unsigned int i = 0; i < cha.size(); ++i) {
+ for(unsigned int j = 0; j < cha.size(); ++j) {
+ if(i>1&&i==j) addToList(cha[i], cha[j], 23);
+ }
+ }
+ }
+ Helicity::FFVVertex::doinit();
+ // weak mixing
+ double sw2 = sin2ThetaW();
+ _cw = sqrt(1. - sw2);
+ _sw = sqrt( sw2 );
+ // Standard Model couplings
+ for(int ix=1;ix<4;++ix) {
+ _gl[2*ix-1] = -0.25*(model->vd() + model->ad() );
+ _gl[2*ix ] = -0.25*(model->vu() + model->au() );
+ _gl[2*ix+9 ] = -0.25*(model->ve() + model->ae() );
+ _gl[2*ix+10] = -0.25*(model->vnu() + model->anu());
+ _gr[2*ix-1] = -0.25*(model->vd() - model->ad() );
+ _gr[2*ix ] = -0.25*(model->vu() - model->au() );
+ _gr[2*ix+9 ] = -0.25*(model->ve() - model->ae() );
+ _gr[2*ix+10] = -0.25*(model->vnu() - model->anu());
+ }
+}
+
+void RPVFFZVertex::setCoupling(Energy2 q2,tcPDPtr part1,
+ tcPDPtr part2,tcPDPtr part3) {
+ // first the overall normalisation
+ if(q2!=_q2last||_couplast==0.) {
+ _couplast = electroMagneticCoupling(q2);
+ _q2last=q2;
+ }
+ long iferm1(part1->id()), iferm2(part2->id()), boson(part3->id());
+ long iferm = abs(iferm1);
+ // chargino coupling to the photon
+ if(part3->id()==ParticleID::gamma) {
+ assert(iferm == abs(iferm2));
+ _gblast = boson;
+ _id1last = iferm1;
+ _id2last = iferm2;
+ _leftlast = -1.;
+ _rightlast = -1.;
+ }
+ // coupling to the Z
+ else {
+ assert(part3->id()==ParticleID::Z0);
+ _couplast /= _sw*_cw;
+ // quarks
+ if(iferm<=6) {
+ _leftlast = _gl[iferm];
+ _rightlast = _gr[iferm];
+ }
+ // charged leptons and charginos
+ else if(part1->iCharge()!=0) {
+ if(boson != _gblast || iferm1 != _id1last || iferm2 != _id2last) {
+ _gblast = boson;
+ _id1last = iferm1;
+ _id2last = iferm2;
+ if(_theV->size().first==2&&iferm<=16) {
+ _leftlast = -_gr[iferm];
+ _rightlast = -_gl[iferm];
+ }
+ else {
+ unsigned int ic1 = abs(iferm1) < 16 ? (abs(iferm1)-11)/2+2 :
+ (abs(iferm1)-1000024)/13;
+ unsigned int ic2 = abs(iferm2) < 16 ? (abs(iferm2)-11)/2+2 :
+ (abs(iferm2)-1000024)/13;
+ _leftlast = -(*_theV)(ic1, 0)*conj((*_theV)(ic2, 0)) -
+ 0.5*(*_theV)(ic1, 1)*conj((*_theV)(ic2, 1));
+ _rightlast = -conj((*_theU)(ic1, 0))*(*_theU)(ic2, 0) -
+ 0.5*conj((*_theU)(ic1, 1))*(*_theU)(ic2, 1);
+ if(abs(iferm1) == abs(iferm2)) {
+ _leftlast += sqr(_sw);
+ _rightlast += sqr(_sw);
+ }
+ if(_theV->size().first==5) {
+ for(unsigned int ix=0;ix<3;++ix) {
+ _leftlast += -0.5*(*_theV)(ic1, 2+ix)*conj((*_theV)(ic2, 2+ix));
+ }
+ }
+ }
+ if(iferm1>0) {
+ Complex temp = _leftlast;
+ _leftlast = -_rightlast;
+ _rightlast = -temp;
+ }
+ }
+ }
+ // neutrinos and neutralinos
+ else {
+ // case where only 4x4 matrix and neutrino
+ if(_theN->size().first==4&&iferm<=16) {
+ assert(iferm==12||iferm==14||iferm==16);
+ _leftlast = _gl[iferm];
+ _rightlast = _gr[iferm];
+ }
+ // neutralino
+ else {
+ long ic1 = part2->id();
+ long ic2 = part1->id();
+ assert(ic1 == ParticleID::SUSY_chi_10 || ic1 == ParticleID::SUSY_chi_20 ||
+ ic1 == ParticleID::SUSY_chi_30 || ic1 == ParticleID::SUSY_chi_40 ||
+ abs(ic1) == 12 || abs(ic1) == 14 || abs(ic1) == 16 );
+ assert(ic2 == ParticleID::SUSY_chi_10 || ic2 == ParticleID::SUSY_chi_20 ||
+ ic2 == ParticleID::SUSY_chi_30 || ic2 == ParticleID::SUSY_chi_40 ||
+ abs(ic2) == 12 || abs(ic2) == 14 || abs(ic2) == 16 );
+ if(ic1 != _id1last || ic2 != _id2last) {
+ _id1last = ic1;
+ _id2last = ic2;
+ unsigned int neu1 = ic1<=16 ? 4+(abs(ic1)-12)/2 :
+ (ic1<1000025 ? ic1 - 1000022 : (ic1-1000025)/10+2);
+ unsigned int neu2 = ic2<=16 ? 4+(abs(ic2)-12)/2 :
+ (ic2<1000025 ? ic2 - 1000022 : (ic2-1000025)/10+2);
+ _leftlast = 0.5*( (*_theN)(neu1, 3)*conj((*_theN)(neu2, 3)) -
+ (*_theN)(neu1, 2)*conj((*_theN)(neu2, 2)) );
+ if(_theN->size().first>4) {
+ for(unsigned int k=0;k<3;++k)
+ _leftlast -= 0.5*(*_theN)(neu1, 4+k)*conj((*_theN)(neu2, 4+k));
+ }
+ _rightlast = -conj(_leftlast);
+ }
+ }
+ }
+ }
+ norm ( _couplast);
+ left ( _leftlast);
+ right(_rightlast);
+}
diff --git a/Models/Susy/RPV/RPVFFZVertex.h b/Models/Susy/RPV/RPVFFZVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVFFZVertex.h
@@ -0,0 +1,183 @@
+// -*- C++ -*-
+#ifndef Herwig_RPVFFZVertex_H
+#define Herwig_RPVFFZVertex_H
+//
+// This is the declaration of the RPVFFZVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
+#include "Herwig++/Models/Susy/SusyBase.h"
+#include "Herwig++/Models/Susy/MixingMatrix.fh"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the RPVFFZVertex class.
+ *
+ * @see \ref RPVFFZVertexInterfaces "The interfaces"
+ * defined for RPVFFZVertex.
+ */
+class RPVFFZVertex: public Helicity::FFVVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ RPVFFZVertex();
+
+ /**
+ * Calculate the couplings.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ */
+ virtual void setCoupling(Energy2 q2, tcPDPtr part1,
+ tcPDPtr part2, tcPDPtr part3);
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr fullclone() const;
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ * @throws InitException if object could not be initialized properly.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ RPVFFZVertex & operator=(const RPVFFZVertex &);
+
+private:
+
+ /**
+ * The value of \f$sin(\theta_W)\f$
+ */
+ double _sw;
+
+ /**
+ * The value of \f$cos(\theta_W)\f$
+ */
+ double _cw;
+
+ /**
+ * Store the neutralino mixing matrix
+ */
+ tMixingMatrixPtr _theN;
+
+ /**
+ * The U mixing matrix
+ */
+ tMixingMatrixPtr _theU;
+
+ /**
+ * The V mixing matrix
+ */
+ tMixingMatrixPtr _theV;
+
+ /**
+ * Store the id of the 1st neutralino when the coupling was last evaluated
+ */
+ long _id1last;
+
+ /**
+ * Store the id of the 2nd neutralino when the coupling was last evaluated
+ */
+ long _id2last;
+
+ /**
+ * Store the value at which the coupling when it was last evalutated
+ */
+ Energy2 _q2last;
+
+ /**
+ * Store the value of the coupling when it was last evalutated
+ */
+ Complex _couplast;
+
+ /**
+ * Store the value of the left-coupling when it was last evalutated
+ */
+ Complex _leftlast;
+
+ /**
+ * Store the value of the right-coupling when it was last evalutated
+ */
+ Complex _rightlast;
+
+ /**
+ * The left couplings of the Standard Model fermions.
+ */
+ vector<double> _gl;
+
+ /**
+ * The right couplings of the Standard Model fermions.
+ */
+ vector<double> _gr;
+
+ /**
+ * The ID of the gauge boson when the vertex was last evaluated
+ */
+ long _gblast;
+
+ /**
+ * Switch controlling which subset of interactions to include
+ */
+ unsigned int _interactions;
+
+};
+
+}
+
+#endif /* Herwig_RPVFFZVertex_H */
diff --git a/Models/Susy/RPV/RPVWSSVertex.cc b/Models/Susy/RPV/RPVWSSVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVWSSVertex.cc
@@ -0,0 +1,373 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the RPVWSSVertex class.
+//
+
+#include "RPVWSSVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "RPV.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
+using namespace Herwig;
+
+RPVWSSVertex::RPVWSSVertex() :_sw(0.), _cw(0.),
+ _s2w(0.), _c2w(0.), _sbma(0.), _cbma(0.),
+ _interactions(0), _q2last(),
+ _ulast(0), _dlast(0), _gblast(0),
+ _factlast(0.), _couplast(0.) {
+ orderInGs(0);
+ orderInGem(1);
+}
+
+IBPtr RPVWSSVertex::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr RPVWSSVertex::fullclone() const {
+ return new_ptr(*this);
+}
+
+void RPVWSSVertex::doinit() {
+ Helicity::VSSVertex::doinit();
+ // extract the model
+ tRPVPtr model = dynamic_ptr_cast<tRPVPtr>(generator()->standardModel());
+ if( !model ) throw InitException() << "RPVWSSVertex::doinit() - The"
+ << " pointer to the RPV object is null!"
+ << Exception::abortnow;
+ // get the mixing matrices
+ _mixH = model->CPevenHiggsMix() ;
+ _mixP = model->CPoddHiggsMix() ;
+ _mixC = model->ChargedHiggsMix();
+ // sfermion interactions
+ if(_interactions==0||_interactions==1) {
+ // squarks
+ //W-
+ //LL-squarks
+ for(long ix=1000001;ix<1000006;ix+=2) {
+ addToList(-24,ix+1,-ix);
+ }
+ //1-2 stop sbottom
+ addToList(-24,1000006,-2000005);
+ //2-1 stop sbottom
+ addToList(-24,2000006,-1000005);
+ //2-2 stop sbottom
+ addToList(-24,2000006,-2000005);
+ //W+
+ for(long ix=1000001;ix<1000006;ix+=2) {
+ addToList(24,-(ix+1),ix);
+ }
+ //1-2 stop sbottom
+ addToList(24,-1000006,2000005);
+ //2-1 stop sbottom
+ addToList(24,-2000006,1000005);
+ //2-2 stop sbottom
+ addToList(24,-2000006,2000005);
+ //LL squarks
+ for(long ix=1000001;ix<1000007;++ix) {
+ addToList(23,ix,-ix);
+ }
+ //RR squarks
+ for(long ix=2000001;ix<2000007;++ix) {
+ addToList(23,ix,-ix);
+ }
+ //L-Rbar stop
+ addToList(23,1000006,-2000006);
+ //Lbar-R stop
+ addToList(23,-1000006,2000006);
+ //L-Rbar sbottom
+ addToList(23,1000005,-2000005);
+ //Lbar-R sbottom
+ addToList(23,-1000005,2000005);
+ // gamma
+ //squarks
+ for(long ix=1000001;ix<1000007;++ix) {
+ addToList(22,ix,-ix);
+ }
+ for(long ix=2000001;ix<2000007;++ix) {
+ addToList(22,ix,-ix);
+ }
+ //sleptons
+ for(long ix=1000011;ix<1000016;ix+=2) {
+ addToList(22,ix,-ix);
+ }
+ for(long ix=2000011;ix<2000016;ix+=2) {
+ addToList(22,ix,-ix);
+ }
+ //LL-sleptons
+ for(long ix=1000011;ix<1000017;++ix) {
+ addToList(23,ix,-ix);
+ }
+ //RR-sleptons
+ for(long ix=2000011;ix<2000016;ix+=2) {
+ addToList(23,ix,-ix);
+ }
+ //L-Rbar stau
+ addToList(23,1000015,-2000015);
+ //Lbar-R stau
+ addToList(23,-1000015,2000015);
+ //LL-sleptons
+ for(long ix=1000011;ix<1000016;ix+=2) {
+ addToList(-24,-ix,ix+1);
+ }
+ //2-L stau
+ addToList(-24,-2000015,1000016);
+ //LL-sleptons
+ for(long ix=1000011;ix<1000016;ix+=2) {
+ addToList(24,ix,-ix-1);
+ }
+ //2-L stau
+ addToList(24,2000015,-1000016);
+ }
+ if(_interactions==0||_interactions==2) {
+ vector<long> pseudo(1,36);
+ vector<long> scalar(2);
+ scalar[0] = 25; scalar[1] = 35;
+ vector<long> charged(1,37);
+ if(_mixH&&_mixH->size().first>2) {
+ scalar.push_back(1000012);
+ scalar.push_back(1000014);
+ scalar.push_back(1000016);
+ }
+ if(_mixP&&_mixP->size().first>1) {
+ pseudo.push_back(1000017);
+ pseudo.push_back(1000018);
+ pseudo.push_back(1000019);
+ }
+ if(_mixC&&_mixC->size().first>2) {
+ charged.push_back(-1000011);
+ charged.push_back(-1000013);
+ charged.push_back(-1000015);
+ charged.push_back(-2000011);
+ charged.push_back(-2000013);
+ charged.push_back(-2000015);
+ }
+ // charged Higgs and photon
+ addToList(22,37,-37);
+ // charged higgs and Z
+ for(unsigned int ix=0;ix<charged.size();++ix) {
+ for(unsigned int iy=0;iy<charged.size();++iy) {
+ if(abs(charged[ix])>1000000&&abs(charged[iy])>100000 &&
+ ( charged[ix]==charged[iy] ||
+ (abs(charged[ix])/1000000==15&&abs(charged[iy])/1000000==15)))
+ continue;
+ addToList(23,charged[ix],-charged[iy]);
+ }
+ }
+ // charged higss, scalar higgs and W
+ for(unsigned int ix=0;ix<charged.size();++ix) {
+ for(unsigned int iy=0;iy<scalar.size();++iy) {
+ if(abs(charged[ix])>1000000&&abs(scalar[iy])>100000 &&
+ abs(scalar[iy])/1000000-1==abs(charged[iy])/1000000)
+ continue;
+ addToList( 24,-charged[ix],scalar[iy]);
+ addToList(-24, charged[ix],scalar[iy]);
+ }
+ }
+ // charged higss, pseudoscalar higgs and W
+ for(unsigned int ix=0;ix<charged.size();++ix) {
+ for(unsigned int iy=0;iy<pseudo.size();++iy) {
+ addToList( 24,-charged[ix],pseudo[iy]);
+ addToList(-24, charged[ix],pseudo[iy]);
+ }
+ }
+ }
+ VSSVertex::doinit();
+ // sfermion mixing
+ _stop = model->stopMix();
+ _sbottom = model->sbottomMix();
+ if(!_stop || !_sbottom)
+ throw InitException() << "RPVWSSVertex::doinit() - "
+ << "A mixing matrix pointer is null."
+ << " stop: " << _stop << " sbottom: " << _sbottom
+ << Exception::abortnow;
+ _stau = model->stauMix();
+ cerr << "testing mixing " << _mixC << "\n";
+ if(!_stau && (!_mixC || _mixC->size().first<2))
+ throw InitException() << "RPVWSSVertex::doinit() either the stau"
+ << " mixing matrix must be set or the stau"
+ << " included in mixing with the"
+ << " charged Higgs bosons" << Exception::abortnow;
+ // weak mixing
+ _sw = sqrt(sin2ThetaW());
+ _cw = sqrt( 1. - sqr(_sw) );
+ _s2w = 2.*_cw*_sw;
+ _c2w = sqr(_cw) - sqr(_sw);
+ double sina = sin(model->higgsMixingAngle());
+ double cosa = sqrt(1. - sqr(sina));
+ double tanb = model->tanBeta();
+ double sinb = tanb/sqrt(1. + sqr(tanb));
+ double cosb = sqrt( 1. - sqr(sinb) );
+ _sbma = sinb*cosa - sina*cosb;
+ _cbma = cosa*cosb + sina*sinb;
+}
+
+void RPVWSSVertex::persistentOutput(PersistentOStream & os) const {
+ os << _interactions << _mixH << _mixP << _mixC << _sw << _cw
+ << _stau << _stop << _sbottom << _s2w << _c2w << _sbma << _cbma;
+}
+
+void RPVWSSVertex::persistentInput(PersistentIStream & is, int) {
+ is >> _interactions >> _mixH >> _mixP >> _mixC >> _sw >> _cw
+ >> _stau >> _stop >> _sbottom >> _s2w >> _c2w >> _sbma >> _cbma;
+}
+
+
+// *** Attention *** The following static variable is needed for the type
+// description system in ThePEG. Please check that the template arguments
+// are correct (the class and its base class), and that the constructor
+// arguments are correct (the class name and the name of the dynamically
+// loadable library where the class implementation can be found).
+DescribeClass<RPVWSSVertex,Helicity::VSSVertex>
+ describeHerwigRPVWSSVertex("Herwig::RPVWSSVertex", "RPVWSSVertex.so");
+
+void RPVWSSVertex::Init() {
+
+ static ClassDocumentation<RPVWSSVertex> documentation
+ ("There is no documentation for the RPVWSSVertex class");
+
+ static Switch<RPVWSSVertex,unsigned int> interfaceInteractions
+ ("Interactions",
+ "Which interactions to include",
+ &RPVWSSVertex::_interactions, 0, false, false);
+ static SwitchOption interfaceInteractionsAll
+ (interfaceInteractions,
+ "All",
+ "Include both the interactions which would have been sfermion"
+ " and Higgs bosons with the gauge bosons in the MSSM",
+ 0);
+ static SwitchOption interfaceInteractionsSfermions
+ (interfaceInteractions,
+ "Sfermions",
+ "Include the sfermion interactions",
+ 1);
+ static SwitchOption interfaceInteractionsHiggs
+ (interfaceInteractions,
+ "Higgs",
+ "Include the Higgs boson interactions",
+ 2);
+
+}
+
+
+void RPVWSSVertex::setCoupling(Energy2 q2,tcPDPtr part1,
+ tcPDPtr part2,tcPDPtr part3){
+ long gboson = part1->id();
+ assert( gboson == ParticleID::Z0 ||
+ gboson == ParticleID::gamma ||
+ abs(gboson) == ParticleID::Wplus );
+ long h1ID = part2->id();
+ long h2ID = part3->id();
+ // squarks and sleptons
+ if(((abs(h1ID)>=1000001&&abs(h1ID)<=1000006)||
+ (abs(h1ID)>=2000001&&abs(h1ID)<=2000006)) ||
+ (((abs(h1ID)>=1000011&&abs(h1ID)<=1000016)||
+ (abs(h1ID)>=2000011&&abs(h1ID)<=2000016))&&
+ (!_mixP||!_mixC||!_mixH||_mixP->size().first==1||
+ _mixC->size().first==1||_mixH->size().first<=2))) {
+ long sf1(abs(part2->id())),sf2(abs(part3->id()));
+ if( sf1 % 2 != 0 ) swap(sf1, sf2);
+ if( sf1 != _ulast || sf2 != _dlast || gboson != _gblast) {
+ _gblast = gboson;
+ _ulast = sf1;
+ _dlast = sf2;
+ //photon is simplest
+ if( gboson == ParticleID::gamma )
+ _factlast = getParticleData(sf1)->charge()/eplus;
+ else {
+ //determine which helicity state
+ unsigned int alpha(sf1/1000000 - 1), beta(sf2/1000000 - 1);
+ //mixing factors
+ Complex m1a(0.), m1b(0.);
+ if( sf1 == ParticleID::SUSY_t_1 || sf1 == ParticleID::SUSY_t_2 )
+ m1a = (*_stop)(alpha, 0);
+ else if( sf1 == ParticleID::SUSY_b_1 || sf1 == ParticleID::SUSY_b_2 )
+ m1a = (*_sbottom)(alpha, 0);
+ else if( sf1 == ParticleID::SUSY_tau_1minus ||
+ sf1 == ParticleID::SUSY_tau_2minus )
+ m1a = (*_stau)(alpha, 0);
+ else
+ m1a = (alpha == 0) ? Complex(1.) : Complex(0.);
+ if( sf2 == ParticleID::SUSY_t_1 || sf2 == ParticleID::SUSY_t_2 )
+ m1b = (*_stop)(beta, 0);
+ else if( sf2 == ParticleID::SUSY_b_1 || sf2 == ParticleID::SUSY_b_2 )
+ m1b = (*_sbottom)(beta, 0);
+ else if( sf2 == ParticleID::SUSY_tau_1minus ||
+ sf2 == ParticleID::SUSY_tau_2minus )
+ m1b = (*_stau)(beta, 0);
+ else
+ m1b = (beta == 0) ? Complex(1.) : Complex(0.);
+ //W boson
+ if( gboson == ParticleID::Wplus ) {
+ _factlast = m1a*m1b/sqrt(2)/_sw;
+ }
+ //Z boson
+ else {
+ if( sf1 == ParticleID::SUSY_nu_eL || sf1 == ParticleID::SUSY_nu_muL ||
+ sf1 == ParticleID::SUSY_nu_tauL ) {
+ _factlast = 1./_cw/2./_sw;
+ }
+ else {
+ double lmda(1.);
+ if( sf2 % 2 == 0 ) lmda = -1.;
+ _factlast = lmda*m1a*m1b;
+ if( alpha == beta) {
+ double ef = getParticleData(sf1)->charge()/eplus;
+ _factlast += 2.*ef*sqr(_sw);
+ }
+ _factlast *= -1./2./_cw/_sw;
+ }
+ }
+ }
+ }
+ norm(_factlast);
+ }
+ // Higgs bosons
+ else {
+ _gblast = gboson;
+ _ulast = h1ID;
+ _dlast = h2ID;
+ _factlast = 0.;
+ if( gboson == ParticleID::Z0 ) {
+// if( abs(h1ID) == ParticleID::Hplus ) {
+// coup = theC2w/theS2w;
+// if(h1ID<0) coup *= -1.;
+// }
+// else if( h1ID == ParticleID::h0 ||
+// h2ID == ParticleID::h0 ) {
+// coup = Complex(0., 1.)*thecbma/theS2w;
+// }
+// else {
+// coup =-Complex(0., 1.)*thesbma/theS2w;
+// }
+// if(h2ID==ParticleID::A0) coup *=-1.;
+ }
+ else if( gboson == ParticleID::gamma ) {
+ _factlast = part2->iCharge()/3;
+ }
+ else {
+// long higgs = abs(h1ID) == ParticleID::Hplus ? h2ID : h1ID;
+// if( higgs == ParticleID::h0 ) {
+// coup = 0.5*thecbma/theSw;
+// }
+// else if( higgs == ParticleID::H0)
+// coup = -0.5*thesbma/theSw;
+// else
+// coup = Complex(0., 0.5)/theSw;
+// if(abs(h2ID) == ParticleID::Hplus ) coup *= -1.;
+// if(gboson<0&&higgs!=ParticleID::A0) coup *= -1.;
+ }
+ norm(_factlast);
+ }
+ if( q2 != _q2last || _couplast==0. ) {
+ _q2last = q2;
+ _couplast = electroMagneticCoupling(q2);
+ }
+ norm(_couplast*norm());
+}
diff --git a/Models/Susy/RPV/RPVWSSVertex.h b/Models/Susy/RPV/RPVWSSVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVWSSVertex.h
@@ -0,0 +1,204 @@
+// -*- C++ -*-
+#ifndef Herwig_RPVWSSVertex_H
+#define Herwig_RPVWSSVertex_H
+//
+// This is the declaration of the RPVWSSVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h"
+#include "Herwig++/Models/Susy/MixingMatrix.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the RPVWSSVertex class.
+ *
+ * @see \ref RPVWSSVertexInterfaces "The interfaces"
+ * defined for RPVWSSVertex.
+ */
+class RPVWSSVertex: public Helicity::VSSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ RPVWSSVertex();
+
+ /**
+ * Calculate the couplings.
+ * @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
+ * @param part1 The ParticleData pointer for the first particle.
+ * @param part2 The ParticleData pointer for the second particle.
+ * @param part3 The ParticleData pointer for the third particle.
+ */
+ virtual void setCoupling(Energy2 q2,tcPDPtr part1,
+ tcPDPtr part2,tcPDPtr part3);
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr fullclone() const;
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ * @throws InitException if object could not be initialized properly.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ RPVWSSVertex & operator=(const RPVWSSVertex &);
+
+private:
+
+ /**
+ * Value of \f$sin(\theta_w)\f$
+ */
+ double _sw;
+
+ /**
+ * Value of \f$cos(\theta_w)\f$
+ */
+ double _cw;
+
+ /**
+ * Neutral scalar Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixH;
+
+ /**
+ * Neutral pseudoscalar Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixP;
+
+ /**
+ * Charged Higgs mixing matrix
+ */
+ MixingMatrixPtr _mixC;
+
+ /**
+ * Stau mixing matrix
+ */
+ tMixingMatrixPtr _stau;
+
+ /**
+ * Stop mixing matrix
+ */
+ tMixingMatrixPtr _stop;
+
+ /**
+ * Sbottom mixing matrix
+ */
+ tMixingMatrixPtr _sbottom;
+
+ /**
+ * The value of \f$\sin2\theta_W\f$
+ */
+ double _s2w;
+
+ /**
+ * The value of \f$\cos2\theta_W\f$
+ */
+ double _c2w;
+
+ /**
+ * The value of \f$\sin(\beta - \alpha)\f$
+ */
+ double _sbma;
+
+ /**
+ * The value of \f$\cos(\beta - \alpha)\f$
+ */
+ double _cbma;
+
+ /**
+ * Which interactions to include
+ */
+ unsigned int _interactions;
+
+ /**
+ * Scale at which the coupling was last evaluated
+ */
+ Energy2 _q2last;
+
+ /**
+ * The up type sfermion present when the vertex was evaluated.
+ */
+ long _ulast;
+
+ /**
+ * The down type sfermion present when the vertex was evaluated.
+ */
+ long _dlast;
+
+ /**
+ * The gauge boson present when the vertex was last evaluated.
+ */
+ long _gblast;
+
+ /**
+ * The value of the mixing matrix dependent part when the vertex was
+ * last evaluated
+ */
+ Complex _factlast;
+
+ /**
+ * Value of coupling when last evaluated
+ */
+ Complex _couplast;
+};
+
+}
+
+#endif /* Herwig_RPVWSSVertex_H */
diff --git a/Models/Susy/RPV/RPVWWHVertex.cc b/Models/Susy/RPV/RPVWWHVertex.cc
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVWWHVertex.cc
@@ -0,0 +1,114 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the RPVWWHVertex class.
+//
+
+#include "RPVWWHVertex.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/EventRecord/Particle.h"
+#include "ThePEG/Repository/UseRandom.h"
+#include "ThePEG/Repository/EventGenerator.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "RPV.h"
+
+using namespace Herwig;
+
+RPVWWHVertex::RPVWWHVertex() : coupLast_(ZERO), q2Last_(ZERO) {
+ orderInGem(1);
+ orderInGs(0);
+}
+
+void RPVWWHVertex::doinit() {
+ // extract some models parameters to decide if we need sneutrinos
+ tRPVPtr model = dynamic_ptr_cast<tRPVPtr>(generator()->standardModel());
+ if( !model ) throw InitException() << "RPVWWHVertex::doinit() - "
+ << "The pointer to the RPV object is null!"
+ << Exception::abortnow;
+ // get the Higgs mixing matrix
+ MixingMatrixPtr mix = model->CPevenHiggsMix();
+ // possible interactions
+ vector<long> higgs(2);
+ higgs[0] = 25; higgs[1] = 35;
+ if(mix->size().first>2) {
+ higgs.push_back(1000012);
+ higgs.push_back(1000014);
+ higgs.push_back(1000016);
+ }
+ for(unsigned int ix=0;ix<higgs.size();++ix) {
+ addToList( 23, 23,higgs[ix]);
+ addToList(-24, 24,higgs[ix]);
+ }
+ VVSVertex::doinit();
+ // SM parameters
+ Energy mw = getParticleData(ParticleID::Wplus)->mass();
+ Energy mz = getParticleData(ParticleID::Z0)->mass();
+ double sw = sqrt(sin2ThetaW());
+ double cw = sqrt(1.-sqr(sw));
+ vector<Energy> vnu = model->sneutrinoVEVs();
+ Energy v = 2.*mw/electroMagneticCoupling(sqr(mw))*sw;
+ double tanb = model->tanBeta();
+ Energy vd = sqrt((sqr(v)-sqr(vnu[0])-sqr(vnu[1])-sqr(vnu[2]))/
+ (1.+sqr(tanb)));
+ Energy vu = vd*tanb;
+ for(unsigned int ix=0;ix<higgs.size();++ix) {
+ complex<Energy> c = vd*(*mix)(ix,0)+vu*(*mix)(ix,1);
+ for(unsigned int iy=0;iy<3;++iy) c += vnu[iy]*(*mix)(ix,2+iy);
+ vector<complex<Energy> > coup(2);
+ coup[0] = c/v*mw;
+ coup[1] = c/v*mz/cw;
+ couplings_.push_back(coup);
+ }
+ cerr << "testing vev " << v/GeV << "\n";
+}
+
+IBPtr RPVWWHVertex::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr RPVWWHVertex::fullclone() const {
+ return new_ptr(*this);
+}
+
+
+void RPVWWHVertex::persistentOutput(PersistentOStream & os) const {
+ os << ounit(couplings_,GeV);
+}
+
+void RPVWWHVertex::persistentInput(PersistentIStream & is, int) {
+ is >> iunit(couplings_,GeV);
+}
+
+
+// *** Attention *** The following static variable is needed for the type
+// description system in ThePEG. Please check that the template arguments
+// are correct (the class and its base class), and that the constructor
+// arguments are correct (the class name and the name of the dynamically
+// loadable library where the class implementation can be found).
+DescribeClass<RPVWWHVertex,Helicity::VVSVertex>
+ describeHerwigRPVWWHVertex("Herwig::RPVWWHVertex", "RPVWWHVertex.so");
+
+void RPVWWHVertex::Init() {
+
+ static ClassDocumentation<RPVWWHVertex> documentation
+ ("There is no documentation for the RPVWWHVertex class");
+
+}
+
+void RPVWWHVertex::setCoupling(Energy2 q2, tcPDPtr particle1, tcPDPtr,
+ tcPDPtr particle3) {
+ long bosonID = abs(particle1->id());
+ long higgsID = particle3->id();
+ assert( bosonID == ParticleID::Wplus || bosonID == ParticleID::Z0 );
+ int ihiggs = higgsID>1000000 ? (higgsID-1000008)/2 : (higgsID-25)/10;
+ assert(ihiggs>=0 && ihiggs<=4);
+ complex<Energy> fact = bosonID==ParticleID::Wplus ?
+ couplings_[ihiggs][0] : couplings_[ihiggs][1];
+ if( q2 != q2Last_ ) {
+ q2Last_ = q2;
+ coupLast_ = weakCoupling(q2);
+ }
+ norm(coupLast_*fact*UnitRemoval::InvE);
+}
diff --git a/Models/Susy/RPV/RPVWWHVertex.h b/Models/Susy/RPV/RPVWWHVertex.h
new file mode 100644
--- /dev/null
+++ b/Models/Susy/RPV/RPVWWHVertex.h
@@ -0,0 +1,123 @@
+// -*- C++ -*-
+#ifndef Herwig_RPVWWHVertex_H
+#define Herwig_RPVWWHVertex_H
+//
+// This is the declaration of the RPVWWHVertex class.
+//
+
+#include "ThePEG/Helicity/Vertex/Scalar/VVSVertex.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * Here is the documentation of the RPVWWHVertex class.
+ *
+ * @see \ref RPVWWHVertexInterfaces "The interfaces"
+ * defined for RPVWWHVertex.
+ */
+class RPVWWHVertex: public Helicity::VVSVertex {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ RPVWWHVertex();
+
+ /**
+ * Calculate the coupling for the vertex
+ * @param q2 The scale to at which evaluate the coupling.
+ * @param particle1 The first particle in the vertex.
+ * @param particle2 The second particle in the vertex.
+ * @param particle3 The third particle in the vertex.
+ */
+ virtual void setCoupling(Energy2 q2, tcPDPtr particle1, tcPDPtr particle2,
+ tcPDPtr particle3);
+
+public:
+
+ /** @name Functions used by the persistent I/O system. */
+ //@{
+ /**
+ * Function used to write out object persistently.
+ * @param os the persistent output stream written to.
+ */
+ void persistentOutput(PersistentOStream & os) const;
+
+ /**
+ * Function used to read in object persistently.
+ * @param is the persistent input stream read from.
+ * @param version the version number of the object when written.
+ */
+ void persistentInput(PersistentIStream & is, int version);
+ //@}
+
+ /**
+ * The standard Init function used to initialize the interfaces.
+ * Called exactly once for each class by the class description system
+ * before the main function starts or
+ * when this class is dynamically loaded.
+ */
+ static void Init();
+
+protected:
+
+ /** @name Clone Methods. */
+ //@{
+ /**
+ * Make a simple clone of this object.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr clone() const;
+
+ /** Make a clone of this object, possibly modifying the cloned object
+ * to make it sane.
+ * @return a pointer to the new object.
+ */
+ virtual IBPtr fullclone() const;
+ //@}
+
+protected:
+
+ /** @name Standard Interfaced functions. */
+ //@{
+ /**
+ * Initialize this object after the setup phase before saving an
+ * EventGenerator to disk.
+ * @throws InitException if object could not be initialized properly.
+ */
+ virtual void doinit();
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ RPVWWHVertex & operator=(const RPVWWHVertex &);
+
+private:
+
+ /**
+ * The factors for the various interactions
+ */
+ vector<vector<complex<Energy> > > couplings_;
+
+ /**
+ * last value of the coupling
+ */
+ double coupLast_;
+
+ /**
+ * last value of the scale
+ */
+ Energy2 q2Last_;
+
+};
+
+}
+
+#endif /* Herwig_RPVWWHVertex_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:58 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805947
Default Alt Text
(101 KB)

Event Timeline