Page MenuHomeHEPForge
Files F10881281

No OneTemporary
Actions

View Options

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<DMMediatorDecayer,DecayIntegrator>
-describeHerwigDMMediatorDecayer("Herwig::DMMediatorDecayer", "Herwig.so");
+DescribeClass<VectorCurrentDecayer,DecayIntegrator>
+describeHerwigVectorCurrentDecayer("Herwig::VectorCurrentDecayer", "Herwig.so");
-void DMMediatorDecayer::Init() {
+void VectorCurrentDecayer::Init() {
- static ClassDocumentation<DMMediatorDecayer> documentation
- ("The DMMediatorDecayer class is designed for the decays of low mass vector bosons");
+ static ClassDocumentation<VectorCurrentDecayer> documentation
+ ("The VectorCurrentDecayer class is designed for the decays of low mass vector bosons");
}
-void DMMediatorDecayer::setDecayInfo(PDPtr in, const vector<tPDPtr> & outCurrent, WeakCurrentPtr current) {
+void VectorCurrentDecayer::setDecayInfo(PDPtr in, const vector<tPDPtr> & outCurrent, WeakCurrentPtr current) {
inpart_ = in;
currentOut_ = outCurrent;
current_ = current;
// cast the model
Ptr<BSMModel>::ptr model = dynamic_ptr_cast<Ptr<BSMModel>::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<long> ext = vertex->search(iloc, in->id());
if(ext.empty()) continue;
for(unsigned int ioff=0;ioff<ext.size();ioff+=3) {
if(iloc!=2) assert(false);
if(abs(ext[ioff])==1 && abs(ext[ioff+1])==1 && ext[ioff]==-ext[ioff+1]) {
foundD = true;
vertex->setCoupling(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<long> id;
id.push_back(parent->id());
for(unsigned int ix=0;ix<children.size();++ix) id.push_back(children[ix]->id());
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;ix<current_->numberOfModes();++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<bool> matched(ptemp.size(),false);
bool match = true;
for(unsigned int iy=0;iy<currentOut_.size();++iy) {
bool found = false;
for(unsigned int iz=0;iz<ptemp.size();++iz) {
if(!matched[iz]&&ptemp[iz]==currentOut_[iy]) {
found = true;
matched[iz] = true;
break;
}
}
if(!found) {
match = false;
break;
}
}
if(!match) continue;
tPDVector out = {};
out.insert(std::end(out), std::begin(ptemp), std::end(ptemp));
// create the mode
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(inpart_,out,1.));
PhaseSpaceChannel channel(mode,true);
bool done=current_->createMode(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<long> id) const {
+int VectorCurrentDecayer::modeNumber(bool & cc, vector<long> 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<int> idout;
for(vector<long>::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<tPPtr>(&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<Lorentz5Momentum> & momenta,
MEOption meopt) const {
using namespace ThePEG::Helicity;
// polarization vectors for the incoming particle
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(vectors_,rho_,
const_ptr_cast<tPPtr>(&part),
incoming,false);
// fix rho if no correlations
fixRho(rho_);
}
// work out the mapping for the hadron vector
int nOut = momenta.size();
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] = 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<LorentzPolarizationVectorE>
hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero),
mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero),
mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
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<unsigned int> 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;hhel<maxsize;++hhel) {
// map the index for the hadrons to a helicity state
for(int ix=nOut;ix>0;--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<tPDPtr> out) {
+Energy VectorCurrentDecayer::partialWidth(tPDPtr part, vector<tPDPtr> out) {
vector<long> id;
id.push_back(part->id());
for(unsigned int ix=0;ix<out.size();++ix) id.push_back(out[ix]->id());
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<Lorentz5Momentum> & 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<tPDPtr> out);
//@}
/**
* set up the decay
*/
void setDecayInfo(PDPtr in, const vector<tPDPtr> & 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<long> 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<tPDPtr> currentOut_;
/**
* DM coupling to the dark mediator
*/
Complex cDMmed_;
/**
* SM couplings to the dark mediator
*/
vector<Complex> 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<double> weights_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Polarization vectors for the decaying particle
*/
mutable vector<VectorWaveFunction> 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<VectorCurrentDecayConstructor,NBodyDecayConstructorBase>
describeHerwigVectorCurrentDecayConstructor("Herwig::VectorCurrentDecayConstructor", "Herwig.so");
void VectorCurrentDecayConstructor::Init() {
static ClassDocumentation<VectorCurrentDecayConstructor> documentation
("The VectorCurrentDecayConstructor class constructs the decays of low mass vector bosons"
" to hadrons via the weak current");
static RefVector<VectorCurrentDecayConstructor,WeakCurrent> interfaceCurrent
("Current",
"The current for the decay mode",
&VectorCurrentDecayConstructor::current_, -1, false, false, true, false, false);
static Parameter<VectorCurrentDecayConstructor,Energy> 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<Ptr<Herwig::StandardModel>::pointer>(generator()->standardModel());
}
void VectorCurrentDecayConstructor::DecayList(const set<PDPtr> & 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<long> ext = vertex->search(iloc, part->id());
if(ext.empty()) continue;
for(unsigned int ioff=0;ioff<ext.size();ioff+=3) {
if(iloc!=2) assert(false);
if(abs(ext[ioff])==1 && abs(ext[ioff+1])==1 && ext[ioff]==-ext[ioff+1]) {
foundD = true;
}
else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) {
foundU = true;
}
else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) {
foundS = true;
}
}
}
}
if(!foundD && !foundU && !foundS) continue;
for(tWeakCurrentPtr current : current_) {
current->init();
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;
// order the particles
bool skip=false;
for(unsigned int ix=0;ix<out.size();++ix) {
if(!out[ix]) {
skip=true;
break;
}
}
if(skip) continue;
multiset<tcPDPtr,ParticleOrdering> 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<DMMediatorDecayerPtr>
+ string classname = "Herwig::VectorCurrentDecayer";
+ VectorCurrentDecayerPtr decayer = dynamic_ptr_cast<VectorCurrentDecayerPtr>
(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);
}
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:06 AM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4971993
Default Alt Text
(27 KB)

Event Timeline