diff --git a/Decay/WeakCurrents/EtaPhiCurrent.cc b/Decay/WeakCurrents/EtaPhiCurrent.cc
--- a/Decay/WeakCurrents/EtaPhiCurrent.cc
+++ b/Decay/WeakCurrents/EtaPhiCurrent.cc
@@ -1,266 +1,266 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the EtaPhiCurrent class.
//
#include "EtaPhiCurrent.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"
using namespace Herwig;
EtaPhiCurrent::EtaPhiCurrent() {
addDecayMode(3,-3);
setInitialModes(3);
// Masses for the resonances
resMasses_ = {1.670*GeV,2.14*GeV};
// widths for the resonances
resWidths_ = {122*MeV,43.5*MeV};
// amplitudes
amp_ = {0.175/GeV,0.00409/GeV};
// phases
phase_ = {0.,2.19};
}
IBPtr EtaPhiCurrent::clone() const {
return new_ptr(*this);
}
IBPtr EtaPhiCurrent::fullclone() const {
return new_ptr(*this);
}
void EtaPhiCurrent::doinit() {
WeakCurrent::doinit();
assert(phase_.size()==amp_.size());
couplings_.clear();
Complex ii(0.,1.);
for(unsigned int ix=0;ix<amp_.size();++ix) {
couplings_.push_back(amp_[ix]*(cos(phase_[ix])+ii*sin(phase_[ix])));
}
}
void EtaPhiCurrent::persistentOutput(PersistentOStream & os) const {
os << ounit(resMasses_,GeV) << ounit(resWidths_,GeV)
<< ounit(amp_,1./GeV) << phase_ << ounit(couplings_,1./GeV);
}
void EtaPhiCurrent::persistentInput(PersistentIStream & is, int) {
is >> iunit(resMasses_,GeV) >> iunit(resWidths_,GeV)
>> iunit(amp_,1./GeV) >> phase_ >> iunit(couplings_,1./GeV);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<EtaPhiCurrent,WeakCurrent>
describeHerwigEtaPhiCurrent("Herwig::EtaPhiCurrent",
"HwWeakCurrents.so");
void EtaPhiCurrent::Init() {
static ClassDocumentation<EtaPhiCurrent> documentation
("The EtaPhiCurrent class implements a current based"
" on the model of SND for eta + phi "
"The current based on the model of \\cite{Achasov:2018ygm}"
" for eta and phi was used.",
"\\bibitem{Achasov:2018ygm}\n"
"M.~N.~Achasov {\\it et al.},\n"
"%``Measurement of the $e^+e^−\\to\\eta K^+K^−$ Cross Section by Means of the SND Detector,''\n"
"Phys.\\ Atom.\\ Nucl.\\ {\\bf 81} (2018) no.2, 205\n"
" [Yad.\\ Fiz.\\ {\\bf 81} (2018) no.2, 195].\n"
"doi:10.1134/S1063778818020023\n"
"%%CITATION = doi:10.1134/S1063778818020023;%%\n");
static ParVector<EtaPhiCurrent,Energy> interfaceResonanceMasses
("ResonanceMasses",
"The masses of the resonances for the form factor",
&EtaPhiCurrent::resMasses_, GeV, 1, 1680*MeV, 0.5*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<EtaPhiCurrent,Energy> interfaceResonanceWidths
("ResonanceWidths",
"The widths of the resonances for the form factor",
- &EtaPhiCurrent::resWidths_, GeV, 1, 150*MeV, 0.5*GeV, 10.0*GeV,
+ &EtaPhiCurrent::resWidths_, GeV, 1, 150*MeV, ZERO, 10.0*GeV,
false, false, Interface::limited);
static ParVector<EtaPhiCurrent,InvEnergy> interfaceAmplitude
("Amplitude",
"The amplitudes of the couplings",
&EtaPhiCurrent::amp_, 0.00115/GeV, 1, 1./GeV, 0.0/GeV, 10/GeV,
false, false, Interface::limited);
static ParVector<EtaPhiCurrent,double> interfacePhase
("Phase",
"The phases of the couplings in radians",
&EtaPhiCurrent::phase_, 1, 0., 0.0, 2.*Constants::pi,
false, false, Interface::limited);
}
// complete the construction of the decay mode for integration
bool EtaPhiCurrent::createMode(int icharge, tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
unsigned int, PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp ) {
// check the charge
if(icharge!=0) return false;
// check the total isospin
if(Itotal!=IsoSpin::IUnknown) {
if(Itotal!=IsoSpin::IZero) return false;
}
// check I_3
if(i3!=IsoSpin::I3Unknown) {
if(i3!=IsoSpin::I3Zero) return false;
}
// check that the mode is are kinematical allowed
Energy min = getParticleData(ParticleID::eta)->mass()+
getParticleData(ParticleID::phi)->massMin();
if(min>upp) return false;
// resonances for the intermediate channels
tPDVector res = {getParticleData(100333)};
// set up the integration channels;
for(unsigned int ix=0;ix<res.size();++ix) {
if(resonance && resonance!=res[ix]) continue;
mode->addChannel((PhaseSpaceChannel(phase),ires,res[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<res.size();++ix) {
mode->resetIntermediate(res[ix],resMasses_[ix],resWidths_[ix]);
}
return true;
}
// the particles produced by the current
tPDVector EtaPhiCurrent::particles(int icharge, unsigned int imode,int,int) {
assert(icharge==0 && imode<=1);
return {getParticleData(ParticleID::eta),getParticleData(ParticleID::phi)};
}
void EtaPhiCurrent::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>
EtaPhiCurrent::current(tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
const int, const int ichan, Energy & scale,
const tPDVector & ,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption) const {
// check the total isospin
if(Itotal!=IsoSpin::IUnknown) {
if(Itotal!=IsoSpin::IZero) return vector<LorentzPolarizationVectorE>();
}
// check I_3
if(i3!=IsoSpin::I3Unknown) {
if(i3!=IsoSpin::I3Zero) return vector<LorentzPolarizationVectorE>();
}
useMe();
// polarization vectors of the photon
vector<LorentzPolarizationVector> temp(3);
for(unsigned int ix=0;ix<3;++ix) {
temp[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing);
}
// total momentum of the system
Lorentz5Momentum q(momenta[0]+momenta[1]);
// overall hadronic mass
q.rescaleMass();
scale=q.mass();
Energy2 q2(q.m2());
unsigned int imin = 0;
unsigned int imax = couplings_.size();
if(ichan>0) {
imin = ichan;
imax = imin+1;
}
if(resonance) {
switch(abs(resonance->id())) {
case 100333:
imin=0;
break;
default:
assert(false);
}
imax=imin+1;
}
// compute the form factor
complex<InvEnergy> formFactor(ZERO);
// loop over the resonances
for(unsigned int ix=imin;ix<imax;++ix) {
Energy2 mR2(sqr(resMasses_[ix]));
// compute the width
Energy width = resWidths_[ix];
formFactor += couplings_[ix]*mR2/(mR2-q2-Complex(0.,1.)*q.mass()*width);
}
// calculate the current
vector<LorentzPolarizationVectorE> ret(3);
for(unsigned int ix=0;ix<3;++ix) {
ret[ix] += formFactor*Helicity::epsilon(q,temp[ix],momenta[1]);
}
return ret;
}
bool EtaPhiCurrent::accept(vector<int> id) {
if(id.size()!=2) return false;
unsigned int neta(0),nphi(0);
for(unsigned int ix=0;ix<id.size();++ix) {
if(abs(id[ix])==ParticleID::eta) ++neta;
else if(id[ix]==ParticleID::phi) ++nphi;
}
return nphi == 1 && neta==1;
}
unsigned int EtaPhiCurrent::decayMode(vector<int>) {
return 0;
}
// output the information for the database
void EtaPhiCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::EtaPhiCurrent " << name()
<< " HwWeakCurrents.so\n";
for(unsigned int ix=0;ix<resMasses_.size();++ix) {
if(ix<1) output << "newdef " << name() << ":ResonanceMasses " << ix
<< " " << resMasses_[ix]/GeV << "\n";
else output << "insert " << name() << ":ResonanceMasses " << ix
<< " " << resMasses_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<resWidths_.size();++ix) {
if(ix<1) output << "newdef " << name() << ":ResonanceWidths " << ix
<< " " << resWidths_[ix]/GeV << "\n";
else output << "insert " << name() << ":ResonanceWidths " << ix
<< " " << resWidths_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<amp_.size();++ix) {
if(ix<1) output << "newdef " << name() << ":Amplitude " << ix
<< " " << amp_[ix]*GeV << "\n";
else output << "insert " << name() << ":Amplitude " << ix
<< " " << amp_[ix]*GeV << "\n";
}
for(unsigned int ix=0;ix<phase_.size();++ix) {
if(ix<1) output << "newdef " << name() << ":Phase " << ix
<< " " << phase_[ix] << "\n";
else output << "insert " << name() << ":Phase " << ix
<< " " << phase_[ix] << "\n";
}
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/WeakCurrents/OmegaPionSNDCurrent.cc b/Decay/WeakCurrents/OmegaPionSNDCurrent.cc
--- a/Decay/WeakCurrents/OmegaPionSNDCurrent.cc
+++ b/Decay/WeakCurrents/OmegaPionSNDCurrent.cc
@@ -1,360 +1,360 @@
// -*- 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,
+ 0.0*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.,
+ &OmegaPionSNDCurrent::phase_, -1, 0., -360., 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;
break;
default:
return false;
}
}
// check that the mode is are kinematical allowed
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>();
break;
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/PionPhotonCurrent.cc b/Decay/WeakCurrents/PionPhotonCurrent.cc
--- a/Decay/WeakCurrents/PionPhotonCurrent.cc
+++ b/Decay/WeakCurrents/PionPhotonCurrent.cc
@@ -1,390 +1,390 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PionPhotonCurrent class.
//
#include "PionPhotonCurrent.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 "ThePEG/Helicity/epsilon.h"
#include "ThePEG/Helicity/HelicityFunctions.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig/Utilities/Kinematics.h"
using namespace Herwig;
PionPhotonCurrent::PionPhotonCurrent() {
// modes handled
addDecayMode(2,-1);
addDecayMode(1,-1);
addDecayMode(2,-2);
setInitialModes(3);
// Masses for the resonances
resMasses_ = {0.77526*GeV,0.78284*GeV,1.01952*GeV,1.45*GeV,1.70*GeV};
// widths for the resonances
resWidths_ = {0.1491 *GeV,0.00868*GeV,0.00421*GeV,0.40*GeV,0.30*GeV};
// amplitudes
amp_ = {0.0426/GeV,0.0434/GeV,0.00303/GeV,0.00523/GeV,ZERO};
// phases
phase_ = {-12.7,0.,158.,180.,0.};
}
IBPtr PionPhotonCurrent::clone() const {
return new_ptr(*this);
}
IBPtr PionPhotonCurrent::fullclone() const {
return new_ptr(*this);
}
void PionPhotonCurrent::doinit() {
WeakCurrent::doinit();
assert(phase_.size()==amp_.size());
couplings_.clear();
Complex ii(0.,1.);
for(unsigned int ix=0;ix<amp_.size();++ix) {
double phi = phase_[ix]/180.*Constants::pi;
couplings_.push_back(amp_[ix]*(cos(phi)+ii*sin(phi)));
}
mpi_ = getParticleData(ParticleID::piplus)->mass();
}
void PionPhotonCurrent::persistentOutput(PersistentOStream & os) const {
os << ounit(resMasses_,GeV) << ounit(resWidths_,GeV)
<< ounit(amp_,1./GeV) << phase_ << ounit(couplings_,1./GeV)
<< ounit(mpi_,GeV);
}
void PionPhotonCurrent::persistentInput(PersistentIStream & is, int) {
is >> iunit(resMasses_,GeV) >> iunit(resWidths_,GeV)
>> iunit(amp_,1./GeV) >> phase_ >> iunit(couplings_,1./GeV)
>> iunit(mpi_,GeV);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<PionPhotonCurrent,WeakCurrent>
describeHerwigPionPhotonCurrent("Herwig::PionPhotonCurrent",
"HwWeakCurrents.so");
void PionPhotonCurrent::Init() {
static ClassDocumentation<PionPhotonCurrent> documentation
("The PionPhotonCurrent class implements a current based"
" on the model of SND for pion+photon",
"The current based on the model of \\cite{Achasov:2016bfr}"
" for pion and photon was used.",
"\\bibitem{Achasov:2016bfr}\n"
"M.~N.~Achasov {\\it et al.} [SND Collaboration],\n"
"%``Study of the reaction $e^+e^- \\to \\pi^0\\gamma$ with the SND detector at the VEPP-2M collider,''\n"
"Phys.\\ Rev.\\ D {\\bf 93} (2016) no.9, 092001\n"
"doi:10.1103/PhysRevD.93.092001\n"
"[arXiv:1601.08061 [hep-ex]].\n"
"%%CITATION = doi:10.1103/PhysRevD.93.092001;%%\n"
"%20 citations counted in INSPIRE as of 23 Aug 2018\n");
static ParVector<PionPhotonCurrent,Energy> interfaceResonanceMasses
("ResonanceMasses",
"The masses of the resonances for the form factor",
&PionPhotonCurrent::resMasses_, GeV, 5, 775.26*MeV, 0.5*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<PionPhotonCurrent,Energy> interfaceResonanceWidths
("ResonanceWidths",
"The widths of the resonances for the form factor",
&PionPhotonCurrent::resWidths_, GeV, 5, 149.1*MeV, 0.5*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<PionPhotonCurrent,InvEnergy> interfaceAmplitude
("Amplitude",
"The amplitudes of the couplings",
&PionPhotonCurrent::amp_, 1./GeV, 5, 1./GeV, 0.0/GeV, 100./GeV,
false, false, Interface::limited);
static ParVector<PionPhotonCurrent,double> interfacePhase
("Phase",
"The phases of the couplings in degrees",
- &PionPhotonCurrent::phase_, 5, 0., 0.0, 360.0,
+ &PionPhotonCurrent::phase_, 5, 0., -360.0, 360.0,
false, false, Interface::limited);
}
// complete the construction of the decay mode for integration
bool PionPhotonCurrent::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(imode==0) {
if(Itotal!=IsoSpin::IOne) return false;
}
else {
if(Itotal!=IsoSpin::IOne &&
Itotal!=IsoSpin::IZero) 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;
break;
default:
return false;
}
}
// check that the mode is are kinematical allowed
Energy min = imode==0 ?
getParticleData(ParticleID::piplus)->mass() :
getParticleData(ParticleID::pi0 )->mass();
if(min>upp) return false;
// resonances for the intermediate channels
tPDVector res;
if(imode==0) {
if(icharge==-3) res.push_back(getParticleData(-213));
else res.push_back(getParticleData( 213));
}
else {
if(Itotal==IsoSpin::IUnknown||Itotal==IsoSpin::IOne)
res.push_back(getParticleData(113));
if(Itotal==IsoSpin::IUnknown||Itotal==IsoSpin::IZero) {
res.push_back(getParticleData( 223));
res.push_back(getParticleData( 333));
res.push_back(getParticleData(100223));
res.push_back(getParticleData( 30223));
}
}
// set up the integration channels;
for(unsigned int ix=0;ix<res.size();++ix) {
if(resonance && resonance!=res[ix]) continue;
mode->addChannel((PhaseSpaceChannel(phase),ires,res[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<res.size();++ix) {
int ires(0);
if(res[ix]->id()==213) ires=1;
else if(res[ix]->id()== 223) ires=2;
else if(res[ix]->id()==100223) ires=3;
else if(res[ix]->id()== 30223) ires=4;
mode->resetIntermediate(res[ix],resMasses_[ires],resWidths_[ires]);
}
return true;
}
// the particles produced by the current
tPDVector PionPhotonCurrent::particles(int icharge, unsigned int imode,int,int) {
tPDVector extpart = {tPDPtr(),
getParticleData(ParticleID::gamma)};
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 PionPhotonCurrent::constructSpinInfo(ParticleVector decay) const {
vector<LorentzPolarizationVector> temp(3);
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) ++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>
PionPhotonCurrent::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(imode==0) {
if(Itotal!=IsoSpin::IOne) return vector<LorentzPolarizationVectorE>();
}
else {
if(Itotal!=IsoSpin::IOne &&
Itotal!=IsoSpin::IZero) 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>();
break;
default:
return vector<LorentzPolarizationVectorE>();
}
}
useMe();
// polarization vectors of the photon
vector<LorentzPolarizationVector> temp(3);
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) ++ix;
temp[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing);
}
// total momentum of the system
Lorentz5Momentum q(momenta[0]+momenta[1]);
// overall hadronic mass
q.rescaleMass();
scale=q.mass();
Energy2 q2(q.m2());
unsigned int imin = 0;
unsigned int imax = imode==0 ? 1 : 5;
if(Itotal==IsoSpin::IOne)
imax = 1;
else if(Itotal==IsoSpin::IZero) {
imin = 1;
}
if(ichan>0) {
if(Itotal==IsoSpin::IZero)
imin = ichan+1;
else
imin = ichan;
imax=imin+1;
}
if(resonance) {
switch(abs(resonance->id())) {
case 113: case 213 :
imin=0;
break;
case 223:
imin=1;
break;
case 333:
imin=2;
break;
case 100223:
imin = 3;
break;
case 30223 :
imin = 4;
break;
default:
assert(false);
}
imax=imin+1;
}
// compute the form factor
complex<InvEnergy> formFactor(ZERO);
// loop over the resonances
for(unsigned int ix=imin;ix<imax;++ix) {
Energy2 mR2(sqr(resMasses_[ix]));
// compute the width
Energy width(ZERO);
// rho
if(ix==0) {
width = resWidths_[0]*mR2/q2*pow(max(double((q2-4.*sqr(mpi_))/(mR2-4.*sqr(mpi_))),0.),1.5);
}
else {
width = resWidths_[ix];
}
formFactor += couplings_[ix]*mR2/(mR2-q2-Complex(0.,1.)*q.mass()*width);
}
// calculate the current
vector<LorentzPolarizationVectorE> ret(3);
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) continue;
ret[ix] += formFactor*Helicity::epsilon(q,temp[ix],momenta[1]);
}
return ret;
}
bool PionPhotonCurrent::accept(vector<int> id) {
if(id.size()!=2) return false;
unsigned int npiplus(0),npi0(0),ngamma(0);
for(unsigned int ix=0;ix<id.size();++ix) {
if(abs(id[ix])==ParticleID::piplus) ++npiplus;
else if(id[ix]==ParticleID::gamma) ++ngamma;
else if(id[ix]==ParticleID::pi0) ++npi0;
}
return ngamma == 1 && (npiplus==1 || npi0==1);
}
unsigned int PionPhotonCurrent::decayMode(vector<int> id) {
int npip(0),npim(0),npi0(0),ngamma(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::gamma) ++ngamma;
}
if((npip==1 || npim == 1) && ngamma==1)
return 0;
else
return 1;
}
// output the information for the database
void PionPhotonCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::PionPhotonCurrent " << name()
<< " HwWeakCurrents.so\n";
for(unsigned int ix=0;ix<resMasses_.size();++ix) {
if(ix<5) output << "newdef " << name() << ":ResonanceMasses " << ix
<< " " << resMasses_[ix]/GeV << "\n";
else output << "insert " << name() << ":ResonanceMasses " << ix
<< " " << resMasses_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<resWidths_.size();++ix) {
if(ix<5) output << "newdef " << name() << ":ResonanceWidths " << ix
<< " " << resWidths_[ix]/GeV << "\n";
else output << "insert " << name() << ":ResonanceWidths " << ix
<< " " << resWidths_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<amp_.size();++ix) {
if(ix<5) output << "newdef " << name() << ":Amplitude " << ix
<< " " << amp_[ix]*GeV << "\n";
else output << "insert " << name() << ":Amplitude " << ix
<< " " << amp_[ix]*GeV << "\n";
}
for(unsigned int ix=0;ix<phase_.size();++ix) {
if(ix<5) output << "newdef " << name() << ":Phase " << ix
<< " " << phase_[ix] << "\n";
else output << "insert " << name() << ":Phase " << ix
<< " " << phase_[ix] << "\n";
}
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}