No OneTemporary
Actions

Size

36 KB

Subscribers

None

View Options

	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,354 @@
	// -- 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.26,0.06};
	- phase_ = {0.,168.,10.};
	+ amp_ = {1.,0.175,0.014};
	+ phase_ = {0.,124.,-63.};
	// rho masses and widths
	- rhoMasses_ = {0.77526GeV,1.491GeV,1.708*GeV};
	- rhoWidths_ = {0.1491 GeV,0.4 GeV,0.25 *GeV};
	+ rhoMasses_ = {0.77526GeV,1.510GeV,1.720*GeV};
	+ rhoWidths_ = {0.1491 GeV,0.44 GeV,0.25 *GeV};
	// coupling
	- gRhoOmegaPi_ = 15.6/GeV;
	+ 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)+iisin(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:2013btb} for $\\omega\\pi$ was used.\n",
	- "\\bibitem{Achasov:2013btb}\n"
	+ "The current based on \\cite{Achasov:2016zvn} for $\\omega\\pi$ was used.\n",
	+ "\\bibitem{Achasov:2016zvn}"
	"M.~N.~Achasov {\\it et al.},\n"
	- "%``Study of $e^+e^- \\to \\omega\\pi^0 \\to \\pi^0\\pi^0\\gamma$ in the energy range $1.05-2.00$ GeV with SND,''\n"
	- "Phys.\\ Rev.\\ D {\\bf 88} (2013) no.5, 054013\n"
	- "doi:10.1103/PhysRevD.88.054013\n"
	- "[arXiv:1303.5198 [hep-ex]].\n"
	- "%%CITATION = doi:10.1103/PhysRevD.88.054013;%%\n"
	- "%50 citations counted in INSPIRE as of 22 Aug 2018\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.5GeV, 10.0GeV,
	false, false, Interface::limited);

	static ParVector<OmegaPionSNDCurrent,Energy> interfaceRhoWidths
	("RhoWidths",
	"The widths of the rho mesons",
	&OmegaPionSNDCurrent::rhoWidths_, GeV, -1, 0.1491*GeV,
	0.5GeV, 10.0GeV,
	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.6/GeV, 0./GeV, 1000./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());
	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/q2pow((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 += mR2wgts_[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]);
	}
	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/TwoPionPhotonSNDCurrent.cc b/Decay/WeakCurrents/TwoPionPhotonSNDCurrent.cc
	--- a/Decay/WeakCurrents/TwoPionPhotonSNDCurrent.cc
	+++ b/Decay/WeakCurrents/TwoPionPhotonSNDCurrent.cc
	@@ -1,419 +1,420 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the TwoPionPhotonSNDCurrent class.
	//

	#include "TwoPionPhotonSNDCurrent.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;

	TwoPionPhotonSNDCurrent::TwoPionPhotonSNDCurrent() {
	// modes handled
	addDecayMode(2,-1);
	addDecayMode(1,-1);
	addDecayMode(2,-2);
	setInitialModes(3);
	// amplitudes for the weights in the current
	- amp_ = {1.,0.26,0.06};
	- phase_ = {0.,168.,10.};
	+ amp_ = {1.,0.175,0.014};
	+ phase_ = {0.,124.,-63.};
	// rho masses and widths
	- rhoMasses_ = {0.77526GeV,1.491GeV,1.708*GeV};
	- rhoWidths_ = {0.1491 GeV,0.4 GeV,0.25 *GeV};
	+ rhoMasses_ = {0.77526GeV,1.510GeV,1.720*GeV};
	+ rhoWidths_ = {0.1491 GeV,0.44 GeV,0.25 *GeV};
	// coupling
	- gRhoOmegaPi_ = 15.6/GeV;
	+ gRhoOmegaPi_ = 15.9/GeV;
	+ fRho_ = 4.9583;
	gGammaOmegaPi_ = 0.695821538653/GeV;
	fRho_ = 4.9583;
	// omega parameters
	omegaMass_ = 782.65*MeV;
	omegaWidth_ = 8.49 *MeV;
	}

	IBPtr TwoPionPhotonSNDCurrent::clone() const {
	return new_ptr(*this);
	}

	IBPtr TwoPionPhotonSNDCurrent::fullclone() const {
	return new_ptr(*this);
	}

	void TwoPionPhotonSNDCurrent::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)+iisin(phi)));
	}
	mpi_ = getParticleData(ParticleID::piplus)->mass();
	}

	void TwoPionPhotonSNDCurrent::persistentOutput(PersistentOStream & os) const {
	os << ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV)
	<< amp_ << phase_ << wgts_ << fRho_
	<< ounit(gRhoOmegaPi_,1./GeV) << ounit(gGammaOmegaPi_,1./GeV)
	<< ounit(omegaMass_,GeV) << ounit(omegaWidth_,GeV) << ounit(mpi_,GeV);
	}

	void TwoPionPhotonSNDCurrent::persistentInput(PersistentIStream & is, int) {
	is >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV)
	>> amp_ >> phase_ >> wgts_ >> fRho_
	>> iunit(gRhoOmegaPi_,1./GeV) >> iunit(gGammaOmegaPi_,1./GeV)
	>> iunit(omegaMass_,GeV) >> iunit(omegaWidth_,GeV) >> iunit(mpi_,GeV);
	}


	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<TwoPionPhotonSNDCurrent,WeakCurrent>
	describeHerwigTwoPionPhotonSNDCurrent("Herwig::TwoPionPhotonSNDCurrent",
	"HwWeakCurrents.so");

	void TwoPionPhotonSNDCurrent::Init() {

	static ClassDocumentation<TwoPionPhotonSNDCurrent> documentation
	("The TwoPionPhotonSNDCurrent class provides the weka current for"
	"pi pi gamma using the model of SND",
	- "The current based on \\cite{Achasov:2013btb} for $\\pi\\pi^0\\gamma$ was used.\n",
	- "\\bibitem{Achasov:2013btb}\n"
	+ "The current based on \\cite{Achasov:2016zvn} for $\\pi\\pi^0\\gamma$ was used.\n",
	+ "\\bibitem{Achasov:2016zvn}"
	"M.~N.~Achasov {\\it et al.},\n"
	- "%``Study of $e^+e^- \\to \\omega\\pi^0 \\to \\pi^0\\pi^0\\gamma$ in the energy range $1.05-2.00$ GeV with SND,''\n"
	- "Phys.\\ Rev.\\ D {\\bf 88} (2013) no.5, 054013\n"
	- "doi:10.1103/PhysRevD.88.054013\n"
	- "[arXiv:1303.5198 [hep-ex]].\n"
	- "%%CITATION = doi:10.1103/PhysRevD.88.054013;%%\n"
	- "%50 citations counted in INSPIRE as of 22 Aug 2018\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<TwoPionPhotonSNDCurrent,Energy> interfaceRhoMasses
	("RhoMasses",
	"The masses of the rho mesons",
	&TwoPionPhotonSNDCurrent::rhoMasses_, GeV, -1, 775.26*MeV,
	0.5GeV, 10.0GeV,
	false, false, Interface::limited);

	static ParVector<TwoPionPhotonSNDCurrent,Energy> interfaceRhoWidths
	("RhoWidths",
	"The widths of the rho mesons",
	&TwoPionPhotonSNDCurrent::rhoWidths_, GeV, -1, 0.1491*GeV,
	0.5GeV, 10.0GeV,
	false, false, Interface::limited);

	static ParVector<TwoPionPhotonSNDCurrent,double> interfaceAmplitudes
	("Amplitudes",
	"THe amplitudes for the different rho resonances",
	&TwoPionPhotonSNDCurrent::amp_, -1, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static ParVector<TwoPionPhotonSNDCurrent,double> interfacePhase
	("Phase",
	"The phases for the different rho resonances in degrees",
	&TwoPionPhotonSNDCurrent::phase_, -1, 0., 0.0, 360.,
	false, false, Interface::limited);

	static Parameter<TwoPionPhotonSNDCurrent,double> interfacefRho
	("fRho",
	"The coupling of the photon and the rho meson",
	&TwoPionPhotonSNDCurrent::fRho_, 4.9583, 0.0, 100.0,
	false, false, Interface::limited);

	static Parameter<TwoPionPhotonSNDCurrent,InvEnergy> interfacegRhoOmegaPi
	("gRhoOmegaPi",
	"The coupling rho-omega-pi",
	&TwoPionPhotonSNDCurrent::gRhoOmegaPi_, 1./GeV,
	- 15.6/GeV, 0./GeV, 1000./GeV,
	+ 15.9/GeV, 0./GeV, 1000./GeV,
	false, false, Interface::limited);

	static Parameter<TwoPionPhotonSNDCurrent,InvEnergy> interfacegGammaOmegaPi
	("gGammaOmegaPi",
	"The coupling gamma-omega-pi",
	&TwoPionPhotonSNDCurrent::gGammaOmegaPi_, 1./GeV,
	0.695821538653/GeV, 0./GeV, 1000./GeV,
	false, false, Interface::limited);

	static Parameter<TwoPionPhotonSNDCurrent,Energy> interfaceOmegaMass
	("OmegaMass",
	"The mass of the omega meson",
	&TwoPionPhotonSNDCurrent::omegaMass_, GeV, 0.78265GeV, 0.0GeV, 10.0*GeV,
	false, false, Interface::limited);

	static Parameter<TwoPionPhotonSNDCurrent,Energy> interfaceOmegaWidth
	("OmegaWidth",
	"The width of the omega meson",
	&TwoPionPhotonSNDCurrent::omegaWidth_, GeV, 8.49MeV, 0.0GeV, 10.0*GeV,
	false, false, Interface::limited);
	}

	// complete the construction of the decay mode for integration
	bool TwoPionPhotonSNDCurrent::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());
	if(min>upp) return false;
	// set up the integration channels;
	tPDPtr omega(getParticleData(ParticleID::omega));
	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,omega,ires+1,iloc+1,
	ires+2,iloc+2,ires+2,iloc+3));
	// channel with the pions exchanged
	if(icharge==0)
	mode->addChannel((PhaseSpaceChannel(phase),ires,rho[ix],
	ires+1,omega,ires+1,iloc+2,
	ires+2,iloc+1,ires+2,iloc+3));
	}
	// 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]);
	}
	// set up the omega masses and widths
	mode->resetIntermediate(omega,omegaMass_,omegaWidth_);
	return true;
	}

	// the particles produced by the current
	tPDVector TwoPionPhotonSNDCurrent::particles(int icharge, unsigned int imode,int,int) {
	tPDVector extpart = {tPDPtr(),
	getParticleData(ParticleID::pi0),
	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 TwoPionPhotonSNDCurrent::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[2]->momentum()
	,ix,Helicity::outgoing);
	}
	for(unsigned int ix=0;ix<2;++ix)
	ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true);
	VectorWaveFunction::constructSpinInfo(temp,decay[2],
	outgoing,true,true);
	}

	// the hadronic currents
	vector<LorentzPolarizationVectorE>
	TwoPionPhotonSNDCurrent::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()+outgoing[2]->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();
	// 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[2],ix,Helicity::outgoing);
	}
	// total momentum of the system
	Lorentz5Momentum q(momenta[0]+momenta[1]+momenta[2]);
	// overall hadronic mass
	q.rescaleMass();
	scale=q.mass();
	Energy2 q2(q.m2());
	unsigned int imin=0, imax = wgts_.size();
	if(ichan>0) {
	if(outgoing[0]!=outgoing[1])
	imin = ichan;
	else
	imin = ichan/2;
	imax = imin+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;
	}
	vector<LorentzPolarizationVectorE> ret(3);
	// need to include exchange of identical particles for the I_3=0 case
	for(int iorder=0;iorder<2;++iorder) {
	Lorentz5Momentum pout(momenta[2]);
	if(outgoing[0]==outgoing[1]) {
	if(ichan>=0&& ichan%2!=iorder) continue;
	}
	else if(iorder==1) continue;
	// add pion momentum
	if(iorder==0) pout += momenta[1];
	else pout += momenta[0];
	// mass of the omega
	pout.rescaleMass();
	Energy2 s2(pout.m2());
	// compute the rho width
	Energy2 mr2(sqr(rhoMasses_[0]));
	Energy grho = rhoWidths_[0]mr2/q2pow(max(double((q2-4.sqr(mpi_))/(mr2-4.sqr(mpi_))),0.),1.5);
	Energy qw = Kinematics::pstarTwoBodyDecay(q.mass(),pout.mass(),mpi_);
	grho += pow<3,1>(qw)*sqr(gRhoOmegaPi_)/12./Constants::pi;
	// compute the prefactor
	complex<InvEnergy4> pre = gRhoOmegaPi_gGammaOmegaPi_/fRho_
	Resonance::BreitWignerFW(s2,omegaMass_,omegaWidth_)/sqr(omegaMass_);
	if(imode==0) pre *=sqrt(2.);
	Complex bw(0.);
	for(unsigned int ix=imin;ix<imax;++ix) {
	Energy wid = ix==0 ? grho : rhoWidths_[ix];
	Energy2 mR2 = sqr(rhoMasses_[ix]);
	bw += mR2wgts_[ix]/(mR2-q2-Complex(0.,1.)q.mass()*wid);
	}
	pre *=bw;
	for(unsigned int ix=0;ix<3;++ix) {
	if(ix==1) continue;
	LorentzVector<complex<Energy2> > v2 = Helicity::epsilon(pout,temp[ix],momenta[2]);
	ret[ix] += prescaleHelicity::epsilon(q,v2,pout);
	}
	}
	return ret;
	}

	bool TwoPionPhotonSNDCurrent::accept(vector<int> id) {
	if(id.size()!=3){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 (npiplus==1&&ngamma==1&&npi0==1) \|\|
	(npi0==2&&ngamma==1);
	}

	unsigned int TwoPionPhotonSNDCurrent::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) && npi0==1 && ngamma==1)
	return 0;
	else
	return 1;
	}

	// output the information for the database
	void TwoPionPhotonSNDCurrent::dataBaseOutput(ofstream & output,bool header,
	bool create) const {
	if(header) output << "update decayers set parameters=\"";
	if(create) output << "create Herwig::TwoPionPhotonSNDCurrent " << 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";
	output << "newdef " << name() << ":gGammaOmegaPi " << gGammaOmegaPi_*GeV << "\n";
	output << "newdef " << name() << ":OmegaMass " << omegaMass_/GeV << "\n";
	output << "newdef " << name() << ":OmegaWidth " << omegaWidth_/GeV << "\n";
	WeakCurrent::dataBaseOutput(output,false,false);
	if(header) output << "\n\" where BINARY ThePEGName=\""
	<< fullName() << "\";" << endl;
	}
	diff --git a/Decay/WeakCurrents/TwoPionPhotonSNDCurrent.h b/Decay/WeakCurrents/TwoPionPhotonSNDCurrent.h
	--- a/Decay/WeakCurrents/TwoPionPhotonSNDCurrent.h
	+++ b/Decay/WeakCurrents/TwoPionPhotonSNDCurrent.h
	@@ -1,244 +1,243 @@
	// -- C++ --
	#ifndef Herwig_TwoPionPhotonSNDCurrent_H
	#define Herwig_TwoPionPhotonSNDCurrent_H
	//
	// This is the declaration of the TwoPionPhotonSNDCurrent class.
	//

	#include "WeakCurrent.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* The TwoPionPhotonSNDCurrent 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 TwoPionPhotonSNDCurrentInterfaces "The interfaces"
	* defined for TwoPionPhotonSNDCurrent.
	*/
	class TwoPionPhotonSNDCurrent: public WeakCurrent {

	public:

	/**
	* The default constructor.
	*/
	TwoPionPhotonSNDCurrent();

	-
	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.
	*/
	TwoPionPhotonSNDCurrent & operator=(const TwoPionPhotonSNDCurrent &);

	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_;

	/**
	* Coupling of the photon to the omega and a pion, \f$g_{\gamma\omega\pi}\f$.
	*/
	InvEnergy gGammaOmegaPi_;

	/**
	* The \f$\omega\f$ mass.
	*/
	Energy omegaMass_;

	/**
	* The \f$\omega\f$ width.
	*/
	Energy omegaWidth_;

	/**
	* The pion mass
	*/
	Energy mpi_;
	};

	}

	#endif /* Herwig_TwoPionPhotonSNDCurrent_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 2:30 PM (1 d, 11 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804744
Default Alt Text: (36 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions