Page MenuHomeHEPForge
Files F10881472

No OneTemporary
Actions

View Options

diff --git a/Decay/WeakCurrents/EtaPiPiCurrent.cc b/Decay/WeakCurrents/EtaPiPiCurrent.cc
new file mode 100644
--- /dev/null
+++ b/Decay/WeakCurrents/EtaPiPiCurrent.cc
@@ -0,0 +1,328 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the EtaPiPiCurrent class.
+//
+
+#include "EtaPiPiCurrent.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"
+
+using namespace Herwig;
+
+EtaPiPiCurrent::EtaPiPiCurrent() : fpi_(93.3*MeV) {
+ rhoMasses_ = {0.77549*GeV,1.470*GeV,1.76*GeV};
+ rhoWidths_ = {0.1494 *GeV,0.28 *GeV,0.3 *GeV};
+ amp_ = {1. , -0.385, -0.08};
+ phase_ = {0. , -5.6, -3.9};
+ // set up for the modes in the base class
+ addDecayMode(2,-1);
+ addDecayMode(1,-1);
+ addDecayMode(2,-2);
+ setInitialModes(3);
+}
+
+IBPtr EtaPiPiCurrent::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr EtaPiPiCurrent::fullclone() const {
+ return new_ptr(*this);
+}
+
+void EtaPiPiCurrent::doinit() {
+ WeakCurrent::doinit();
+ // check consistency of parametrers
+ if(rhoMasses_.size() != rhoWidths_.size())
+ throw InitException() << "Inconsistent parameters in EtaPiPiCurrent"
+ << "::doinit()" << Exception::abortnow;
+ // weights for the rho channels
+ if(amp_.size()!=phase_.size())
+ throw InitException() << "The vectors containing the weights and phase for the"
+ << " rho channel must be the same size in "
+ << "EtaPiPiCurrent::doinit()" << Exception::runerror;
+ // combine mags and phase
+ weights_.clear();
+ for(unsigned int ix=0;ix<amp_.size();++ix) {
+ weights_.push_back(amp_[ix]*(cos(phase_[ix])+Complex(0.,1.)*sin(phase_[ix])));
+ }
+}
+
+void EtaPiPiCurrent::persistentOutput(PersistentOStream & os) const {
+ os << weights_ << amp_ << phase_ << ounit(fpi_,MeV)
+ << ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV);
+}
+
+void EtaPiPiCurrent::persistentInput(PersistentIStream & is, int) {
+ is >> weights_ >> amp_ >> phase_ >> iunit(fpi_,MeV)
+ >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV);
+}
+
+// The following static variable is needed for the type
+// description system in ThePEG.
+DescribeClass<EtaPiPiCurrent,WeakCurrent>
+describeHerwigEtaPiPiCurrent("Herwig::EtaPiPiCurrent", "HwWeakCurrents.so");
+
+void EtaPiPiCurrent::Init() {
+
+ static ClassDocumentation<EtaPiPiCurrent> documentation
+ ("There is no documentation for the EtaPiPiCurrent class");
+
+ static ParVector<EtaPiPiCurrent,Energy> interfaceRhoMasses
+ ("RhoMasses",
+ "The masses of the different rho resonances for the pi pi channel",
+ &EtaPiPiCurrent::rhoMasses_, MeV, -1, 775.8*MeV, ZERO, 10000.*MeV,
+ false, false, true);
+
+ static ParVector<EtaPiPiCurrent,Energy> interfaceRhoWidths
+ ("RhoWidths",
+ "The widths of the different rho resonances for the pi pi channel",
+ &EtaPiPiCurrent::rhoWidths_, MeV, -1, 150.3*MeV, ZERO, 1000.*MeV,
+ false, false, true);
+
+ static ParVector<EtaPiPiCurrent,double> interfaceRhoMagnitude
+ ("RhoMagnitude",
+ "Magnitude of the weight of the different resonances for the pi pi channel",
+ &EtaPiPiCurrent::amp_, -1, 0., 0, 0,
+ false, false, Interface::nolimits);
+
+ static ParVector<EtaPiPiCurrent,double> interfaceRhoPhase
+ ("RhoPhase",
+ "Phase of the weight of the different resonances for the pi pi channel",
+ &EtaPiPiCurrent::phase_, -1, 0., 0, 0,
+ false, false, Interface::nolimits);
+
+ static Parameter<EtaPiPiCurrent,Energy> interfaceFPi
+ ("FPi",
+ "The pion decay constant",
+ &EtaPiPiCurrent::fpi_, MeV, 93.3*MeV, ZERO, 200.0*MeV,
+ false, false, true);
+}
+
+// complete the construction of the decay mode for integration
+bool EtaPiPiCurrent::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((imode==0 && abs(icharge)!=3) ||
+ (imode>0 && icharge !=0)) 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==0) 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;
+ }
+ }
+ // make sure that the decays are kinematically allowed
+ int iq(0),ia(0);
+ tPDVector part = particles(icharge,imode,iq,ia);
+ Energy min=ZERO;
+ for(tPDPtr p : part) min += p->massMin();
+ if(min>upp) return false;
+ // set up the resonances
+ tPDPtr res[3];
+ if(icharge==0) {
+ res[0] =getParticleData(113);
+ res[1] =getParticleData(100113);
+ res[2] =getParticleData(30113);
+ }
+ else {
+ res[0] =getParticleData(213);
+ res[1] =getParticleData(100213);
+ res[2] =getParticleData(30213);
+ if(icharge==-3) {
+ for(unsigned int ix=0;ix<3;++ix) {
+ if(res[ix]&&res[ix]->CC()) res[ix]=res[ix]->CC();
+ }
+ }
+ }
+ // create the channels
+ for(unsigned int ix=0;ix<3;++ix) {
+ if(!res[ix]) continue;
+ if(resonance && resonance != res[ix]) continue;
+ mode->addChannel((PhaseSpaceChannel(phase),ires,res[ix],ires+1,res[0],ires+1,iloc+3,
+ ires+2,iloc+1,ires+2,iloc+2));
+ }
+ // reset the masses in the intergrators
+ for(unsigned int ix=0;ix<3;++ix) {
+ if(ix<rhoMasses_.size()&&res[ix]) {
+ mode->resetIntermediate(res[ix],rhoMasses_[ix],rhoWidths_[ix]);
+ }
+ }
+ return true;
+}
+
+// the particles produced by the current
+tPDVector EtaPiPiCurrent::particles(int icharge, unsigned int imode,
+ int,int) {
+ tPDVector output(3);
+ output[2]=getParticleData(ParticleID::eta);
+ if(imode==0) {
+ output[0]=getParticleData(ParticleID::piplus);
+ output[1]=getParticleData(ParticleID::pi0);
+ if(icharge==-3) {
+ for(unsigned int ix=0;ix<output.size();++ix) {
+ if(output[ix]->CC()) output[ix]=output[ix]->CC();
+ }
+ }
+ }
+ else {
+ output[0]=getParticleData(ParticleID::piplus);
+ output[1]=getParticleData(ParticleID::piminus);
+ }
+ return output;
+}
+
+// hadronic current
+vector<LorentzPolarizationVectorE>
+EtaPiPiCurrent::current(tcPDPtr resonance,
+ IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
+ const int imode, const int ichan,Energy & scale,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ DecayIntegrator2::MEOption) const {
+ useMe();
+ // check the isospin
+ if(Itotal!=IsoSpin::IUnknown && Itotal!=IsoSpin::IOne)
+ return vector<LorentzPolarizationVectorE>();
+ int icharge = outgoing[0]->iCharge()+outgoing[1]->iCharge();
+ // check I_3
+ if(i3!=IsoSpin::I3Unknown) {
+ switch(i3) {
+ case IsoSpin::I3Zero:
+ if(imode==0) 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>();
+ }
+ }
+ Lorentz5Momentum Q=momenta[0]+momenta[1]+momenta[2];
+ Q.rescaleMass();
+ Energy2 s1 = (momenta[0]+momenta[1]).m2();
+ Complex BW1 = Resonance::BreitWignerPWave(s1,rhoMasses_[0],rhoWidths_[0],
+ momenta[0].mass(),momenta[1].mass());
+ vector<Complex> BWs = {Resonance::BreitWignerPWave(Q.mass2(),rhoMasses_[0],rhoWidths_[0],
+ momenta[0].mass(),momenta[1].mass()),
+ Resonance::BreitWignerFW(Q.mass2(),rhoMasses_[1],
+ rhoWidths_[1]*pow(Q.mass()/rhoMasses_[1],3)),
+ Resonance::BreitWignerFW(Q.mass2(),rhoMasses_[2],
+ rhoWidths_[2]*pow(Q.mass()/rhoMasses_[2],3))};
+ Complex denom = std::accumulate(weights_.begin(),weights_.end(),Complex(0.));
+ unsigned int imin=0,imax=3;
+ if(resonance) {
+ switch(resonance->id()/1000) {
+ case 0:
+ imax = 1;
+ break;
+ case 100:
+ imin = 1;
+ imax = 2;
+ break;
+ case 30 :
+ imin = 2;
+ imax = 3;
+ break;
+ default:
+ assert(false);
+ }
+ }
+ if(ichan>0) {
+ imin = ichan;
+ imax = ichan+1;
+ }
+ // form factor
+ Complex fact(0.);
+ for(unsigned int ix=imin;ix<imax;++ix) {
+ fact += weights_[ix]*BWs[ix];
+ }
+ fact *= -0.25*Complex(0.,1.)/sqr(Constants::pi)/sqrt(3.)/denom*BW1;
+ if(imode==0) fact *=sqrt(2.);
+ scale=Q.mass();
+ LorentzPolarizationVectorE output = fact/pow<3,1>(fpi_)*Q.mass()*
+ Helicity::epsilon(momenta[0],momenta[1],momenta[2]);
+ return vector<LorentzPolarizationVectorE>(1,output);
+}
+
+bool EtaPiPiCurrent::accept(vector<int> id) {
+ // check there are only three particles
+ if(id.size()!=3) return false;
+ unsigned int npip(0),npim(0),npi0(0),neta(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::eta) ++neta;
+ }
+ if( (npip==1&&npim==1&&neta==1) ||
+ (npi0==1&&npim+npip==1&&neta==1))
+ return true;
+ else
+ return false;
+}
+
+// the decay mode
+unsigned int EtaPiPiCurrent::decayMode(vector<int> idout) {
+ unsigned int npi(0);
+ for(unsigned int ix=0;ix<idout.size();++ix) {
+ if(abs(idout[ix])==ParticleID::piplus) ++npi;
+ }
+ if(npi==2) return 1;
+ else return 0;
+}
+
+// output the information for the database
+void EtaPiPiCurrent::dataBaseOutput(ofstream & output,bool header,
+ bool create) const {
+ if(header) output << "update decayers set parameters=\"";
+ if(create) output << "create Herwig::EtaPiPiCurrent "
+ << name() << " HwWeakCurrents.so\n";
+ unsigned int ix;
+ for(ix=0;ix<rhoMasses_.size();++ix) {
+ if(ix<3) output << "newdef ";
+ else output << "insert ";
+ output << name() << ":RhoMasses " << ix << " " << rhoMasses_[ix]/MeV << "\n";
+ }
+ for(ix=0;ix<rhoWidths_.size();++ix) {
+ if(ix<3) output << "newdef ";
+ else output << "insert ";
+ output << name() << ":RhoWidths " << ix << " " << rhoWidths_[ix]/MeV << "\n";
+ }
+ for(ix=0;ix<weights_.size();++ix) {
+ if(ix<3) output << "newdef ";
+ else output << "insert ";
+ output << name() << ":RhoMagnitude " << ix << " " << amp_[ix] << "\n";
+ if(ix<3) output << "newdef ";
+ else output << "insert ";
+ output << name() << ":RhoPhase " << ix << " " << phase_[ix] << "\n";
+ }
+ output << "newdef " << name() << ":FPi " << fpi_/MeV << "\n";
+ WeakCurrent::dataBaseOutput(output,false,false);
+ if(header) output << "\n\" where BINARY ThePEGName=\""
+ << fullName() << "\";" << endl;
+}
+
diff --git a/Decay/WeakCurrents/EtaPiPiCurrent.h b/Decay/WeakCurrents/EtaPiPiCurrent.h
new file mode 100644
--- /dev/null
+++ b/Decay/WeakCurrents/EtaPiPiCurrent.h
@@ -0,0 +1,213 @@
+// -*- C++ -*-
+#ifndef Herwig_EtaPiPiCurrent_H
+#define Herwig_EtaPiPiCurrent_H
+//
+// This is the declaration of the EtaPiPiCurrent class.
+//
+
+#include "WeakCurrent.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The EtaPiPiCurrent class implements the weak current for
+ * two pions and the \f$\eta\f$ meson.
+ *
+ * @see \ref EtaPiPiCurrentInterfaces "The interfaces"
+ * defined for EtaPiPiCurrent.
+ */
+class EtaPiPiCurrent: public WeakCurrent {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ EtaPiPiCurrent();
+
+
+ /** @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 two pseudoscalar
+ * mesons and the photon.
+ * @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,
+ DecayIntegrator2::MEOption meopt) const;
+
+ /**
+ * Accept the decay. Checks the particles are the allowed mode.
+ * @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.
+ * @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:
+
+ /** @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();
+ //@}
+
+private:
+
+ /**
+ * The assignment operator is private and must never be called.
+ * In fact, it should not even be implemented.
+ */
+ EtaPiPiCurrent & operator=(const EtaPiPiCurrent &);
+
+private:
+
+ /**
+ * The weights for the form factor
+ */
+ vector<Complex> weights_;
+
+ /**
+ * The amplitudes of the weights
+ */
+ vector<double> amp_;
+
+ /**
+ * The phases of the weights
+ */
+ vector<double> phase_;
+
+ /**
+ * The masses of the \f$\rho\f$ resonances.
+ */
+ vector<Energy> rhoMasses_;
+
+ /**
+ * The widths of the \f$\rho\f$ resonances.
+ */
+ vector<Energy> rhoWidths_;
+
+ /**
+ * Pion decay constant
+ */
+ Energy fpi_;
+
+};
+
+}
+
+#endif /* Herwig_EtaPiPiCurrent_H */
diff --git a/Decay/WeakCurrents/Makefile.am b/Decay/WeakCurrents/Makefile.am
--- a/Decay/WeakCurrents/Makefile.am
+++ b/Decay/WeakCurrents/Makefile.am
@@ -1,46 +1,48 @@
BUILT_SOURCES = WeakCurrents__all.cc
CLEANFILES = WeakCurrents__all.cc
WeakCurrents__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 = \
FourPionNovosibirskCurrent.h \
ScalarMesonCurrent.h\
ThreeMesonCurrentBase.h \
ThreeMesonDefaultCurrent.h\
ThreePionCLEOCurrent.h\
TwoPionRhoCurrent.h\
KPiKStarCurrent.h\
TwoPionPhotonCurrent.h\
VectorMesonCurrent.h\
FivePionCurrent.h \
KPiCurrent.h\
KaonThreeMesonCurrent.h\
TwoPionCzyzCurrent.h\
TwoKaonCzyzCurrent.h\
ThreePionCzyzCurrent.h\
-FourPionCzyzCurrent.h
+FourPionCzyzCurrent.h\
+EtaPiPiCurrent.h
DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES))
ALL_CC_FILES = \
FourPionNovosibirskCurrent.cc \
ScalarMesonCurrent.cc \
ThreeMesonCurrentBase.cc \
ThreeMesonDefaultCurrent.cc \
ThreePionCLEOCurrent.cc \
TwoPionRhoCurrent.cc \
KPiKStarCurrent.cc \
TwoPionPhotonCurrent.cc \
VectorMesonCurrent.cc \
FivePionCurrent.cc \
KPiCurrent.cc \
KaonThreeMesonCurrent.cc \
TwoPionCzyzCurrent.cc \
TwoKaonCzyzCurrent.cc \
ThreePionCzyzCurrent.cc \
-FourPionCzyzCurrent.cc
+FourPionCzyzCurrent.cc\
+EtaPiPiCurrent.cc
diff --git a/Decay/WeakCurrents/TwoPionCzyzCurrent.cc b/Decay/WeakCurrents/TwoPionCzyzCurrent.cc
--- a/Decay/WeakCurrents/TwoPionCzyzCurrent.cc
+++ b/Decay/WeakCurrents/TwoPionCzyzCurrent.cc
@@ -1,431 +1,430 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TwoPionCzyzCurrent class.
//
#include "TwoPionCzyzCurrent.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 "Herwig/Decay/ResonanceHelpers.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
TwoPionCzyzCurrent::TwoPionCzyzCurrent()
: omegaMag_(18.7e-4), omegaPhase_(0.106),
omegaMass_(782.4*MeV),omegaWidth_(8.33*MeV), beta_(2.148),
nMax_(1000) {
// various parameters
rhoMag_ = {1.,1.,0.59,0.048,0.40,0.43};
rhoPhase_ = {0.,0.,-2.20,-2.0,-2.9,1.19};
rhoMasses_ = {773.37*MeV,1490*MeV, 1870*MeV,2120*MeV,2321*MeV,2567*MeV};
rhoWidths_ = { 147.1*MeV, 429*MeV, 357*MeV, 300*MeV, 444*MeV, 491*MeV};
// set up for the modes in the base class
addDecayMode(2,-1);
addDecayMode(1,-1);
addDecayMode(2,-2);
setInitialModes(3);
}
IBPtr TwoPionCzyzCurrent::clone() const {
return new_ptr(*this);
}
IBPtr TwoPionCzyzCurrent::fullclone() const {
return new_ptr(*this);
}
void TwoPionCzyzCurrent::persistentOutput(PersistentOStream & os) const {
os << beta_ << omegaWgt_ << omegaMag_ << omegaPhase_
<< ounit(omegaMass_,GeV) << ounit(omegaWidth_,GeV)
<< rhoWgt_ << rhoMag_ << rhoPhase_
<< ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV)
<< ounit(mass_,GeV) << ounit(width_,GeV) << coup_
<< dh_ << ounit(hres_,GeV2) << ounit(h0_,GeV2) << nMax_;
}
void TwoPionCzyzCurrent::persistentInput(PersistentIStream & is, int) {
is >> beta_ >> omegaWgt_ >> omegaMag_ >> omegaPhase_
>> iunit(omegaMass_,GeV) >> iunit(omegaWidth_,GeV)
>> rhoWgt_ >> rhoMag_ >> rhoPhase_
>> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV)
>> iunit(mass_,GeV) >> iunit(width_,GeV) >> coup_
>> dh_ >> iunit(hres_,GeV2) >> iunit(h0_,GeV2) >> nMax_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<TwoPionCzyzCurrent,WeakCurrent>
describeHerwigTwoPionCzyzCurrent("Herwig::TwoPionCzyzCurrent", "HwWeakCurrents.so");
void TwoPionCzyzCurrent::Init() {
static ClassDocumentation<TwoPionCzyzCurrent> documentation
("The TwoPionCzyzCurrent class uses the currents from "
"PRD 81 094014 for the weak current with two pions",
"The current for two pions from \\cite{Czyz:2010hj} was used.",
"%\\cite{Czyz:2010hj}\n"
"\\bibitem{Czyz:2010hj}\n"
"H.~Czyz, A.~Grzelinska and J.~H.~Kuhn,\n"
"%``Narrow resonances studies with the radiative return method,''\n"
"Phys.\\ Rev.\\ D {\\bf 81} (2010) 094014\n"
"doi:10.1103/PhysRevD.81.094014\n"
"[arXiv:1002.0279 [hep-ph]].\n"
"%%CITATION = doi:10.1103/PhysRevD.81.094014;%%\n"
"%28 citations counted in INSPIRE as of 30 Jul 2018\n");
static ParVector<TwoPionCzyzCurrent,Energy> interfaceRhoMasses
("RhoMasses",
"The masses of the different rho resonances for the pi pi channel",
&TwoPionCzyzCurrent::rhoMasses_, MeV, -1, 775.8*MeV, ZERO, 10000.*MeV,
false, false, true);
static ParVector<TwoPionCzyzCurrent,Energy> interfaceRhoWidths
("RhoWidths",
"The widths of the different rho resonances for the pi pi channel",
&TwoPionCzyzCurrent::rhoWidths_, MeV, -1, 150.3*MeV, ZERO, 1000.*MeV,
false, false, true);
static ParVector<TwoPionCzyzCurrent,double> interfaceRhoMagnitude
("RhoMagnitude",
"Magnitude of the weight of the different resonances for the pi pi channel",
&TwoPionCzyzCurrent::rhoMag_, -1, 0., 0, 0,
false, false, Interface::nolimits);
static ParVector<TwoPionCzyzCurrent,double> interfaceRhoPhase
("RhoPhase",
"Phase of the weight of the different resonances for the pi pi channel",
&TwoPionCzyzCurrent::rhoPhase_, -1, 0., 0, 0,
false, false, Interface::nolimits);
static Parameter<TwoPionCzyzCurrent,unsigned int> interfacenMax
("nMax",
"The maximum number of resonances to include in the sum,"
" should be approx infinity",
&TwoPionCzyzCurrent::nMax_, 1000, 10, 10000,
false, false, Interface::limited);
static Parameter<TwoPionCzyzCurrent,double> interfacebeta
("beta",
"The beta parameter for the couplings",
&TwoPionCzyzCurrent::beta_, 2.148, 0.0, 100.,
false, false, Interface::limited);
static Parameter<TwoPionCzyzCurrent,Energy> interfaceOmegaMass
("OmegaMass",
"The mass of the omega meson",
&TwoPionCzyzCurrent::omegaMass_, MeV,782.4*MeV, 0.0*MeV, 100.0*MeV,
false, false, Interface::limited);
static Parameter<TwoPionCzyzCurrent,Energy> interfaceOmegaWidth
("OmegaWidth",
"The mass of the omega meson",
&TwoPionCzyzCurrent::omegaWidth_, MeV, 8.33*MeV, 0.0*MeV, 100.0*MeV,
false, false, Interface::limited);
static Parameter<TwoPionCzyzCurrent,double> interfaceOmegaMagnitude
("OmegaMagnitude",
"The magnitude of the omega couplings",
&TwoPionCzyzCurrent::omegaMag_, 18.7e-4, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<TwoPionCzyzCurrent,double> interfaceOmegaPhase
("OmegaPhase",
"The magnitude of the omega couplings",
&TwoPionCzyzCurrent::omegaPhase_, 0.106, 0.0, 2.*Constants::pi,
false, false, Interface::limited);
}
void TwoPionCzyzCurrent::doinit() {
WeakCurrent::doinit();
// check consistency of parametrers
if(rhoMasses_.size()!=rhoWidths_.size())
throw InitException() << "Inconsistent parameters in TwoPionCzyzCurrent"
<< "::doinit()" << Exception::abortnow;
// weights for the rho channels
if(rhoMag_.size()!=rhoPhase_.size())
throw InitException() << "The vectors containing the weights and phase for the"
<< " rho channel must be the same size in "
<< "TwoPionCzyzCurrent::doinit()" << Exception::runerror;
Complex rhoSum(0.);
for(unsigned int ix=0;ix<rhoMag_.size();++ix) {
rhoWgt_.push_back(rhoMag_[ix]*(cos(rhoPhase_[ix])+Complex(0.,1.)*sin(rhoPhase_[ix])));
if(ix>0) rhoSum +=rhoWgt_.back();
}
omegaWgt_ = omegaMag_*(cos(omegaPhase_)+Complex(0.,1.)*sin(omegaPhase_));
// set up the masses and widths of the rho resonances
double gamB(tgamma(2.-beta_));
Complex cwgt(0.);
Energy mpi(getParticleData(ParticleID::piplus)->mass());
for(unsigned int ix=0;ix<nMax_;++ix) {
// this is gam(2-beta+n)/gam(n+1)
if(ix>0) {
gamB *= ((1.-beta_+double(ix)))/double(ix);
}
Complex c_n = tgamma(beta_-0.5) /(0.5+double(ix)) / sqrt(Constants::pi) *
sin(Constants::pi*(beta_-1.-double(ix)))/Constants::pi*gamB;
if(ix%2!=0) c_n *= -1.;
// set the masses and widths
// calc for higher resonances
if(ix>=rhoMasses_.size()) {
mass_ .push_back(rhoMasses_[0]*sqrt(1.+2.*double(ix)));
width_.push_back(rhoWidths_[0]/rhoMasses_[0]*mass_.back());
}
// input for lower ones
else {
mass_ .push_back(rhoMasses_[ix]);
width_.push_back(rhoWidths_[ix]);
if(ix>0) cwgt += c_n;
}
// parameters for the gs propagators
hres_.push_back(Resonance::Hhat(sqr(mass_.back()),mass_.back(),width_.back(),mpi,mpi));
dh_ .push_back(Resonance::dHhatds(mass_.back(),width_.back(),mpi,mpi));
h0_.push_back(Resonance::H(ZERO,mass_.back(),width_.back(),mpi,mpi,dh_.back(),hres_.back()));
coup_.push_back(c_n);
}
// fix up the early weights
for(unsigned int ix=1;ix<rhoMasses_.size();++ix) {
coup_[ix] = rhoWgt_[ix]*cwgt/rhoSum;
}
}
// complete the construction of the decay mode for integration
bool TwoPionCzyzCurrent::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((imode==0 && abs(icharge)!=3) ||
(imode>0 && icharge !=0)) 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==0) 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;
}
}
// make sure that the decays are kinematically allowed
tPDPtr part[2];
if(imode==0) {
part[0]=getParticleData(ParticleID::piplus);
part[1]=getParticleData(ParticleID::pi0);
}
else {
part[0]=getParticleData(ParticleID::piplus);
part[1]=getParticleData(ParticleID::piminus);
}
Energy min(part[0]->massMin()+part[1]->massMin());
if(min>upp) return false;
// set up the resonances
tPDPtr res[3];
if(icharge==0) {
res[0] =getParticleData(113);
res[1] =getParticleData(100113);
res[2] =getParticleData(30113);
}
else {
res[0] =getParticleData(213);
res[1] =getParticleData(100213);
res[2] =getParticleData(30213);
if(icharge==-3) {
for(unsigned int ix=0;ix<3;++ix) {
if(res[ix]&&res[ix]->CC()) res[ix]=res[ix]->CC();
}
}
}
// create the channels
for(unsigned int ix=0;ix<3;++ix) {
if(!res[ix]) continue;
if(resonance && resonance != res[ix]) continue;
- PhaseSpaceChannel newChannel((PhaseSpaceChannel(phase),ires,res[ix],ires+1,iloc+1,ires+1,iloc+2));
- mode->addChannel(newChannel);
+ mode->addChannel((PhaseSpaceChannel(phase),ires,res[ix],ires+1,iloc+1,ires+1,iloc+2));
}
// reset the masses in the intergrators
for(unsigned int ix=0;ix<3;++ix) {
if(ix<rhoMasses_.size()&&res[ix]) {
mode->resetIntermediate(res[ix],rhoMasses_[ix],rhoWidths_[ix]);
}
}
return true;
}
// the particles produced by the current
tPDVector TwoPionCzyzCurrent::particles(int icharge, unsigned int imode,
int,int) {
tPDVector output(2);
if(imode==0) {
output[0]=getParticleData(ParticleID::piplus);
output[1]=getParticleData(ParticleID::pi0);
if(icharge==-3) {
for(unsigned int ix=0;ix<output.size();++ix) {
if(output[ix]->CC()) output[ix]=output[ix]->CC();
}
}
}
else {
output[0]=getParticleData(ParticleID::piplus);
output[1]=getParticleData(ParticleID::piminus);
}
return output;
}
// hadronic current
vector<LorentzPolarizationVectorE>
TwoPionCzyzCurrent::current(tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
const int imode, const int ichan,Energy & scale,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator2::MEOption) const {
useMe();
// check the isospin
if(Itotal!=IsoSpin::IUnknown && Itotal!=IsoSpin::IOne)
return vector<LorentzPolarizationVectorE>();
int icharge = outgoing[0]->iCharge()+outgoing[1]->iCharge();
// check I_3
if(i3!=IsoSpin::I3Unknown) {
switch(i3) {
case IsoSpin::I3Zero:
if(imode==0) 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>();
}
}
// momentum difference and sum of the mesons
Lorentz5Momentum pdiff(momenta[0]-momenta[1]);
Lorentz5Momentum psum (momenta[0]+momenta[1]);
psum.rescaleMass();
scale=psum.mass();
// mass2 of vector intermediate state
Energy2 q2(psum.m2());
double dot(psum*pdiff/q2);
psum *=dot;
// compute the form factor
Complex FPI=Fpi(q2,imode,ichan,resonance,momenta[0].mass(),momenta[1].mass());
// calculate the current
pdiff -= psum;
return vector<LorentzPolarizationVectorE>(1,FPI*pdiff);
}
bool TwoPionCzyzCurrent::accept(vector<int> id) {
// check there are only two particles
if(id.size()!=2) return false;
// pion modes
if((abs(id[0])==ParticleID::piplus && id[1] ==ParticleID::pi0 ) ||
( id[0] ==ParticleID::pi0 && abs(id[1])==ParticleID::piplus))
return true;
else if((id[0]==ParticleID::piminus && id[1]==ParticleID::piplus) ||
(id[0]==ParticleID::piplus && id[1]==ParticleID::piminus))
return true;
else
return false;
}
// the decay mode
unsigned int TwoPionCzyzCurrent::decayMode(vector<int> idout) {
unsigned int npi(0);
for(unsigned int ix=0;ix<idout.size();++ix) {
if(abs(idout[ix])==ParticleID::piplus) ++npi;
}
if(npi==2) return 1;
else return 0;
}
// output the information for the database
void TwoPionCzyzCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::TwoPionCzyzCurrent "
<< name() << " HwWeakCurrents.so\n";
unsigned int ix;
for(ix=0;ix<rhoMasses_.size();++ix) {
if(ix<6) output << "newdef ";
else output << "insert ";
output << name() << ":RhoMasses " << ix << " " << rhoMasses_[ix]/MeV << "\n";
}
for(ix=0;ix<rhoWidths_.size();++ix) {
if(ix<6) output << "newdef ";
else output << "insert ";
output << name() << ":RhoWidths " << ix << " " << rhoWidths_[ix]/MeV << "\n";
}
for(ix=0;ix<rhoWgt_.size();++ix) {
if(ix<6) output << "newdef ";
else output << "insert ";
output << name() << ":RhoMagnitude " << ix << " " << rhoMag_[ix] << "\n";
if(ix<6) output << "newdef ";
else output << "insert ";
output << name() << ":RhoPhase " << ix << " " << rhoPhase_[ix] << "\n";
}
output << "newdef " << name() << ":OmegaMass " << omegaMass_/MeV << "\n";
output << "newdef " << name() << ":OmegaWidth " << omegaWidth_/MeV << "\n";
output << "newdef " << name() << ":OmegaMagnitude " << omegaMag_ << "\n";
output << "newdef " << name() << ":OmegaPhase " << omegaPhase_ << "\n";
output << "newdef " << name() << ":nMax " << nMax_ << "\n";
output << "newdef " << name() << ":beta " << beta_ << "\n";
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
Complex TwoPionCzyzCurrent::Fpi(Energy2 q2,const int imode, const int ichan,
tcPDPtr resonance, Energy ma, Energy mb) const {
Complex FPI(0.);
unsigned int imin=0, imax = coup_.size();
if(ichan>0) {
imin = ichan;
imax = ichan+1;
}
if(resonance) {
switch(resonance->id()/1000) {
case 0:
imax = 1;
break;
case 100:
imin = 1;
imax = 2;
break;
case 30 :
imin = 2;
imax = 3;
break;
default:
assert(false);
}
}
for(unsigned int ix=imin;ix<imax;++ix) {
Complex term = coup_[ix]*Resonance::BreitWignerGS(q2,mass_[ix],width_[ix],
ma,mb,h0_[ix],dh_[ix],hres_[ix]);
// include rho-omega if needed
if(ix==0&&imode!=0)
term *= 1./(1.+omegaWgt_)*(1.+omegaWgt_*Resonance::BreitWignerFW(q2,omegaMass_,omegaWidth_));
FPI += term;
}
// factor for cc mode
if(imode==0) FPI *= sqrt(2.0);
return FPI;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:25 AM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4968096
Default Alt Text
(34 KB)

Event Timeline