diff --git a/Decay/WeakCurrents/OmegaPionSNDCurrent.cc b/Decay/WeakCurrents/OmegaPionSNDCurrent.cc
--- a/Decay/WeakCurrents/OmegaPionSNDCurrent.cc
+++ b/Decay/WeakCurrents/OmegaPionSNDCurrent.cc
@@ -1,354 +1,358 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the OmegaPionSNDCurrent class.
//
#include "OmegaPionSNDCurrent.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/Helicity/epsilon.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Utilities/Kinematics.h"
using namespace Herwig;
OmegaPionSNDCurrent::OmegaPionSNDCurrent() {
// modes handled
addDecayMode(2,-1);
addDecayMode(1,-1);
addDecayMode(2,-2);
setInitialModes(3);
// amplitudes for the weights in the current
amp_ = {1.,0.175,0.014};
phase_ = {0.,124.,-63.};
// rho masses and widths
rhoMasses_ = {0.77526*GeV,1.510*GeV,1.720*GeV};
rhoWidths_ = {0.1491 *GeV,0.44 *GeV,0.25 *GeV};
// coupling
gRhoOmegaPi_ = 15.9/GeV;
fRho_ = 4.9583;
}
IBPtr OmegaPionSNDCurrent::clone() const {
return new_ptr(*this);
}
IBPtr OmegaPionSNDCurrent::fullclone() const {
return new_ptr(*this);
}
void OmegaPionSNDCurrent::doinit() {
WeakCurrent::doinit();
assert(phase_.size()==amp_.size());
wgts_.clear();
Complex ii(0.,1.);
for(unsigned int ix=0;ix<amp_.size();++ix) {
double phi = phase_[ix]/180.*Constants::pi;
wgts_.push_back(amp_[ix]*(cos(phi)+ii*sin(phi)));
}
mpi_ = getParticleData(ParticleID::piplus)->mass();
}
void OmegaPionSNDCurrent::persistentOutput(PersistentOStream & os) const {
os << ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV)
<< amp_ << phase_ << wgts_ << fRho_
<< ounit(gRhoOmegaPi_,1./GeV) << ounit(mpi_,GeV);
}
void OmegaPionSNDCurrent::persistentInput(PersistentIStream & is, int) {
is >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV)
>> amp_ >> phase_ >> wgts_ >> fRho_
>> iunit(gRhoOmegaPi_,1./GeV) >> iunit(mpi_,GeV);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<OmegaPionSNDCurrent,WeakCurrent>
describeHerwigOmegaPionSNDCurrent("Herwig::OmegaPionSNDCurrent",
"HwWeakCurrents.so");
void OmegaPionSNDCurrent::Init() {
static ClassDocumentation<OmegaPionSNDCurrent> documentation
("The OmegaPionSNDCurrent class provides a current for omega pi"
" using the model of SND",
"The current based on \\cite{Achasov:2016zvn} for $\\omega\\pi$ was used.\n",
"\\bibitem{Achasov:2016zvn}"
"M.~N.~Achasov {\\it et al.},\n"
"%``Updated measurement of the $e^+e^- \\to \\omega \\pi^0 \\to \\pi^0\\pi^0\\gamma$ cross section with the SND detector,''\n"
"Phys.\\ Rev.\\ D {\\bf 94} (2016) no.11, 112001\n"
"doi:10.1103/PhysRevD.94.112001\n"
"[arXiv:1610.00235 [hep-ex]].\n"
"%%CITATION = doi:10.1103/PhysRevD.94.112001;%%\n"
"%12 citations counted in INSPIRE as of 22 Aug 2018\n");
static ParVector<OmegaPionSNDCurrent,Energy> interfaceRhoMasses
("RhoMasses",
"The masses of the rho mesons",
&OmegaPionSNDCurrent::rhoMasses_, GeV, -1, 775.26*MeV,
0.5*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<OmegaPionSNDCurrent,Energy> interfaceRhoWidths
("RhoWidths",
"The widths of the rho mesons",
&OmegaPionSNDCurrent::rhoWidths_, GeV, -1, 0.1491*GeV,
0.5*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<OmegaPionSNDCurrent,double> interfaceAmplitudes
("Amplitudes",
"THe amplitudes for the different rho resonances",
&OmegaPionSNDCurrent::amp_, -1, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static ParVector<OmegaPionSNDCurrent,double> interfacePhase
("Phase",
"The phases for the different rho resonances in degrees",
&OmegaPionSNDCurrent::phase_, -1, 0., 0.0, 360.,
false, false, Interface::limited);
static Parameter<OmegaPionSNDCurrent,double> interfacefRho
("fRho",
"The coupling of the photon and the rho meson",
&OmegaPionSNDCurrent::fRho_, 4.9583, 0.0, 100.0,
false, false, Interface::limited);
static Parameter<OmegaPionSNDCurrent,InvEnergy> interfacegRhoOmegaPi
("gRhoOmegaPi",
"The coupling rho-omega-pi",
&OmegaPionSNDCurrent::gRhoOmegaPi_, 1./GeV,
15.9/GeV, 0./GeV, 1000./GeV,
false, false, Interface::limited);
}
// complete the construction of the decay mode for integration
bool OmegaPionSNDCurrent::createMode(int icharge, tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
unsigned int imode,PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp ) {
// check the charge
if((abs(icharge)!=3 && imode == 0) ||
( icharge!=0 && imode >= 1))
return false;
// check the total isospin
if(Itotal!=IsoSpin::IUnknown) {
if(Itotal!=IsoSpin::IOne) return false;
}
// check I_3
if(i3!=IsoSpin::I3Unknown) {
switch(i3) {
case IsoSpin::I3Zero:
if(imode!=1) return false;
break;
case IsoSpin::I3One:
if(imode>1 || icharge ==-3) return false;
break;
case IsoSpin::I3MinusOne:
if(imode>1 || icharge ==3) return false;
default:
return false;
}
}
// check that the mode is are kinematical allowed
- Energy min(getParticleData(ParticleID::piplus)->mass()+
- getParticleData(ParticleID::pi0 )->mass());
+ Energy min = getParticleData(ParticleID::omega)->massMin();
+ if(imode==0)
+ min += getParticleData(ParticleID::piplus)->mass();
+ else
+ min += getParticleData(ParticleID::pi0 )->mass();
if(min>upp) return false;
// set up the integration channels;
vector<tPDPtr> rho;
if(icharge==-3)
rho = {getParticleData(-213),getParticleData(-100213),getParticleData(-30213)};
else if(icharge==0)
rho = {getParticleData( 113),getParticleData( 100113),getParticleData( 30113)};
else if(icharge==3)
rho = {getParticleData( 213),getParticleData( 100213),getParticleData( 30213)};
for(unsigned int ix=0;ix<3;++ix) {
if(resonance && resonance!=rho[ix]) continue;
mode->addChannel((PhaseSpaceChannel(phase),ires,rho[ix],
ires+1,iloc+1,ires+1,iloc+2));
}
// reset the masses and widths of the resonances if needed
for(unsigned int ix=0;ix<3;++ix) {
mode->resetIntermediate(rho[ix],rhoMasses_[ix],rhoWidths_[ix]);
}
return true;
}
// the particles produced by the current
tPDVector OmegaPionSNDCurrent::particles(int icharge, unsigned int imode,int,int) {
tPDVector extpart = {tPDPtr(),
getParticleData(ParticleID::omega)};
if(imode==0) {
if(icharge==3) extpart[0] = getParticleData(ParticleID::piplus );
else if(icharge==-3) extpart[0] = getParticleData(ParticleID::piminus);
}
else {
extpart[0] = getParticleData(ParticleID::pi0);
}
return extpart;
}
void OmegaPionSNDCurrent::constructSpinInfo(ParticleVector decay) const {
vector<LorentzPolarizationVector> temp(3);
for(unsigned int ix=0;ix<3;++ix) {
temp[ix] = HelicityFunctions::polarizationVector(-decay[1]->momentum()
,ix,Helicity::outgoing);
}
ScalarWaveFunction::constructSpinInfo(decay[0],outgoing,true);
VectorWaveFunction::constructSpinInfo(temp,decay[1],
outgoing,true,true);
}
// the hadronic currents
vector<LorentzPolarizationVectorE>
OmegaPionSNDCurrent::current(tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
const int imode, const int ichan, Energy & scale,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption) const {
int icharge = outgoing[0]->iCharge()+outgoing[1]->iCharge();
// check the charge
if((abs(icharge)!=3 && imode == 0) ||
( icharge!=0 && imode == 1))
return vector<LorentzPolarizationVectorE>();
// check the total isospin
if(Itotal!=IsoSpin::IUnknown) {
if(Itotal!=IsoSpin::IOne) return vector<LorentzPolarizationVectorE>();
}
// check I_3
if(i3!=IsoSpin::I3Unknown) {
switch(i3) {
case IsoSpin::I3Zero:
if(imode!=1) return vector<LorentzPolarizationVectorE>();
break;
case IsoSpin::I3One:
if(imode>1 || icharge ==-3) return vector<LorentzPolarizationVectorE>();
break;
case IsoSpin::I3MinusOne:
if(imode>1 || icharge ==3) return vector<LorentzPolarizationVectorE>();
default:
return vector<LorentzPolarizationVectorE>();
}
}
useMe();
vector<LorentzPolarizationVector> temp(3);
for(unsigned int ix=0;ix<3;++ix)
temp[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing);
// locate the particles
Lorentz5Momentum q(momenta[0]+momenta[1]);
// overall hadronic mass
q.rescaleMass();
scale=q.mass();
Energy2 q2(q.m2());
// compute the rho width
Energy2 mr2(sqr(rhoMasses_[0]));
Energy grho = rhoWidths_[0]*mr2/q2*pow((q2-4.*sqr(mpi_))/(mr2-4.*sqr(mpi_)),1.5);
Energy qw = Kinematics::pstarTwoBodyDecay(q.mass(),momenta[0].mass(),momenta[1].mass());
grho += pow<3,1>(qw)*sqr(gRhoOmegaPi_)/12./Constants::pi;
unsigned int imin=0, imax = wgts_.size();
if(ichan>0) {
imin = ichan;
imax = ichan+1;
}
if(resonance) {
switch(resonance->id()/1000) {
case 0:
imin = 0;
break;
case 100:
imin = 1;
break;
case 30 :
imin = 2;
break;
default:
assert(false);
}
imax=imin+1;
}
// compute the prefactor
complex<InvEnergy> pre = gRhoOmegaPi_/fRho_;
Complex bw(0.);
for(unsigned int ix=imin;ix<imax;++ix) {
Energy wid = ix==0 ? grho : rhoWidths_[ix];
Energy2 mR2 = sqr(rhoMasses_[ix]);
bw += mR2*wgts_[ix]/(mR2-q2-Complex(0.,1.)*q.mass()*wid);
}
pre *=bw;
vector<LorentzPolarizationVectorE> ret(3);
for(unsigned int ix=0;ix<3;++ix) {
ret[ix] = pre*Helicity::epsilon(q,temp[ix],momenta[1]);
}
+ if(imode==0) pre *=sqrt(2.);
return ret;
}
bool OmegaPionSNDCurrent::accept(vector<int> id) {
if(id.size()!=2){return false;}
unsigned int npiplus(0),npi0(0),nomega(0);
for(unsigned int ix=0;ix<id.size();++ix) {
if(abs(id[ix])==ParticleID::piplus) ++npiplus;
else if(id[ix]==ParticleID::omega) ++nomega;
else if(id[ix]==ParticleID::pi0) ++npi0;
}
return nomega==1 && (npiplus==1||npi0==1);
}
unsigned int OmegaPionSNDCurrent::decayMode(vector<int> id) {
int npip(0),npim(0),npi0(0),nomega(0);
for(unsigned int ix=0;ix<id.size();++ix) {
if(id[ix]==ParticleID::piplus) ++npip;
else if(id[ix]==ParticleID::piminus) ++npim;
else if(id[ix]==ParticleID::pi0) ++npi0;
else if(id[ix]==ParticleID::omega) ++nomega;
}
if((npip==1 || npim == 1) && nomega==1)
return 0;
else
return 1;
}
// output the information for the database
void OmegaPionSNDCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::OmegaPionSNDCurrent " << name()
<< " HwWeakCurrents.so\n";
for(unsigned int ix=0;ix<rhoMasses_.size();++ix) {
if(ix<3) output << "newdef " << name() << ":RhoMasses " << ix
<< " " << rhoMasses_[ix]/GeV << "\n";
else output << "insert " << name() << ":RhoMasses " << ix
<< " " << rhoMasses_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<rhoWidths_.size();++ix) {
if(ix<3) output << "newdef " << name() << ":RhoWidths " << ix
<< " " << rhoWidths_[ix]/GeV << "\n";
else output << "insert " << name() << ":RhoWidths " << ix
<< " " << rhoWidths_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<amp_.size();++ix) {
if(ix<3) output << "newdef " << name() << ":Amplitudes " << ix
<< " " << amp_[ix] << "\n";
else output << "insert " << name() << ":Amplitudes " << ix
<< " " << amp_[ix] << "\n";
}
for(unsigned int ix=0;ix<phase_.size();++ix) {
if(ix<3) output << "newdef " << name() << ":Phases " << ix
<< " " << phase_[ix] << "\n";
else output << "insert " << name() << ":Phases " << ix
<< " " << phase_[ix] << "\n";
}
output << "newdef " << name() << ":fRho " << fRho_ << "\n";
output << "newdef " << name() << ":gRhoOmegaPi " << gRhoOmegaPi_*GeV << "\n";
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/WeakCurrents/OmegaPionSNDCurrent.h b/Decay/WeakCurrents/OmegaPionSNDCurrent.h
--- a/Decay/WeakCurrents/OmegaPionSNDCurrent.h
+++ b/Decay/WeakCurrents/OmegaPionSNDCurrent.h
@@ -1,229 +1,228 @@
// -*- C++ -*-
#ifndef Herwig_OmegaPionSNDCurrent_H
#define Herwig_OmegaPionSNDCurrent_H
//
// This is the declaration of the OmegaPionSNDCurrent class.
//
#include "WeakCurrent.h"
namespace Herwig {
using namespace ThePEG;
/**
* The OmegaPionSNDCurrent class implements the decay current for \f$\pi^\pm\pi^0 \gamma\f$ via
* an intermediate \f$\omega\f$. It inherits from the <code>WeakCurrent</code>
* class and implements the hadronic current.
*
* The model is based on the one from Phys.Rev. D88 (2013) no.5, 054013.
*
* @see \ref OmegaPionSNDCurrentInterfaces "The interfaces"
* defined for OmegaPionSNDCurrent.
*/
class OmegaPionSNDCurrent: public WeakCurrent {
public:
/**
* The default constructor.
*/
OmegaPionSNDCurrent();
-
public:
/** @name Methods for the construction of the phase space integrator. */
//@{
/**
* Complete the construction of the decay mode for integration.classes inheriting
* from this one.
* This method is purely virtual and must be implemented in the classes inheriting
* from WeakCurrent.
* @param icharge The total charge of the outgoing particles in the current.
* @param resonance If specified only include terms with this particle
* @param Itotal If specified the total isospin of the current
* @param I3 If specified the thrid component of isospin
* @param imode The mode in the current being asked for.
* @param mode The phase space mode for the integration
* @param iloc The location of the of the first particle from the current in
* the list of outgoing particles.
* @param ires The location of the first intermediate for the current.
* @param phase The prototype phase space channel for the integration.
* @param upp The maximum possible mass the particles in the current are
* allowed to have.
* @return Whether the current was sucessfully constructed.
*/
virtual bool createMode(int icharge, tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
unsigned int imode,PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp );
/**
* The particles produced by the current. This just returns the pseudoscalar
* meson.
* @param icharge The total charge of the particles in the current.
* @param imode The mode for which the particles are being requested
* @param iq The PDG code for the quark
* @param ia The PDG code for the antiquark
* @return The external particles for the current.
*/
virtual tPDVector particles(int icharge, unsigned int imode, int iq, int ia);
//@}
/**
* Hadronic current. This method is purely virtual and must be implemented in
* all classes inheriting from this one.
* @param resonance If specified only include terms with this particle
* @param Itotal If specified the total isospin of the current
* @param I3 If specified the thrid component of isospin
* @param imode The mode
* @param ichan The phase-space channel the current is needed for.
* @param scale The invariant mass of the particles in the current.
* @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 current.
*/
virtual vector<LorentzPolarizationVectorE>
current(tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
const int imode, const int ichan,Energy & scale,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption meopt) const;
/**
* Construct the SpinInfo for the decay products
*/
virtual void constructSpinInfo(ParticleVector decay) const;
/**
* Accept the decay. Checks the meson against the list
* @param id The id's of the particles in the current.
* @return Can this current have the external particles specified.
*/
virtual bool accept(vector<int> id);
/**
* Return the decay mode number for a given set of particles in the current.
* Checks the meson against the list
* @param id The id's of the particles in the current.
* @return The number of the mode
*/
virtual unsigned int decayMode(vector<int> id);
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
* @param create Whether or not to add a statement creating the object
*/
virtual void dataBaseOutput(ofstream & os,bool header,bool create) const;
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:
/**
* 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.
*/
OmegaPionSNDCurrent & operator=(const OmegaPionSNDCurrent &);
private :
/**
* Masses of the \f$\rho\f$ resonances
*/
vector<Energy> rhoMasses_;
/**
* Widths of the \f$\rho\f$ resonances
*/
vector<Energy> rhoWidths_;
/**
* Ampltitudes for the different rhos in the current
*/
vector<double> amp_;
/**
* Phases for the different rhos in the current
*/
vector<double> phase_;
/**
* Weights of the different rho resonances in the current
*/
vector<Complex> wgts_;
/**
* Coupling of the rho to the photon, \f$f_\rho\f$.
*/
double fRho_;
/**
* Coupling of the rho to the omega and a pion, \f$g_{\rho\omega\pi}\f$.
*/
InvEnergy gRhoOmegaPi_;
/**
* The pion mass
*/
Energy mpi_;
};
}
#endif /* Herwig_OmegaPionSNDCurrent_H */