diff --git a/Decay/General/Makefile.am b/Decay/General/Makefile.am --- a/Decay/General/Makefile.am +++ b/Decay/General/Makefile.am @@ -1,75 +1,75 @@ noinst_LTLIBRARIES = libHwGeneralDecay.la libHwGeneralDecay_la_SOURCES = \ GeneralTwoBodyDecayer.h GeneralTwoBodyDecayer.fh GeneralTwoBodyDecayer.cc \ GeneralThreeBodyDecayer.h GeneralThreeBodyDecayer.fh \ GeneralThreeBodyDecayer.cc \ GeneralCurrentDecayer.h GeneralCurrentDecayer.fh GeneralCurrentDecayer.cc \ -DMMediatorDecayer.h DMMediatorDecayer.fh DMMediatorDecayer.cc \ +VectorCurrentDecayer.h VectorCurrentDecayer.fh VectorCurrentDecayer.cc \ GeneralFourBodyDecayer.h GeneralFourBodyDecayer.fh \ GeneralFourBodyDecayer.cc \ StoFFFFDecayer.h StoFFFFDecayer.cc # check the gauge invariance of corrections (for testing ONLY) #libHwGeneralDecay_la_CPPFLAGS= $(AM_CPPFLAGS) -DGAUGE_CHECK nodist_libHwGeneralDecay_la_SOURCES = \ GeneralDecayer__all.cc BUILT_SOURCES = GeneralDecayer__all.cc CLEANFILES = GeneralDecayer__all.cc GeneralDecayer__all.cc : $(DIR_H_FILES) $(DIR_CC_FILES) Makefile @echo "Concatenating .cc files into $@" @$(top_srcdir)/cat_with_cpplines $(DIR_CC_FILES) > $@ EXTRA_DIST = $(ALL_H_FILES) $(ALL_CC_FILES) DIR_H_FILES = $(addprefix $(srcdir)/,$(ALL_H_FILES)) ALL_H_FILES = \ FFSDecayer.h \ FFVDecayer.h \ SFFDecayer.h \ SSSDecayer.h \ SSVDecayer.h \ SVVDecayer.h \ TFFDecayer.h \ TSSDecayer.h \ TVVDecayer.h \ VFFDecayer.h \ VSSDecayer.h \ VVSDecayer.h \ VVVDecayer.h \ SRFDecayer.h \ FRSDecayer.h \ FRVDecayer.h \ FtoFFFDecayer.h \ StoSFFDecayer.h \ StoFFVDecayer.h \ VtoFFVDecayer.h \ FtoFVVDecayer.h \ FFVCurrentDecayer.h DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES)) ALL_CC_FILES = \ FFSDecayer.cc \ FFVDecayer.cc \ SFFDecayer.cc \ SSSDecayer.cc \ SSVDecayer.cc \ SVVDecayer.cc \ TFFDecayer.cc \ TSSDecayer.cc \ TVVDecayer.cc \ VFFDecayer.cc \ VSSDecayer.cc \ VVSDecayer.cc \ VVVDecayer.cc \ SRFDecayer.cc \ FRSDecayer.cc \ FRVDecayer.cc \ FtoFFFDecayer.cc \ StoSFFDecayer.cc \ StoFFVDecayer.cc \ VtoFFVDecayer.cc \ FtoFVVDecayer.cc \ FFVCurrentDecayer.cc diff --git a/Decay/General/DMMediatorDecayer.cc b/Decay/General/VectorCurrentDecayer.cc rename from Decay/General/DMMediatorDecayer.cc rename to Decay/General/VectorCurrentDecayer.cc --- a/Decay/General/DMMediatorDecayer.cc +++ b/Decay/General/VectorCurrentDecayer.cc @@ -1,273 +1,273 @@ // -*- C++ -*- // -// DMMediatorDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator +// VectorCurrentDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member -// functions of the DMMediatorDecayer class. +// functions of the VectorCurrentDecayer class. // -#include "DMMediatorDecayer.h" +#include "VectorCurrentDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include "Herwig/Models/General/BSMModel.h" using namespace Herwig; -IBPtr DMMediatorDecayer::clone() const { +IBPtr VectorCurrentDecayer::clone() const { return new_ptr(*this); } -IBPtr DMMediatorDecayer::fullclone() const { +IBPtr VectorCurrentDecayer::fullclone() const { return new_ptr(*this); } -void DMMediatorDecayer::persistentOutput(PersistentOStream & os) const { +void VectorCurrentDecayer::persistentOutput(PersistentOStream & os) const { os << inpart_ << currentOut_ << current_ << mode_ << wgtloc_ << wgtmax_ << weights_ << cSMmed_; } -void DMMediatorDecayer::persistentInput(PersistentIStream & is, int) { +void VectorCurrentDecayer::persistentInput(PersistentIStream & is, int) { is >> inpart_ >> currentOut_ >> current_ >> mode_ >> wgtloc_ >> wgtmax_ >> weights_ >> cSMmed_; } // The following static variable is needed for the type // description system in ThePEG. -DescribeClass -describeHerwigDMMediatorDecayer("Herwig::DMMediatorDecayer", "Herwig.so"); +DescribeClass +describeHerwigVectorCurrentDecayer("Herwig::VectorCurrentDecayer", "Herwig.so"); -void DMMediatorDecayer::Init() { +void VectorCurrentDecayer::Init() { - static ClassDocumentation documentation - ("The DMMediatorDecayer class is designed for the decays of low mass vector bosons"); + static ClassDocumentation documentation + ("The VectorCurrentDecayer class is designed for the decays of low mass vector bosons"); } -void DMMediatorDecayer::setDecayInfo(PDPtr in, const vector & outCurrent, WeakCurrentPtr current) { +void VectorCurrentDecayer::setDecayInfo(PDPtr in, const vector & outCurrent, WeakCurrentPtr current) { inpart_ = in; currentOut_ = outCurrent; current_ = current; // cast the model Ptr::ptr model = dynamic_ptr_cast::ptr>(generator()->standardModel()); bool 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 ext = vertex->search(iloc, in->id()); if(ext.empty()) continue; for(unsigned int ioff=0;ioffsetCoupling(sqr(in->mass()),getParticleData(1),getParticleData(-1),in); 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(in->mass()),getParticleData(2),getParticleData(-2),in); 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(in->mass()),getParticleData(3),getParticleData(-3),in); cSMmed_[2] = vertex->norm(); } } } } if(!foundD) { - throw InitException() << "Cannot find down quark coupling in DMMediatorDecayer::doinit()"; + throw InitException() << "Cannot find down quark coupling in VectorCurrentDecayer::doinit()"; } if(!foundU) { - throw InitException() << "Cannot find up quark coupling in DMMediatorDecayer::doinit()"; + throw InitException() << "Cannot find up quark coupling in VectorCurrentDecayer::doinit()"; } if(!foundS) { - throw InitException() << "Cannot find strange quark coupling in DMMediatorDecayer::doinit()"; + throw InitException() << "Cannot find strange quark coupling in VectorCurrentDecayer::doinit()"; } } -int DMMediatorDecayer::modeNumber(bool & cc, tcPDPtr parent, +int VectorCurrentDecayer::modeNumber(bool & cc, tcPDPtr parent, const tPDVector & children) const { vector id; id.push_back(parent->id()); for(unsigned int ix=0;ixid()); return modeNumber(cc,id); } -void DMMediatorDecayer::doinitrun() { +void VectorCurrentDecayer::doinitrun() { current_->initrun(); DecayIntegrator::doinitrun(); } -void DMMediatorDecayer::doinit() { +void VectorCurrentDecayer::doinit() { DecayIntegrator::doinit(); // make sure the current got initialised current_->init(); // find the mode for(unsigned int ix=0;ixnumberOfModes();++ix) { // get the external particles for this mode int iq(0),ia(0); tPDVector ptemp = current_->particles(inpart_->iCharge(),ix,iq,ia); // check this is the right mode if(ptemp.size()!=currentOut_.size()) continue; vector matched(ptemp.size(),false); bool match = true; for(unsigned int iy=0;iycreateMode(inpart_->iCharge(),tcPDPtr(),FlavourInfo(), ix,mode,0,-1,channel,inpart_->mass()); if(done) { // the maximum weight and the channel weights // the weights for the channel if(weights_.empty()) { weights_.resize(mode->channels().size(),1./(mode->channels().size())); } mode_ = ix; // special for the two body modes if(out.size()==2) { weights_.clear(); mode=new_ptr(PhaseSpaceMode(inpart_,out,1.)); } mode->maxWeight(wgtmax_); mode->setWeights(weights_); addMode(mode); } break; } } -int DMMediatorDecayer::modeNumber(bool & cc, vector id) const { +int VectorCurrentDecayer::modeNumber(bool & cc, vector id) const { // incoming particle long idtemp; tPDPtr p0=getParticleData(id[0]); idtemp = p0->CC() ? -id[0] : id[0]; if( id[0] ==inpart_->id()) cc=false; else if(idtemp==inpart_->id()) cc=true ; else return -1; vector idout; for(vector::iterator it=++id.begin();it!=id.end();++it) { idout.push_back(*it); } unsigned int icurr=current_->decayMode(idout); if(mode_==icurr) return 0; else return -1; } -void DMMediatorDecayer:: +void VectorCurrentDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast(&part), Helicity::incoming,true,false); weakCurrent()->constructSpinInfo(ParticleVector(decay.begin(),decay.end())); } -double DMMediatorDecayer::me2(const int ichan, const Particle & part, +double VectorCurrentDecayer::me2(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { using namespace ThePEG::Helicity; // polarization vectors for the incoming particle if(meopt==Initialize) { VectorWaveFunction::calculateWaveFunctions(vectors_,rho_, const_ptr_cast(&part), incoming,false); // fix rho if no correlations fixRho(rho_); } // work out the mapping for the hadron vector int nOut = momenta.size(); vector constants(nOut+1); vector iSpin(nOut); vector hadrons(nOut); int itemp(1); int ix(nOut); do { --ix; iSpin[ix] = outgoing[ix]->iSpin(); itemp *= iSpin[ix]; constants[ix] = itemp; hadrons[ix] = outgoing[ix]->id(); } while(ix>0); constants[nOut] = 1; Energy2 scale(sqr(part.mass())); // calculate the hadron current Energy q = part.mass(); // currents for the different flavour components vector hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); vector hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); vector hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); // compute the matrix element GeneralDecayMEPtr newME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,iSpin))); vector ihel(momenta.size()+1); 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;hhel0;--ix) { ihel[ix]=(hhel%constants[ix-1])/constants[ix]; } for(ihel[0]=0;ihel[0]<3;++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.)/q*(vectors_[ihel[0]].wave().dot(hadronI0[hhel])); if(hI1_size !=0) amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI1[hhel])); if(hss_size !=0) amp += cSMmed_[2] /q*(vectors_[ihel[0]].wave().dot(hadronssbar[hhel])); (*newME)(ihel) = amp; } } // store the matrix element ME(newME); // return the answer double output = (ME()->contract(rho_)).real(); return output; } -Energy DMMediatorDecayer::partialWidth(tPDPtr part, vector out) { +Energy VectorCurrentDecayer::partialWidth(tPDPtr part, vector out) { vector id; id.push_back(part->id()); for(unsigned int ix=0;ixid()); bool cc; int mode=modeNumber(cc,id); imode(mode); return initializePhaseSpaceMode(mode,true,true); } diff --git a/Decay/General/DMMediatorDecayer.fh b/Decay/General/VectorCurrentDecayer.fh rename from Decay/General/DMMediatorDecayer.fh rename to Decay/General/VectorCurrentDecayer.fh --- a/Decay/General/DMMediatorDecayer.fh +++ b/Decay/General/VectorCurrentDecayer.fh @@ -1,22 +1,22 @@ // -*- C++ -*- // -// This is the forward declaration of the DMMediatorDecayer class. +// This is the forward declaration of the VectorCurrentDecayer class. // -#ifndef HERWIG_DMMediatorDecayer_FH -#define HERWIG_DMMediatorDecayer_FH +#ifndef HERWIG_VectorCurrentDecayer_FH +#define HERWIG_VectorCurrentDecayer_FH #include "ThePEG/Config/Pointers.h" namespace Herwig { -class DMMediatorDecayer; +class VectorCurrentDecayer; } namespace ThePEG { -ThePEG_DECLARE_POINTERS(Herwig::DMMediatorDecayer,DMMediatorDecayerPtr); +ThePEG_DECLARE_POINTERS(Herwig::VectorCurrentDecayer,VectorCurrentDecayerPtr); } #endif diff --git a/Decay/General/DMMediatorDecayer.h b/Decay/General/VectorCurrentDecayer.h rename from Decay/General/DMMediatorDecayer.h rename to Decay/General/VectorCurrentDecayer.h --- a/Decay/General/DMMediatorDecayer.h +++ b/Decay/General/VectorCurrentDecayer.h @@ -1,238 +1,238 @@ // -*- C++ -*- // -// DMMediatorDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator +// VectorCurrentDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // -#ifndef HERWIG_DMMediatorDecayer_H -#define HERWIG_DMMediatorDecayer_H +#ifndef HERWIG_VectorCurrentDecayer_H +#define HERWIG_VectorCurrentDecayer_H // -// This is the declaration of the DMMediatorDecayer class. +// This is the declaration of the VectorCurrentDecayer class. // #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/WeakCurrents/WeakCurrent.h" #include "Herwig/Decay/PhaseSpaceMode.h" -#include "DMMediatorDecayer.fh" +#include "VectorCurrentDecayer.fh" namespace Herwig { using namespace ThePEG; using ThePEG::Helicity::VectorWaveFunction; /** - * Here is the documentation of the DMMediatorDecayer class. + * Here is the documentation of the VectorCurrentDecayer class. * - * @see \ref DMMediatorDecayerInterfaces "The interfaces" - * defined for DMMediatorDecayer. + * @see \ref VectorCurrentDecayerInterfaces "The interfaces" + * defined for VectorCurrentDecayer. */ -class DMMediatorDecayer: public DecayIntegrator { +class VectorCurrentDecayer: public DecayIntegrator { public: /** * The default constructor. */ - DMMediatorDecayer() : cSMmed_({0.,0.,0.}), wgtmax_(0.) + VectorCurrentDecayer() : cSMmed_({0.,0.,0.}), wgtmax_(0.) {} /** @name Virtual functions required by the Decayer class. */ //@{ /** * Which of the possible decays is required * @param cc Is this mode the charge conjugate * @param parent The decaying particle * @param children The decay products */ virtual int modeNumber(bool & cc, tcPDPtr parent,const tPDVector & children) const; //@} /** @name Virtual functions required by the Decayer class. */ //@{ /** * Return the matrix element squared for a given mode and phase-space channel. * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param outgoing The particles produced in the decay * @param momenta The momenta of the particles produced in the decay * @param meopt Option for the calculation of the matrix element * @return The matrix element squared for the phase-space configuration. */ double me2(const int ichan,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const; /** * Construct the SpinInfos for the particles produced in the decay */ virtual void constructSpinInfo(const Particle & part, ParticleVector outgoing) const; /** * Function to return partial Width * @param inpart Pointer to incoming particle data object * @param out The outgoing particles from the current */ virtual Energy partialWidth(tPDPtr inpart, vector out); //@} /** * set up the decay */ void setDecayInfo(PDPtr in, const vector & outCurrent, WeakCurrentPtr current); 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(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); //@} protected: /** * The number of the mode * @param cc Whether of not this is the charge conjugate of the defined mode * @param id The PDG codes of the particles */ int modeNumber(bool & cc, vector id) const; /** * Access to the map between the number of the mode and the modes in * the current */ unsigned int mode() const { return mode_; } /** * Access to the weak current */ WeakCurrentPtr weakCurrent() const { return current_; } private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ - DMMediatorDecayer & operator=(const DMMediatorDecayer &) = delete; + VectorCurrentDecayer & operator=(const VectorCurrentDecayer &) = delete; private: /** * Incoming particle **/ PDPtr inpart_; /** * Outgoing particles from the current */ vector currentOut_; /** * DM coupling to the dark mediator */ Complex cDMmed_; /** * SM couplings to the dark mediator */ vector cSMmed_; /** * Pointer to the current */ WeakCurrentPtr current_; /** * mapping of the modes to the currents */ unsigned int mode_; /** * location of the weights */ int wgtloc_; /** * the maximum weight */ double wgtmax_; /** * The weights for the different channels */ vector weights_; /** * Spin density matrix */ mutable RhoDMatrix rho_; /** * Polarization vectors for the decaying particle */ mutable vector vectors_; }; } -#endif /* HERWIG_DMMediatorDecayer_H */ +#endif /* HERWIG_VectorCurrentDecayer_H */ diff --git a/Models/General/VectorCurrentDecayConstructor.cc b/Models/General/VectorCurrentDecayConstructor.cc --- a/Models/General/VectorCurrentDecayConstructor.cc +++ b/Models/General/VectorCurrentDecayConstructor.cc @@ -1,170 +1,170 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the VectorCurrentDecayConstructor class. // #include "VectorCurrentDecayConstructor.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 "Herwig/Decay/General/DMMediatorDecayer.h" +#include "Herwig/Decay/General/VectorCurrentDecayer.h" namespace { struct ParticleOrdering { /** * Operator for the ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ bool operator() (tcPDPtr p1, tcPDPtr p2) const { return abs(p1->id()) > abs(p2->id()) || ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) || ( p1->id() == p2->id() && p1->fullName() > p2->fullName() ); } }; } using namespace Herwig; IBPtr VectorCurrentDecayConstructor::clone() const { return new_ptr(*this); } IBPtr VectorCurrentDecayConstructor::fullclone() const { return new_ptr(*this); } void VectorCurrentDecayConstructor::persistentOutput(PersistentOStream & os) const { os << ounit(massCut_,GeV) << current_; } void VectorCurrentDecayConstructor::persistentInput(PersistentIStream & is, int) { is >> iunit(massCut_,GeV) >> current_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigVectorCurrentDecayConstructor("Herwig::VectorCurrentDecayConstructor", "Herwig.so"); void VectorCurrentDecayConstructor::Init() { static ClassDocumentation documentation ("The VectorCurrentDecayConstructor class constructs the decays of low mass vector bosons" " to hadrons via the weak current"); static RefVector interfaceCurrent ("Current", "The current for the decay mode", &VectorCurrentDecayConstructor::current_, -1, false, false, true, false, false); static Parameter interfaceMassCut ("MassCut", "The maximum mass difference for the decay", &VectorCurrentDecayConstructor::massCut_, GeV, 2.0*GeV, 1.0*GeV, 10.0*GeV, false, false, Interface::limited); } void VectorCurrentDecayConstructor::doinit() { NBodyDecayConstructorBase::doinit(); model_ = dynamic_ptr_cast::pointer>(generator()->standardModel()); } void VectorCurrentDecayConstructor::DecayList(const set & particles) { if( particles.empty() ) return; unsigned int nv(model_->numberOfVertices()); for(PDPtr part : particles) { if(part->iSpin()!=PDT::Spin1) continue; if(part->iCharge()!=0) continue; bool foundD(false),foundU(false),foundS(false); if(part->mass()>massCut_) continue; for(unsigned int iv = 0; iv < nv; ++iv) { VertexBasePtr vertex = model_->vertex(iv); if( !vertex->isIncoming(part) || vertex->getNpoint() != 3 ) continue; for(unsigned int iloc = 0;iloc < 3; ++iloc) { vector ext = vertex->search(iloc, part->id()); if(ext.empty()) continue; for(unsigned int ioff=0;ioffinit(); for(unsigned int imode=0;imodenumberOfModes();++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; // order the particles bool skip=false; for(unsigned int ix=0;ix outgoing(out.begin(),out.end()); string tag = part->PDGName() + "->"; bool first=false; for(tcPDPtr part : outgoing) { if(!first) first=true; else tag+=","; tag+=part->PDGName(); } tag+=";"; int charge(0); for(tcPDPtr part : outgoing) charge+=part->iCharge(); if(charge!=0) continue; // create the decayer ostringstream fullname; fullname << "/Herwig/Decays/DMMediator_" << part->PDGName(); for(tcPDPtr part : out) fullname << "_" << part->PDGName(); - string classname = "Herwig::DMMediatorDecayer"; - DMMediatorDecayerPtr decayer = dynamic_ptr_cast + string classname = "Herwig::VectorCurrentDecayer"; + VectorCurrentDecayerPtr decayer = dynamic_ptr_cast (generator()->preinitCreate(classname,fullname.str())); decayer->setDecayInfo(part,out,current); // // set decayer options from base class // setDecayerInterfaces(fullname.str()); // initialize the decayer decayer->init(); // calculate the width Energy pWidth = decayer->partialWidth(part,out); if(pWidth<=ZERO) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); continue; } tDMPtr ndm = generator()->preinitCreateDecayMode(tag); generator()->preinitInterface(ndm, "Decayer", "set", decayer->fullName()); part->stable(false); generator()->preinitInterface(ndm, "Active", "set", "Yes"); setBranchingRatio(ndm, pWidth); } } } }