diff --git a/Models/DarkMatter/DMMediatorQuarksVertex.cc b/Models/DarkMatter/DMMediatorQuarksVertex.cc
--- a/Models/DarkMatter/DMMediatorQuarksVertex.cc
+++ b/Models/DarkMatter/DMMediatorQuarksVertex.cc
@@ -1,68 +1,68 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DMMediatorQuarksVertex class.
//
#include "DMMediatorQuarksVertex.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 "DMModel.h"
using namespace Herwig;
DMMediatorQuarksVertex::DMMediatorQuarksVertex() {
orderInGem(1);
orderInGs(0);
colourStructure(ColourStructure::DELTA);
}
IBPtr DMMediatorQuarksVertex::clone() const {
return new_ptr(*this);
}
IBPtr DMMediatorQuarksVertex::fullclone() const {
return new_ptr(*this);
}
void DMMediatorQuarksVertex::persistentOutput(PersistentOStream & os) const {
os << cSMmed_;
}
void DMMediatorQuarksVertex::persistentInput(PersistentIStream & is, int) {
is >> cSMmed_;
}
void DMMediatorQuarksVertex::setCoupling(Energy2 ,tcPDPtr aa,tcPDPtr,tcPDPtr) {
int iferm=abs(aa->id());
assert(iferm>0 && iferm<4);
- norm(cSMmed_[iferm]);
+ norm(cSMmed_[iferm-1]);
left(1.);
right(1.);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<DMMediatorQuarksVertex,FFVVertex>
describeHerwigDMMediatorQuarksVertex("Herwig::DMMediatorQuarksVertex", "HwDMModel.so");
void DMMediatorQuarksVertex::Init() {
static ClassDocumentation<DMMediatorQuarksVertex> documentation
("The DMMediatorQuarksVertex class implements the coupling of the quarks to the mediator.");
}
void DMMediatorQuarksVertex::doinit() {
FFVVertex::doinit();
cDMModelPtr model = dynamic_ptr_cast<cDMModelPtr>(generator()->standardModel());
cSMmed_ = model->cSMmed();
for(int ix=1;ix<4;++ix) {
addToList(-ix, ix, 32);
}
}
diff --git a/Models/DarkMatter/DMModel.cc b/Models/DarkMatter/DMModel.cc
--- a/Models/DarkMatter/DMModel.cc
+++ b/Models/DarkMatter/DMModel.cc
@@ -1,81 +1,87 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DMModel class.
//
#include "DMModel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/ParVector.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;
DMModel::DMModel() : cDMmed_(0.), cSMmed_({1.0,1.0,1.0}) {}
IBPtr DMModel::clone() const {
return new_ptr(*this);
}
IBPtr DMModel::fullclone() const {
return new_ptr(*this);
}
void DMModel::persistentOutput(PersistentOStream & os) const {
os << cDMmed_ << cSMmed_ << QQZpVertex_ << DMDMZpVertex_;
}
void DMModel::persistentInput(PersistentIStream & is, int) {
is >> cDMmed_ >> cSMmed_ >> QQZpVertex_ >> DMDMZpVertex_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<DMModel,BSMModel>
describeHerwigDMModel("Herwig::DMModel", "HwDMModel.so");
void DMModel::Init() {
static ClassDocumentation<DMModel> documentation
("The DMModel class is designed to implement a simple dark matter model"
" with fermionic dark matter and a vector mediator, "
"as described in arXiv:1911.11147",
"The DMModel class is designed to implement a simple dark matter model"
" with fermionic dark matter and a vector mediator, "
"as described in \\cite{Plehn:2019jeo}",
"\\bibitem{Plehn:2019jeo}"
"T.~Plehn, P.~Reimitz and P.~Richardson,"
"%``Hadronic Footprint of GeV-Mass Dark Matter,''"
"arXiv:1911.11147 [hep-ph]."
"%%CITATION = ARXIV:1911.11147;%%");
static Parameter<DMModel,double> interfacecDMmed
("cDMmed",
"coupling of DM to dark mediator",
&DMModel::cDMmed_, 1.0, 0., 10., false, false, Interface::limited);
static ParVector<DMModel,double> interfacecSMmed
("cSMmed",
"coupling of SM to dark mediator",
&DMModel::cSMmed_, -1 , 1.0 , -10. , 10. , false, false, Interface::limited);
static Reference<DMModel,AbstractFFVVertex> interfaceQQZpVertex
("Vertex/QQZpVertex",
"The vertex coupling the quarks and the mediator",
&DMModel::QQZpVertex_, false, false, true, false, false);
static Reference<DMModel,AbstractFFVVertex> interfaceDMDMZpVertex
("Vertex/DMDMZpVertex",
"The vertex coupling the DM to the mediator",
&DMModel::DMDMZpVertex_, false, false, true, false, false);
}
+void DMModel::doinit() {
+ // additional vertices
+ if(QQZpVertex_) addVertex(QQZpVertex_);
+ if(DMDMZpVertex_) addVertex(DMDMZpVertex_);
+ BSMModel::doinit();
+}
diff --git a/Models/DarkMatter/DMModel.h b/Models/DarkMatter/DMModel.h
--- a/Models/DarkMatter/DMModel.h
+++ b/Models/DarkMatter/DMModel.h
@@ -1,147 +1,159 @@
// -*- C++ -*-
#ifndef Herwig_DMModel_H
#define Herwig_DMModel_H
//
// This is the declaration of the DMModel class.
//
#include "Herwig/Models/General/BSMModel.h"
#include "DMModel.fh"
namespace Herwig {
using namespace ThePEG;
/**
* The DMModel class is designed to implement a simple dark matter mode
* with fermionic dark matter and a vector mediator, as described in arXiv:1911.11147
*
* @see \ref DMModelInterfaces "The interfaces"
* defined for DMModel.
*/
class DMModel: public BSMModel {
public:
/**
* The default constructor.
*/
DMModel();
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();
public:
/**
* Access to the couplings
*/
//@{
/**
* DM coupling to the mediator
*/
const double & cDMmed() const {return cDMmed_;}
/**
* The couplings of the quarks to the mediators
*/
const vector<double> & cSMmed() const {return cSMmed_;}
//@}
public:
/**
* Pointers to the objects handling the vertices.
*/
//@{
/**
* Pointer to the quark-quark mediator vertex
*/
virtual tAbstractFFVVertexPtr vertexQQZp() const {
return QQZpVertex_;
}
/**
* Pointer to the DM-DM mediator vertex
*/
virtual tAbstractFFVVertexPtr vertexDMDMZp() const {
return DMDMZpVertex_;
}
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.
*/
DMModel & operator=(const DMModel &) = delete;
private:
/**
* DM coupling to the dark mediator
*/
double cDMmed_;
/**
* SM couplings to the dark mediator
*/
vector<double> cSMmed_;
private:
/**
* Pointer to the quark-quark mediator vertex
*/
AbstractFFVVertexPtr QQZpVertex_;
/**
* Pointer to the DM-DM mediator vertex
*/
AbstractFFVVertexPtr DMDMZpVertex_;
};
}
#endif /* Herwig_DMModel_H */
diff --git a/Models/DarkMatter/MEDM2Mesons.cc b/Models/DarkMatter/MEDM2Mesons.cc
--- a/Models/DarkMatter/MEDM2Mesons.cc
+++ b/Models/DarkMatter/MEDM2Mesons.cc
@@ -1,228 +1,268 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEDM2Mesons class.
//
#include "MEDM2Mesons.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.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 "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
+#include "Herwig/Models/General/BSMModel.h"
using namespace Herwig;
typedef LorentzVector<complex<InvEnergy> > LorentzPolarizationVectorInvE;
-MEDM2Mesons::MEDM2Mesons() : cDMmed_(0.) {
- cSMmed_ = {1.0,1.0,1.0};
-}
+MEDM2Mesons::MEDM2Mesons() : cDMmed_(0.), cSMmed_({0.0,0.0,0.0})
+{}
Energy2 MEDM2Mesons::scale() const {
return sHat();
}
unsigned int MEDM2Mesons::orderInAlphaS() const {
return 0;
}
unsigned int MEDM2Mesons::orderInAlphaEW() const {
return 0;
}
IBPtr MEDM2Mesons::clone() const {
return new_ptr(*this);
}
IBPtr MEDM2Mesons::fullclone() const {
return new_ptr(*this);
}
void MEDM2Mesons::doinit() {
// make sure the current got initialised
current_->init();
// max energy
Energy Emax = generator()->maximumCMEnergy();
// loop over the modes
int nmode=0;
for(unsigned int imode=0;imode<current_->numberOfModes();++imode) {
// get the external particles for this mode
int iq(0),ia(0);
tPDVector out = current_->particles(0,imode,iq,ia);
current_->decayModeInfo(imode,iq,ia);
if(iq==2&&ia==-2) continue;
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(incomingA_,out,1.,
incomingB_,Emax));
PhaseSpaceChannel channel(mode);
if(!current_->createMode(0,tcPDPtr(), FlavourInfo(), imode,mode,0,-1,channel,Emax)) continue;
modeMap_[imode] = nmode;
addMode(mode);
++nmode;
}
+ // cast the model
+ Ptr<BSMModel>::ptr model = dynamic_ptr_cast<Ptr<BSMModel>::ptr>(generator()->standardModel());
+ bool foundDM(false),foundU(false),foundD(false),foundS(false);
+ // find the vertices we need and extract the couplings
+ for(unsigned int ix = 0; ix < model->numberOfVertices(); ++ix ) {
+ VertexBasePtr vertex = model->vertex(ix);
+ if(vertex->getNpoint()!=3) continue;
+ for(unsigned int iloc = 0;iloc < 3; ++iloc) {
+ vector<long> ext = vertex->search(iloc, Mediator_->id());
+ if(ext.empty()) continue;
+ for(unsigned int ioff=0;ioff<ext.size();ioff+=3) {
+ if(iloc!=2) assert(false);
+ if((ext[ioff]==incomingA_->id() && ext[ioff+1]==incomingB_->id()) ||
+ (ext[ioff]==incomingB_->id() && ext[ioff+1]==incomingA_->id())) {
+ foundDM = true;
+ vertex->setCoupling(sqr(Emax),incomingA_,incomingB_,Mediator_);
+ cDMmed_ = vertex->norm();
+ }
+ else if(abs(ext[ioff])==1 && abs(ext[ioff+1])==1 && ext[ioff]==-ext[ioff+1]) {
+ foundD = true;
+ vertex->setCoupling(sqr(Emax),getParticleData(1),getParticleData(-1),Mediator_);
+ cSMmed_[0] = vertex->norm();
+ }
+ else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) {
+ foundU = true;
+ vertex->setCoupling(sqr(Emax),getParticleData(2),getParticleData(-2),Mediator_);
+ cSMmed_[1] = vertex->norm();
+ }
+ else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) {
+ foundS = true;
+ vertex->setCoupling(sqr(Emax),getParticleData(3),getParticleData(-3),Mediator_);
+ cSMmed_[2] = vertex->norm();
+ }
+ }
+ }
+ }
+ if(!foundDM) {
+ throw InitException() << "Cannot find DM coupling in MEDM2Mesons::doinit()";
+ }
+ if(!foundD) {
+ throw InitException() << "Cannot find down quark coupling in MEDM2Mesons::doinit()";
+ }
+ if(!foundU) {
+ throw InitException() << "Cannot find up quark coupling in MEDM2Mesons::doinit()";
+ }
+ if(!foundS) {
+ throw InitException() << "Cannot find strange quark coupling in MEDM2Mesons::doinit()";
+ }
+ cDMmed_ = 1.;
+ cSMmed_ = {1.0,1.0,1.0};
MEMultiChannel::doinit();
}
void MEDM2Mesons::persistentOutput(PersistentOStream & os) const {
- os << current_ << modeMap_ << incomingA_ << incomingB_ << Mediator_ << cDMmed_ << cSMmed_ << ounit(MMed_,GeV);
+ os << current_ << modeMap_ << incomingA_ << incomingB_ << Mediator_ << cDMmed_ << cSMmed_;
}
void MEDM2Mesons::persistentInput(PersistentIStream & is, int) {
- is >> current_ >> modeMap_ >> incomingA_ >> incomingB_ >> Mediator_ >> cDMmed_ >> cSMmed_ >> iunit(MMed_,GeV);
+ is >> current_ >> modeMap_ >> incomingA_ >> incomingB_ >> Mediator_ >> cDMmed_ >> cSMmed_;
}
//The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEDM2Mesons,MEMultiChannel>
describeHerwigMEDM2Mesons("Herwig::MEDM2Mesons", "Herwig.so");
void MEDM2Mesons::Init() {
static ClassDocumentation<MEDM2Mesons> documentation
("The MEDM2Mesons class simulates the annhilation of"
" DM particles to mesons at low energy");
static Reference<MEDM2Mesons,WeakCurrent> interfaceWeakCurrent
("WeakCurrent",
"The reference for the decay current to be used.",
&MEDM2Mesons::current_, false, false, true, false, false);
static Reference<MEDM2Mesons,ParticleData> interfaceIncomingA
("IncomingA",
"First incoming particle",
&MEDM2Mesons::incomingA_, false, false, true, false, false);
static Reference<MEDM2Mesons,ParticleData> interfaceIncomingB
("IncomingB",
"Second incoming particle",
&MEDM2Mesons::incomingB_, false, false, true, false, false);
static Reference<MEDM2Mesons,ParticleData> interfaceMediator_
("Mediator",
"DM mediator",
&MEDM2Mesons::Mediator_, false, false, true, false, false);
- static Parameter<MEDM2Mesons,double> interfacecDMmed
- ("cDMmed",
- "coupling of DM to dark mediator",
- &MEDM2Mesons::cDMmed_, 1.0, 0., 10., false, false, Interface::limited);
-
- static ParVector<MEDM2Mesons,double> interfacecSMmed
- ("cSMmed",
- "coupling of SM to dark mediator",
- &MEDM2Mesons::cSMmed_, -1 , 1.0 , -10. , 10. , false, false, Interface::limited);
-
}
double MEDM2Mesons::me2(const int ichan) const {
// compute the incoming current
LorentzPolarizationVectorInvE lepton[2][2];
if(incomingA_->iSpin()==PDT::Spin1Half && incomingB_->iSpin()==PDT::Spin1Half) {
SpinorWaveFunction em_in( meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction ep_in( meMomenta()[1],mePartonData()[1],incoming);
vector<SpinorWaveFunction> f1;
vector<SpinorBarWaveFunction> a1;
for(unsigned int ix=0;ix<2;++ix) {
em_in.reset(ix);
f1.push_back(em_in);
ep_in.reset(ix);
a1.push_back(ep_in);
}
// this should be coupling of DM to mediator/ mediator propagator
complex<Energy> mmed = Mediator_->mass();
complex<Energy2> mmed2 = sqr(mmed);
complex<Energy> mwid = Mediator_->width();
complex<Energy2> prop = sHat()-mmed2+Complex(0.,1.)*mmed*mwid;
complex<InvEnergy2> pre = cDMmed_/prop;
for(unsigned ix=0;ix<2;++ix) {
for(unsigned iy=0;iy<2;++iy) {
lepton[ix][iy]= pre*f1[ix].dimensionedWave().vectorCurrent(a1[iy].dimensionedWave());
}
}
}
// TODO think about other spins for the DM
else
assert(false);
// work out the mapping for the hadron vector
int nOut = int(meMomenta().size())-2;
vector<unsigned int> constants(nOut+1);
vector<PDT::Spin > iSpin(nOut);
vector<int> hadrons(nOut);
int itemp(1);
int ix(nOut);
do {
--ix;
iSpin[ix] = mePartonData()[ix+2]->iSpin();
itemp *= iSpin[ix];
constants[ix] = itemp;
hadrons[ix] = mePartonData()[ix+2]->id();
}
while(ix>0);
constants[nOut] = 1;
// calculate the matrix element
me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,iSpin));
// calculate the hadron current
unsigned int imode = current_->decayMode(hadrons);
Energy q = sqrt(sHat());
vector<Lorentz5Momentum> momenta(meMomenta() .begin()+2, meMomenta().end());
tPDVector out = mode(modeMap_.at(imode))->outgoing();
if(ichan<0) iMode(modeMap_.at(imode));
// get the hadronic currents for the I=1 and I=0 components
vector<LorentzPolarizationVectorE>
hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero),
imode,ichan,q,out,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero),
imode,ichan,q,out,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar),
imode,ichan,q,out,momenta,DecayIntegrator::Calculate));
// compute the matrix element
vector<unsigned int> ihel(meMomenta().size());
double output(0.);
unsigned int hI0_size = hadronI0.size();
unsigned int hI1_size = hadronI1.size();
unsigned int hss_size = hadronssbar.size();
unsigned int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size);
for(unsigned int hhel=0;hhel<maxsize;++hhel) {
// map the index for the hadrons to a helicity state
for(int ix=nOut;ix>0;--ix) {
ihel[ix+1]=(hhel%constants[ix-1])/constants[ix];
}
// loop over the helicities of the incoming particles
for(ihel[1]=0;ihel[1]<2;++ihel[1]){
for(ihel[0]=0;ihel[0]<2;++ihel[0]) {
Complex amp = 0.;
// work on coefficients for the I1 and I0 bits
if(hI0_size != 0 )
amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel]));
if(hI1_size !=0)
amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI1[hhel]));
if(hss_size !=0)
amp += cSMmed_[2]* (lepton[ihel[0]][ihel[1]].dot(hadronssbar[hhel]));
me_(ihel)= amp;
output += std::norm(amp);
}
}
}
// symmetry factors
map<long,int> ncount;
double symmetry(1.);
for(tPDPtr o : out) ncount[o->id()]+=1;
for(map<long,int>::const_iterator it=ncount.begin();it!=ncount.end();++it) {
symmetry *= it->second;
}
// prefactors
output *= 0.25*sqr(pow(sqrt(sHat())/q,int(momenta.size()-2)));
return output/symmetry;
}
void MEDM2Mesons::constructVertex(tSubProPtr) {
}
diff --git a/Models/DarkMatter/MEDM2Mesons.h b/Models/DarkMatter/MEDM2Mesons.h
--- a/Models/DarkMatter/MEDM2Mesons.h
+++ b/Models/DarkMatter/MEDM2Mesons.h
@@ -1,182 +1,177 @@
// -*- C++ -*-
#ifndef Herwig_MEDM2Mesons_H
#define Herwig_MEDM2Mesons_H
//
// This is the declaration of the MEDM2Mesons class.
//
#include "Herwig/MatrixElement/MEMultiChannel.h"
#include "Herwig/Decay/WeakCurrents/WeakCurrent.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MEDM2Mesons class implements dark matter annhilation via a vector current to
* mesons at low energies
*
* @see \ref MEDM2MesonsInterfaces "The interfaces"
* defined for MEDM2Mesons.
*/
class MEDM2Mesons: public MEMultiChannel {
public:
/**
* The default constructor.
*/
MEDM2Mesons();
public:
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* Return the scale associated with the last set phase space point.
*/
virtual Energy2 scale() const;
//@}
/**
* Construct the vertex of spin correlations.
*/
virtual void constructVertex(tSubProPtr);
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:
/**
* Return the matrix element squared for a given mode and phase-space channel.
* @param ichan The channel we are calculating the matrix element for.
*/
virtual double me2(const int ichan) const;
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;
//@}
private:
/**
* 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.
*/
MEDM2Mesons & operator=(const MEDM2Mesons &);
private :
/**
* Hadronic current etc
*/
//@{
/**
* the hadronic current
*/
WeakCurrentPtr current_;
/**
* The matrix element
*/
mutable ProductionMatrixElement me_;
/**
* Map for the modes
*/
map<int,int> modeMap_;
//@}
/**
* DM
*/
//@{
/**
* Incoming Particles
*/
PDPtr incomingA_, incomingB_;
//@}
/**
* DM coupling to the dark mediator
*/
- double cDMmed_;
+ Complex cDMmed_;
/**
* SM couplings to the dark mediator
*/
- vector<double> cSMmed_;
-
- /**
- * DM mediator mass
- */
- Energy MMed_;
+ vector<Complex> cSMmed_;
/**
* DM vector mediator
*/
PDPtr Mediator_;
};
}
#endif /* Herwig_MEDM2Mesons_H */
diff --git a/src/DM.model b/src/DM.model
--- a/src/DM.model
+++ b/src/DM.model
@@ -1,153 +1,148 @@
# -*- ThePEG-repository -*-
##################################################
# Common setup for Zprime models
#
# This file does not contain anything that
# users need to touch. User settings are in
# DM.in
#
###################################################
#
# Particle Data object for the new resonances
#
###################################################
cd /Herwig/Particles
# Dark model
create ThePEG::ParticleData Zp
setup Zp 32 Zp 2.0 0.0001 0.0 0.0 0 0 3 0
create ThePEG::BeamParticleData chi
setup chi 52 chi 1.0 0.0 0.0 0.0 0 0 2 0
###################################################
#
# Main directory and model object
#
###################################################
mkdir /Herwig/NewPhysics/DM
cd /Herwig/NewPhysics/DM
create Herwig::DMModel Model HwDMModel.so
# SM couplings
set Model:QCD/RunningAlphaS /Herwig/Couplings/NLOAlphaS
set Model:EW/RunningAlphaEM /Herwig/Couplings/AlphaEM
set Model:EW/CKM /Herwig/CKM
set Model:RunningMass /Herwig/RunningMass
# Z prime couplings
# quark-anti-quark
set Model:cDMmed 1.
set Model:cSMmed 0 1.0
set Model:cSMmed 1 1.0
set Model:cSMmed 2 1.0
###################################################
#
# Vertices
#
###################################################
# create model vertices
mkdir /Herwig/Vertices/DM
cd /Herwig/Vertices/DM
library HwDMModel.so
create Herwig::DMMediatorQuarksVertex DM_QQZPVertex
create Herwig::DMDMMediatorVertex DM_DMDMZPVertex
cd /Herwig/NewPhysics/DM
# SM vertices
set Model:Vertex/FFZ /Herwig/Vertices/FFZVertex
set Model:Vertex/FFW /Herwig/Vertices/FFWVertex
set Model:Vertex/FFH /Herwig/Vertices/FFHVertex
set Model:Vertex/FFG /Herwig/Vertices/FFGVertex
set Model:Vertex/FFP /Herwig/Vertices/FFPVertex
set Model:Vertex/GGG /Herwig/Vertices/GGGVertex
set Model:Vertex/GGGG /Herwig/Vertices/GGGGVertex
set Model:Vertex/WWH /Herwig/Vertices/WWHVertex
set Model:Vertex/WWW /Herwig/Vertices/WWWVertex
set Model:Vertex/WWWW /Herwig/Vertices/WWWWVertex
set Model:Vertex/HGG /Herwig/Vertices/HGGVertex
set Model:Vertex/HHH /Herwig/Vertices/HHHVertex
set Model:Vertex/WWHH /Herwig/Vertices/WWHHVertex
set Model:Vertex/HHH /Herwig/Vertices/HHHVertex
set Model:Vertex/HPP /Herwig/Vertices/HPPVertex
# model vertices
set Model:Vertex/QQZpVertex /Herwig/Vertices/DM/DM_QQZPVertex
set Model:Vertex/DMDMZpVertex /Herwig/Vertices/DM/DM_DMDMZPVertex
###################################################
#
# Set up spin correlation Decayers
#
###################################################
cd /Herwig/NewPhysics
set TwoBodyDC:CreateDecayModes Yes
set ThreeBodyDC:CreateDecayModes No
# which particles get the off-shell treatment
set NewModel:WhichOffshell All
# particles for which decays are included
#insert NewModel:DecayParticles 0 /Herwig/Particles/Zp
###################################################
# Set up the model framework
###################################################
set DM/Model:ModelGenerator NewModel
###################################################
#
# Choose Model over SM
#
###################################################
cd /Herwig/Generators
set EventGenerator:StandardModelParameters /Herwig/NewPhysics/DM/Model
###################################################
# Matrix Elements for the low energy annhilation
###################################################
cd /Herwig/MatrixElements
# X X -> mesons
create Herwig::MEDM2Mesons DMMatrixElement
-set DMMatrixElement:IncomingA /Herwig/Particles/chi
+set DMMatrixElement:IncomingA /Herwig/Particles/chi
set DMMatrixElement:IncomingB /Herwig/Particles/chi
-# mediator coupling to u,d,s
-set DMMatrixElement:cSMmed[0] 1.0
-set DMMatrixElement:cSMmed[1] 1.0
-set DMMatrixElement:cSMmed[2] 1.0
-# DM to mediator coupling
-set DMMatrixElement:cDMmed 1.0
+# DM to mediator
set DMMatrixElement:Mediator /Herwig/Particles/Zp
# create the matrix elements for the different processes
cp DMMatrixElement MEDM2Pions
set MEDM2Pions:WeakCurrent /Herwig/Decays/TwoPionCzyzCurrent
cp DMMatrixElement MEDM2Kaons
set MEDM2Kaons:WeakCurrent /Herwig/Decays/TwoKaonCzyzCurrent
cp DMMatrixElement MEDM3Pions
set MEDM3Pions:WeakCurrent /Herwig/Decays/ThreePionCzyzCurrent
cp DMMatrixElement MEDM4Pions
set MEDM4Pions:WeakCurrent /Herwig/Decays/FourPionCzyzCurrent
cp DMMatrixElement MEDM2EtaPiPi
set MEDM2EtaPiPi:WeakCurrent /Herwig/Decays/EtaPiPiCurrent
cp DMMatrixElement MEDM2EtaPrimePiPi
set MEDM2EtaPrimePiPi:WeakCurrent /Herwig/Decays/EtaPrimePiPiCurrent
cp DMMatrixElement MEDM2OmegaPiPi
set MEDM2OmegaPiPi:WeakCurrent /Herwig/Decays/OmegaPiPiCurrent
cp DMMatrixElement MEDM2OmegaPi
set MEDM2OmegaPi:WeakCurrent /Herwig/Decays/OmegaPiCurrent
cp DMMatrixElement MEDM2PiGamma
set MEDM2PiGamma:WeakCurrent /Herwig/Decays/PiGammaCurrent
cp DMMatrixElement MEDM2EtaPhoton
set MEDM2EtaPhoton:WeakCurrent /Herwig/Decays/EtaGammaCurrent
cp DMMatrixElement MEDM2EtaPhi
set MEDM2EtaPhi:WeakCurrent /Herwig/Decays/EtaPhiCurrent
cp DMMatrixElement MEDM2EtaOmega
set MEDM2EtaOmega:WeakCurrent /Herwig/Decays/EtaOmegaCurrent
cp DMMatrixElement MEDM2ppbar
set /Herwig/Decays/CzyzCurrent:FormFactor /Herwig/Decays/CzyzFormFactor
set MEDM2ppbar:WeakCurrent /Herwig/Decays/CzyzCurrent
cp DMMatrixElement MEDM2KKPi
set MEDM2KKPi:WeakCurrent /Herwig/Decays/KKPiCurrent
cp DMMatrixElement MEDM2PhiPi
set MEDM2PhiPi:WeakCurrent /Herwig/Decays/PhiPiCurrent