Page MenuHomeHEPForge
Files F9514773

No OneTemporary
Actions

View Options

diff --git a/Decay/WeakCurrents/FourPionCzyzCurrent.cc b/Decay/WeakCurrents/FourPionCzyzCurrent.cc
--- a/Decay/WeakCurrents/FourPionCzyzCurrent.cc
+++ b/Decay/WeakCurrents/FourPionCzyzCurrent.cc
@@ -1,1075 +1,1111 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FourPionCzyzCurrent class.
//
#include "FourPionCzyzCurrent.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.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;
-FourPionCzyzCurrent::FourPionCzyzCurrent()
- : mpip_(140*MeV), mpi0_(140*MeV) { // /**
+FourPionCzyzCurrent::FourPionCzyzCurrent()
+ : mpip_(140*MeV), mpi0_(140*MeV), channelMap_(6,vector<int>()) {
// Masses and widths of the particles
// rho (PDG for most of current)
rhoMasses_ = {0.7755*GeV,1.459*GeV,1.72*GeV};
rhoWidths_ = {0.1494*GeV,0.4 *GeV,0.25*GeV};
// fitted for F_rho
rhoMasses_Frho_ = {0.7755*GeV,1.437*GeV,1.738*GeV,2.12*GeV};
rhoWidths_Frho_ = {0.1494*GeV,0.520*GeV,0.450*GeV,0.30*GeV};
// omega
omegaMass_ = 0.78265*GeV;
omegaWidth_ = 0.00849*GeV;
// f0
f0Mass_ = 1.35*GeV;
f0Width_ = 0.2 *GeV;
// a1
a1Mass_ = 1.23*GeV;
a1Width_ = 0.2*GeV;
// Coefficents for sum over \f$\rho\f$ resonances
beta_a1_ ={1.,-0.066,-0.021,-0.0043};
beta_f0_ ={1.,7e4,-2.5e3,1.9e3};
beta_omega_ ={1.,-0.33,0.012,-0.0053};
beta_B_ ={1.,-0.145};
beta_bar_ ={1.,0.08,-0.0075};
// couplings for the various termsz
c_a1_ = -225./GeV2;
c_f0_ = 64./GeV2;
c_omega_ = -1.47/GeV;
c_rho_ = -2.46;
// meson meson meson couplings
g_rho_pi_pi_ = 5.997;
g_omega_pi_rho_= 42.3/GeV/GeV2/GeV2;
g_rho_gamma_ = 0.1212*GeV2;
addDecayMode(2,-1);
addDecayMode(2,-1);
addDecayMode(1,-1);
addDecayMode(2,-2);
addDecayMode(1,-1);
addDecayMode(2,-2);
setInitialModes(6);
}
IBPtr FourPionCzyzCurrent::clone() const {
return new_ptr(*this);
}
IBPtr FourPionCzyzCurrent::fullclone() const {
return new_ptr(*this);
}
void FourPionCzyzCurrent::persistentOutput(PersistentOStream & os) const {
os << ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV)
<< ounit(rhoMasses_Frho_,GeV) << ounit(rhoWidths_Frho_,GeV)
<< ounit(omegaMass_,GeV) << ounit(omegaWidth_,GeV)
<< ounit(f0Mass_,GeV) << ounit(f0Width_,GeV)
<< ounit(a1Mass_,GeV) << ounit(a1Width_,GeV)
<< beta_a1_ << beta_f0_ << beta_omega_ << beta_B_ << beta_bar_
<< ounit(c_a1_,1./GeV2) << ounit(c_f0_,1./GeV2) << ounit(c_omega_,1./GeV) << c_rho_
<< g_rho_pi_pi_ << ounit(g_omega_pi_rho_,1/GeV/GeV2/GeV2) << ounit(g_rho_gamma_,GeV2)
- << ounit(mpip_,GeV) << ounit(mpi0_,GeV);
+ << ounit(mpip_,GeV) << ounit(mpi0_,GeV) << channelMap_;
}
void FourPionCzyzCurrent::persistentInput(PersistentIStream & is, int) {
is >> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV)
>> iunit(rhoMasses_Frho_,GeV) >> iunit(rhoWidths_Frho_,GeV)
>> iunit(omegaMass_,GeV) >> iunit(omegaWidth_,GeV)
>> iunit(f0Mass_,GeV) >> iunit(f0Width_,GeV)
>> iunit(a1Mass_,GeV) >> iunit(a1Width_,GeV)
>> beta_a1_ >> beta_f0_ >> beta_omega_ >> beta_B_ >> beta_bar_
>> iunit(c_a1_,1./GeV2) >> iunit(c_f0_,1./GeV2) >> iunit(c_omega_,1./GeV) >> c_rho_
>> g_rho_pi_pi_ >> iunit(g_omega_pi_rho_,1/GeV/GeV2/GeV2) >> iunit(g_rho_gamma_,GeV2)
- >> iunit(mpip_,GeV) >> iunit(mpi0_,GeV);
+ >> iunit(mpip_,GeV) >> iunit(mpi0_,GeV) >> channelMap_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<FourPionCzyzCurrent,WeakCurrent>
describeHerwigFourPionCzyzCurrent("Herwig::FourPionCzyzCurrent",
"HwWeakCurrents.so");
void FourPionCzyzCurrent::Init() {
static ClassDocumentation<FourPionCzyzCurrent> documentation
("The FourPionCzyzCurrent class is designed to implement "
"the four pion current for e+e- collisions from Phys.Rev. D77 (2008) 114005",
"The current from \\cite{Czyz:2008kw} was used for four pions.",
"\\bibitem{Czyz:2008kw}\n"
"H.~Czyz, J.~H.~Kuhn and A.~Wapienik,\n"
"%``Four-pion production in tau decays and e+e- annihilation: An Update,''\n"
"Phys.\\ Rev.\\ D {\\bf 77} (2008) 114005\n"
"doi:10.1103/PhysRevD.77.114005\n"
"[arXiv:0804.0359 [hep-ph]].\n"
"%%CITATION = doi:10.1103/PhysRevD.77.114005;%%\n"
"%35 citations counted in INSPIRE as of 02 Aug 2018\n");
static ParVector<FourPionCzyzCurrent,Energy> interfaceRhoMasses
("RhoMasses",
"The masses of the rho mesons used by default in the current",
&FourPionCzyzCurrent::rhoMasses_, GeV, -1, 0.7755*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<FourPionCzyzCurrent,Energy> interfaceRhoWidths
("RhoWidths",
"The widths of the rho mesons used by default in the current",
&FourPionCzyzCurrent::rhoWidths_, GeV, -1, 0.1494*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<FourPionCzyzCurrent,Energy> interfaceRhoMassesFrho
("RhoMassesFrho",
"The masses of the rho mesons used in the F_rho piece",
&FourPionCzyzCurrent::rhoMasses_Frho_, GeV, -1, 0.7755*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<FourPionCzyzCurrent,Energy> interfaceRhoWidthsFrho
("RhoWidthsFrho",
"The widths of the rho mesons used in the F_rho piece",
&FourPionCzyzCurrent::rhoWidths_Frho_, GeV, -1, 0.1494*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,Energy> interfaceomegaMass
("omegaMass",
"The mass of the omega meson",
&FourPionCzyzCurrent::omegaMass_, GeV, 0.78265*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,Energy> interfaceomegaWidth
("omegaWidth",
"The width of the omega meson",
&FourPionCzyzCurrent::omegaWidth_, GeV, 0.00849*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,Energy> interfaceF0Mass
("f0Mass",
"The mass of the f0 meson",
&FourPionCzyzCurrent::f0Mass_, GeV, 1.35*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,Energy> interfaceF0Width
("f0Width",
"The width of the f0 meson",
&FourPionCzyzCurrent::f0Width_, GeV, 0.2*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,Energy> interfaceA1Mass
("a1Mass",
"The mass of the a1 meson",
&FourPionCzyzCurrent::a1Mass_, GeV, 1.23*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,Energy> interfaceA1Width
("a1Width",
"The width of the a1 meson",
&FourPionCzyzCurrent::a1Width_, GeV, 0.2*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static ParVector<FourPionCzyzCurrent,double> interfacebeta_a1
("beta_a1",
"The coefficients for the sum over rho resonances in the a_1 term",
&FourPionCzyzCurrent::beta_a1_, -1, 1.0, 0, 0,
false, false, Interface::nolimits);
static ParVector<FourPionCzyzCurrent,double> interfacebeta_f0
("beta_f0",
"The coefficients for the sum over rho resonances in the f_0 term",
&FourPionCzyzCurrent::beta_f0_, -1, 1.0, 0, 0,
false, false, Interface::nolimits);
static ParVector<FourPionCzyzCurrent,double> interfacebeta_omega
("beta_omega",
"The coefficients for the sum over rho resonances in the omega term",
&FourPionCzyzCurrent::beta_omega_, -1, 1.0, 0, 0,
false, false, Interface::nolimits);
static ParVector<FourPionCzyzCurrent,double> interfacebeta_B
("beta_B",
"The coefficients for the sum over rho resonances in the B_rho term",
&FourPionCzyzCurrent::beta_B_, -1, 1.0, 0, 0,
false, false, Interface::nolimits);
static ParVector<FourPionCzyzCurrent,double> interfacebeta_bar
("beta_bar",
"The coefficients for the sum over rho resonances in the T_rho term",
&FourPionCzyzCurrent::beta_bar_, -1, 1.0, 0, 0,
false, false, Interface::nolimits);
static Parameter<FourPionCzyzCurrent,InvEnergy2> interfacec_a1
("c_a1",
"The coupling for the a_1 channel",
&FourPionCzyzCurrent::c_a1_, 1./GeV2, -255./GeV2, -1e5/GeV2, 1e5/GeV2,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,InvEnergy2> interfacec_f0
("c_f0",
"The coupling for the f_0 channel",
&FourPionCzyzCurrent::c_f0_, 1./GeV2, 64./GeV2, -1e5/GeV2, 1e5/GeV2,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,InvEnergy> interfacec_omega
("c_omega",
"The coupling for the omega channel",
&FourPionCzyzCurrent::c_omega_, 1./GeV, -1.47/GeV, -1e5/GeV, 1e5/GeV,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,double> interfacec_rho
("c_rho",
"The coupling for the rho channel",
&FourPionCzyzCurrent::c_rho_, -2.46, 0, 0,
false, false, Interface::nolimits);
static Parameter<FourPionCzyzCurrent,double> interfaceg_rho_pi_pi
("g_rho_pi_pi",
"The coupling of rho to two pions",
&FourPionCzyzCurrent::g_rho_pi_pi_, 5.997, 0, 0,
false, false, Interface::nolimits);
static Parameter<FourPionCzyzCurrent,ThePEG::Qty<std::ratio<0,1>, std::ratio<-5,1>, std::ratio<0,1> >> interfaceg_omega_pi_rho
("g_omega_pi_rho",
"The coupling of omega to rho and pi",
&FourPionCzyzCurrent::g_omega_pi_rho_, 1./GeV/GeV2/GeV2, 42.3/GeV/GeV2/GeV2, 0.0/GeV/GeV2/GeV2, 1e5/GeV/GeV2/GeV2,
false, false, Interface::limited);
static Parameter<FourPionCzyzCurrent,Energy2> interfaceg_rho_gamma
("g_rho_gamma",
"The coupling of the rho to the photon",
&FourPionCzyzCurrent::g_rho_gamma_, GeV2, 0.1212*GeV2, 0.0*GeV2, 10.0*GeV2,
false, false, Interface::limited);
}
void FourPionCzyzCurrent::doinit() {
WeakCurrent::doinit();
mpip_ = getParticleData(211)->mass();
mpi0_ = getParticleData(111)->mass();
// test of the current for a fixed momentum configuration
// Lorentz5Momentum q1(0.13061870567796208*GeV,-0.21736300316234394*GeV,
// 0.51725254282500699*GeV,0.59167288008090657*GeV);
// Lorentz5Momentum q2(-1.1388573867484255 *GeV, 0.37727761929037396 *GeV, 0.31336706796993302 *GeV, 1.2472979400305677*GeV);
// Lorentz5Momentum q3(0.11806773412672231 *GeV, 0.17873024832600765 *GeV, 0.10345508827447017 *GeV, 0.27580297667647757*GeV );
// Lorentz5Momentum q4 (7.7487017488620830E-002*GeV, 0.16118198754624435 *GeV, 6.5813182962706746E-002*GeV, 0.23620982448991124*GeV );
// q1.rescaleMass();
// q2.rescaleMass();
// q2.rescaleMass();
// q3.rescaleMass();
// cerr << q1/GeV << " " << q1.mass()/GeV << "\n";
// cerr << q2/GeV << " " << q2.mass()/GeV << "\n";
// cerr << q3/GeV << " " << q3.mass()/GeV << "\n";
// cerr << q4/GeV << " " << q4.mass()/GeV << "\n";
// Lorentz5Momentum Q(q1+q2+q3+q4);
// Q.rescaleMass();
// LorentzVector<complex<InvEnergy> > base = baseCurrent(Q.mass2(),tcPDPtr(),-1,Q,q1,q2,q3,q4);
// LorentzVector<complex<InvEnergy> > test( Complex( 376.35697290395467 , 66.446392015809550 )/GeV,
// Complex( -135.73591401998152 , 112.36912660073307 )/GeV,
// Complex( 83.215302375273723 , -54.986430577097920 )/GeV,
// Complex( -123.56434266557559 , -22.465096431505703 )/GeV);
// cerr << "current test X " << (base.x()-test.x())/(base.x()+test.x()) << "\n";
// cerr << "current test Y " << (base.y()-test.y())/(base.x()+test.y()) << "\n";
// cerr << "current test Z " << (base.z()-test.z())/(base.x()+test.z()) << "\n";
// cerr << "current test T " << (base.t()-test.t())/(base.x()+test.t()) << "\n";
}
void FourPionCzyzCurrent::createChannels(unsigned int imode,
int icharge, tcPDPtr resonance,
unsigned int iloc,int ires,
tPDVector outgoing, PhaseSpaceModePtr mode,
PhaseSpaceChannel phase,
unsigned int j1, unsigned int j2,
- unsigned int j3, unsigned int j4) {
+ unsigned int j3, unsigned int j4,
+ int & nchan) {
tPDPtr rho0[3] = {getParticleData( 113),getParticleData( 100113),getParticleData( 30113)};
tPDPtr rhop[3] = {getParticleData( 213),getParticleData( 100213),getParticleData( 30213)};
tPDPtr rhom[3] = {getParticleData(-213),getParticleData(-100213),getParticleData(-30213)};
tPDPtr omega(getParticleData(ParticleID::omega));
tPDPtr f0(getParticleData(10221));
tPDPtr a1p = getParticleData(ParticleID::a_1plus);
tPDPtr a1m = getParticleData(ParticleID::a_1minus);
tPDPtr a10 = getParticleData(ParticleID::a_10);
if(icharge==3) {
swap(rhop,rhom);
swap(a1p,a1m);
}
int rhoCharge;
tPDPtr rho;
tPDPtr rhoin;
// first the a1 channels
for(unsigned int irho=0;irho<3;++irho) {
tPDPtr rhoin = icharge==0 ? rho0[irho] : rhom[irho];
- if(resonance && rhoin!=resonance) continue;
+ if(resonance && rhoin!=resonance) {
+ nchan+=8;
+ continue;
+ }
int a1Charge;
tPDPtr a1;
for(unsigned int irho2=0;irho2<2;++irho2) {
// // find the a1
a1Charge = outgoing[j1-1]->iCharge()+outgoing[j2-1]->iCharge()+outgoing[j4-1]->iCharge();
a1 = a1Charge==0 ? a10 : a1p;
if(a1->iCharge()!=a1Charge) a1=a1m;
+ assert(abs(a1Charge)<=3);
// find the rho
rhoCharge = outgoing[j1-1]->iCharge()+outgoing[j4-1]->iCharge();
rho = rhoCharge==0 ? rho0[irho2] : rhop[irho2];
if(rho->iCharge()!=rhoCharge) rho = rhom[irho2];
+ assert(abs(rhoCharge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,a1,ires+1,iloc+j3,
ires+2,rho,ires+2,iloc+j2,ires+3,iloc+j1,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
// find the rho
if(imode!=4) {
rhoCharge = outgoing[j2-1]->iCharge()+outgoing[j4-1]->iCharge();
rho = rhoCharge==0 ? rho0[irho2] : rhop[irho2];
if(rho->iCharge()!=rhoCharge) rho = rhom[irho2];
+ assert(abs(rhoCharge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,a1,ires+1,iloc+j3,
ires+2,rho,ires+2,iloc+j1,ires+3,iloc+j2,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
}
+ else ++nchan;
// find the second a1
a1Charge = outgoing[j1-1]->iCharge()+outgoing[j2-1]->iCharge()+outgoing[j3-1]->iCharge();
a1 = a1Charge==0 ? a10 : a1p;
if(a1->iCharge()!=a1Charge) a1=a1m;
+ assert(abs(a1Charge)<=3);
// find the rho
if(imode!=4) {
rhoCharge = outgoing[j1-1]->iCharge()+outgoing[j3-1]->iCharge();
rho = rhoCharge==0 ? rho0[irho2] : rhop[irho2];
if(rho->iCharge()!=rhoCharge) rho = rhom[irho2];
+ assert(abs(rhoCharge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,a1,ires+1,iloc+j4,
ires+2,rho,ires+2,iloc+j2,ires+3,iloc+j1,ires+3,iloc+j3));
+ ++nchan; channelMap_[imode].push_back(nchan);
}
+ else
+ ++nchan;
// find the rho
if(imode!=1) {
rhoCharge = outgoing[j2-1]->iCharge()+outgoing[j3-1]->iCharge();
rho = rhoCharge==0 ? rho0[irho2] : rhop[irho2];
if(rho->iCharge()!=rhoCharge) rho = rhom[irho2];
+ assert(abs(rhoCharge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,a1,ires+1,iloc+j4,
ires+2,rho,ires+2,iloc+j1,ires+3,iloc+j2,ires+3,iloc+j3));
+ ++nchan; channelMap_[imode].push_back(nchan);
}
+ else
+ ++nchan;
}
}
// now the f_0 channel
for(unsigned int irho=0;irho<3;++irho) {
tPDPtr rhoin = icharge==0 ? rho0[irho] : rhom[irho];
if(resonance && rhoin!=resonance) continue;
rhoCharge = outgoing[j1-1]->iCharge()+outgoing[j2-1]->iCharge();
for(unsigned int irho2=0;irho2<3;++irho2) {
rho= rhoCharge==0 ? rho0[irho2] : rhop[irho2];
if(rho->iCharge()!=rhoCharge) rho = rhom[irho2];
+ assert(abs(rhoCharge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,rho,ires+1,f0,
ires+2,iloc+j1,ires+2,iloc+j2,ires+3,iloc+j3,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
}
}
// now the omega channels
if(imode>=1&&imode<=3) {
for(unsigned int irho=0;irho<3;++irho) {
tPDPtr rhoin = icharge==0 ? rho0[irho] : rhom[irho];
if(resonance && rhoin!=resonance) continue;
- rhoCharge = outgoing[j2-1]->iCharge()+outgoing[j3-1]->iCharge();
+ if(imode!=1) {
+ rhoCharge = outgoing[j2-1]->iCharge()+outgoing[j3-1]->iCharge();
+ rho= rhoCharge==0 ? rho0[0] : rhop[0];
+ if(rho->iCharge()!=rhoCharge) rho = rhom[0];
+ assert(abs(rhoCharge)<=3);
+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j1,
+ ires+2,rho,ires+2,iloc+j4,ires+3,iloc+j2,ires+3,iloc+j3));
+ ++nchan; channelMap_[imode].push_back(nchan);
+ rhoCharge = outgoing[j2-1]->iCharge()+outgoing[j4-1]->iCharge();
+ rho= rhoCharge==0 ? rho0[0] : rhop[0];
+ if(rho->iCharge()!=rhoCharge) rho = rhom[0];
+ assert(abs(rhoCharge)<=3);
+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j1,
+ ires+2,rho,ires+2,iloc+j3,ires+3,iloc+j2,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
+ rhoCharge = outgoing[j3-1]->iCharge()+outgoing[j4-1]->iCharge();
+ rho= rhoCharge==0 ? rho0[0] : rhop[0];
+ if(rho->iCharge()!=rhoCharge) rho = rhom[0];
+ assert(abs(rhoCharge)<=3);
+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j1,
+ ires+2,rho,ires+2,iloc+j2,ires+3,iloc+j3,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
+ }
+ else nchan+=3;
+ rhoCharge = outgoing[j1-1]->iCharge()+outgoing[j3-1]->iCharge();
rho= rhoCharge==0 ? rho0[0] : rhop[0];
if(rho->iCharge()!=rhoCharge) rho = rhom[0];
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j1,
- ires+2,rho,ires+2,iloc+j4,ires+3,iloc+j2,ires+3,iloc+j3));
- rhoCharge = outgoing[j2-1]->iCharge()+outgoing[j4-1]->iCharge();
+ assert(abs(rhoCharge)<=3);
+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j2,
+ ires+2,rho,ires+2,iloc+j4,ires+3,iloc+j1,ires+3,iloc+j3));
+ ++nchan; channelMap_[imode].push_back(nchan);
+ rhoCharge = outgoing[j1-1]->iCharge()+outgoing[j4-1]->iCharge();
rho= rhoCharge==0 ? rho0[0] : rhop[0];
if(rho->iCharge()!=rhoCharge) rho = rhom[0];
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j1,
- ires+2,rho,ires+2,iloc+j3,ires+3,iloc+j2,ires+3,iloc+j4));
+ assert(abs(rhoCharge)<=3);
+ mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j2,
+ ires+2,rho,ires+2,iloc+j3,ires+3,iloc+j1,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
rhoCharge = outgoing[j3-1]->iCharge()+outgoing[j4-1]->iCharge();
rho= rhoCharge==0 ? rho0[0] : rhop[0];
if(rho->iCharge()!=rhoCharge) rho = rhom[0];
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j1,
- ires+2,rho,ires+2,iloc+j2,ires+3,iloc+j3,ires+3,iloc+j4));
- rhoCharge = outgoing[j1-1]->iCharge()+outgoing[j3-1]->iCharge();
- rho= rhoCharge==0 ? rho0[0] : rhop[0];
- if(rho->iCharge()!=rhoCharge) rho = rhom[0];
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j2,
- ires+2,rho,ires+2,iloc+j4,ires+3,iloc+j1,ires+3,iloc+j3));
- rhoCharge = outgoing[j1-1]->iCharge()+outgoing[j4-1]->iCharge();
- rho= rhoCharge==0 ? rho0[0] : rhop[0];
- if(rho->iCharge()!=rhoCharge) rho = rhom[0];
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j2,
- ires+2,rho,ires+2,iloc+j3,ires+3,iloc+j1,ires+3,iloc+j4));
- rhoCharge = outgoing[j3-1]->iCharge()+outgoing[j4-1]->iCharge();
- rho= rhoCharge==0 ? rho0[0] : rhop[0];
- if(rho->iCharge()!=rhoCharge) rho = rhom[0];
+ assert(abs(rhoCharge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,omega,ires+1,iloc+j2,
ires+2,rho,ires+2,iloc+j1,ires+3,iloc+j3,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
}
}
+ else
+ nchan+=18;
// the rho channels cancel for -000 and ++--
- if(imode==0 || imode>3) return;
+ if(imode==0 || imode>3) {
+ nchan+=16;
+ return;
+ }
// rho channels
for(unsigned int irho=0;irho<2;++irho) {
tPDPtr rhoin = icharge==0 ? rho0[irho] : rhom[irho];
if(resonance && rhoin!=resonance) continue;
for(unsigned int irho1=0;irho1<2;++irho1) {
for(unsigned int irho2=0;irho2<2;++irho2) {
int rho1Charge = outgoing[j1-1]->iCharge()+outgoing[j3-1]->iCharge();
tPDPtr rho1= rho1Charge==0 ? rho0[irho1] : rhop[irho1];
if(rho1->iCharge()!=rho1Charge) rho1 = rhom[irho1];
assert(abs(rho1Charge)<=3);
int rho2Charge = outgoing[j2-1]->iCharge()+outgoing[j4-1]->iCharge();
tPDPtr rho2= rho2Charge==0 ? rho0[irho2] : rhop[irho2];
if(rho2->iCharge()!=rho2Charge) rho2 = rhom[irho2];
assert(abs(rho2Charge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,rho1,ires+1,rho2,
ires+2,iloc+j1,ires+2,iloc+j3,ires+3,iloc+j2,ires+3,iloc+j4));
+ ++nchan; channelMap_[imode].push_back(nchan);
assert(icharge==rho1Charge+rho2Charge);
if(imode!=1) {
rho1Charge = outgoing[j1-1]->iCharge()+outgoing[j4-1]->iCharge();
rho1= rho1Charge==0 ? rho0[irho1] : rhop[irho1];
if(rho1->iCharge()!=rho1Charge) rho1 = rhom[irho1];
assert(abs(rho1Charge)<=3);
rho2Charge = outgoing[j2-1]->iCharge()+outgoing[j3-1]->iCharge();
rho2= rho2Charge==0 ? rho0[irho2] : rhop[irho2];
if(rho2->iCharge()!=rho2Charge) rho2 = rhom[irho2];
assert(abs(rho2Charge)<=3);
mode->addChannel((PhaseSpaceChannel(phase),ires,rhoin,ires+1,rho1,ires+1,rho2,
ires+2,iloc+j1,ires+2,iloc+j4,ires+3,iloc+j2,ires+3,iloc+j3));
+ ++nchan; channelMap_[imode].push_back(nchan);
assert(icharge==rho1Charge+rho2Charge);
}
+ else
+ ++nchan;
}
}
}
}
// complete the construction of the decay mode for integration
bool FourPionCzyzCurrent::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<2) ||
(imode>=2&&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;
}
}
// check that the modes are kinematical allowed
Energy min(ZERO);
if(imode==0) min = mpip_+3.*mpi0_;
else if(imode==1) min = 3.*mpip_+ mpi0_;
else if(imode==2||imode==3) min = 2.*mpip_+2.*mpi0_;
else min = 4.*mpip_;
if(min>upp) return false;
// resonances we need
// get the external particles
int iq(0),ia(0);
tPDVector outgoing = particles(icharge,imode,iq,ia);
+ int nchan(-1);
+ channelMap_[imode].clear();
if(imode==0) {
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,4,1,2);
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,4,1,3);
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,3,1,4);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,4,1,2,nchan);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,4,1,3,nchan);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,3,1,4,nchan);
}
else if(imode==1) {
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,2,1,4);
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,1,2,4);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,2,1,4,nchan);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,1,2,4,nchan);
}
// pi0 pi0 pi+ pi-
else if(imode==2||imode==3) {
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,4,1,2);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,3,4,1,2,nchan);
}
else {
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,4,1,3);
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,1,4,2,3);
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,3,1,4);
- createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,1,3,2,4);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,4,1,3,nchan);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,1,4,2,3,nchan);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,2,3,1,4,nchan);
+ createChannels(imode,icharge,resonance,iloc,ires,outgoing,mode,phase,1,3,2,4,nchan);
}
return true;
}
// the particles produced by the current
tPDVector FourPionCzyzCurrent::particles(int icharge, unsigned int imode,
int,int) {
tPDVector output(4);
tPDPtr pi0=getParticleData(ParticleID::pi0);
tPDPtr pip=getParticleData(ParticleID::piplus);
tPDPtr pim=pip->CC();
if(imode==0) {
output[0]=pim;
output[1]=pi0;
output[2]=pi0;
output[3]=pi0;
}
else if(imode==1) {
output[0]=pim;
output[1]=pim;
output[2]=pip;
output[3]=pi0;
}
else if(imode==2||imode==3) {
output[0]=pip;
output[1]=pim;
output[2]=pi0;
output[3]=pi0;
}
else {
output[0]=pip;
output[1]=pip;
output[2]=pim;
output[3]=pim;
}
if(icharge==3) {
for(unsigned int ix=0;ix<output.size();++ix) {
if(output[ix]->CC()) output[ix]=output[ix]->CC();
}
}
// return the answer
return output;
}
// hadronic current
vector<LorentzPolarizationVectorE>
FourPionCzyzCurrent::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 {
int icharge(0);
for(tPDPtr out : outgoing) icharge+=out->iCharge();
// 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==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>();
}
}
useMe();
// the momenta of the particles
Lorentz5Momentum q1(momenta[0]);
Lorentz5Momentum q2(momenta[1]);
Lorentz5Momentum q3(momenta[2]);
Lorentz5Momentum q4(momenta[3]);
Lorentz5Momentum Q(q1+q2+q3+q4);
Q.rescaleMass();
scale = Q.mass();
LorentzVector<complex<InvEnergy> > output;
+ assert(ichan<int(channelMap_[imode].size()));
+ int ichannelB = ichan<0 ? -1 : channelMap_[imode][ichan];
if(imode==0) {
- if(ichan<51)
- output += baseCurrent(Q.mass2(),resonance,ichan ,Q,q3,q4,q1,q2);
- if(ichan<0 || (ichan>=51&&ichan<102))
- output += baseCurrent(Q.mass2(),resonance,ichan-51,Q,q2,q4,q1,q3);
- if(ichan<0 || ichan>=102)
- output += baseCurrent(Q.mass2(),resonance,ichan-102,Q,q2,q3,q1,q4);
+ if(ichannelB<51)
+ output += baseCurrent(Q.mass2(),resonance,ichannelB ,Q,q3,q4,q1,q2);
+ if(ichannelB<0 || (ichannelB>=67&&ichannelB<133))
+ output += baseCurrent(Q.mass2(),resonance,ichannelB-67,Q,q2,q4,q1,q3);
+ if(ichannelB<0 || ichannelB>=133)
+ output += baseCurrent(Q.mass2(),resonance,ichannelB-133,Q,q2,q3,q1,q4);
output *= sqrt(1./3.);
}
else if(imode==1) {
- if(ichan<117)
- output += baseCurrent(Q.mass2(),resonance,ichan,Q,q3,q2,q1,q4);
- if(ichan<0||ichan>=117)
- output += baseCurrent(Q.mass2(),resonance,ichan-117,Q,q3,q1,q2,q4);
+ if(ichannelB<117)
+ output += baseCurrent(Q.mass2(),resonance,ichannelB,Q,q3,q2,q1,q4);
+ if(ichannelB<0||ichannelB>=67)
+ output += baseCurrent(Q.mass2(),resonance,ichannelB-67,Q,q3,q1,q2,q4);
}
else if(imode==2||imode==3) {
- output = baseCurrent(Q.mass2(),resonance,ichan,Q,q3,q4,q1,q2);
+ output = baseCurrent(Q.mass2(),resonance,ichannelB,Q,q3,q4,q1,q2);
}
else if(imode==4||imode==5) {
- if(ichan<51)
- output += baseCurrent(Q.mass2(),resonance,ichan ,Q,q2,q4,q1,q3);
- if(ichan<0 || (ichan>=51&&ichan<102))
- output += baseCurrent(Q.mass2(),resonance,ichan-51 ,Q,q1,q4,q2,q3);
- if(ichan<0 || (ichan>=102&&ichan<153))
- output += baseCurrent(Q.mass2(),resonance,ichan-102,Q,q2,q3,q1,q4);
- if(ichan<0 || ichan>=153)
- output += baseCurrent(Q.mass2(),resonance,ichan-153,Q,q1,q3,q2,q4);
+ if(ichannelB<67)
+ output += baseCurrent(Q.mass2(),resonance,ichannelB ,Q,q2,q4,q1,q3);
+ if(ichannelB<0 || (ichannelB>=67&&ichannelB<133))
+ output += baseCurrent(Q.mass2(),resonance,ichannelB-67 ,Q,q1,q4,q2,q3);
+ if(ichannelB<0 || (ichannelB>=133&&ichannelB<200))
+ output += baseCurrent(Q.mass2(),resonance,ichannelB-133,Q,q2,q3,q1,q4);
+ if(ichannelB<0 || ichannelB>=200)
+ output += baseCurrent(Q.mass2(),resonance,ichannelB-200,Q,q1,q3,q2,q4);
}
return vector<LorentzPolarizationVectorE>(1,output*Q.mass2());
}
bool FourPionCzyzCurrent::accept(vector<int> id) {
bool allowed(false);
// check four products
if(id.size()!=4){return false;}
int npiminus=0,npiplus=0,npi0=0;
for(unsigned int ix=0;ix<id.size();++ix) {
if(id[ix]==ParticleID:: piplus) ++npiplus;
else if(id[ix]==ParticleID::piminus) ++npiminus;
else if(id[ix]==ParticleID::pi0) ++npi0;
}
if(npiminus==2&&npiplus==1&&npi0==1) allowed=true;
else if(npiminus==1&&npi0==3) allowed=true;
else if(npiplus==2&&npiminus==1&&npi0==1) allowed=true;
else if(npiplus==1&&npi0==3) allowed=true;
else if(npiplus==2&&npiminus==2) allowed=true;
else if(npiplus==1&&npiminus==1&&npi0==2) allowed=true;
return allowed;
}
// the decay mode
unsigned int FourPionCzyzCurrent::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==1) return 0;
else if(npi==2) return 2;
else if(npi==3) return 1;
else return 4;
}
// output the information for the database
void FourPionCzyzCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::FourPionCzyzCurrent "
<< name() << " HwWeakCurrents.so\n";
for(unsigned int ix=0;ix<rhoMasses_.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":RhoMasses " << ix << " " << rhoMasses_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<rhoWidths_.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":RhoWidths " << ix << " " << rhoWidths_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<rhoMasses_.size();++ix) {
if(ix<4) output << "newdef ";
else output << "insert ";
output << name() << ":RhoMassesFrho " << ix << " " << rhoMasses_Frho_[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<rhoWidths_.size();++ix) {
if(ix<4) output << "newdef ";
else output << "insert ";
output << name() << ":RhoWidthsFrho " << ix << " " << rhoWidths_Frho_[ix]/GeV << "\n";
}
output << "newdef " << name() << ":omegaMass " << omegaMass_/GeV << "\n";
output << "newdef " << name() << ":omegaWidth " << omegaWidth_/GeV << "\n";
output << "newdef " << name() << ":f0Mass " << f0Mass_/GeV << "\n";
output << "newdef " << name() << ":f0Width " << f0Width_/GeV << "\n";
output << "newdef " << name() << ":a1Mass " << a1Mass_/GeV << "\n";
output << "newdef " << name() << ":a1Width " << a1Width_/GeV << "\n";
for(unsigned int ix=0;ix<beta_a1_.size();++ix) {
if(ix<4) output << "newdef ";
else output << "insert ";
output << name() << ":beta_a1 " << ix << " " << beta_a1_[ix] << "\n";
}
for(unsigned int ix=0;ix<beta_f0_.size();++ix) {
if(ix<4) output << "newdef ";
else output << "insert ";
output << name() << ":beta_f0 " << ix << " " << beta_f0_[ix] << "\n";
}
for(unsigned int ix=0;ix<beta_omega_.size();++ix) {
if(ix<4) output << "newdef ";
else output << "insert ";
output << name() << ":beta_omega " << ix << " " << beta_omega_[ix] << "\n";
}
for(unsigned int ix=0;ix<beta_B_.size();++ix) {
if(ix<2) output << "newdef ";
else output << "insert ";
output << name() << ":beta_B " << ix << " " << beta_B_[ix] << "\n";
}
for(unsigned int ix=0;ix<beta_bar_.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":beta_bar " << ix << " " << beta_bar_[ix] << "\n";
}
output << "newdef " << name() << ":c_a1 " << c_a1_*GeV2 << "\n";
output << "newdef " << name() << ":c_f0 " << c_f0_*GeV2 << "\n";
output << "newdef " << name() << ":c_omega " << c_omega_*GeV << "\n";
output << "newdef " << name() << ":c_rho " << c_rho_ << "\n";
output << "newdef " << name() << ":g_rho_pi_pi " << g_rho_pi_pi_ << "\n";
output << "newdef " << name() << ":g_omega_pi_rho "
<< g_omega_pi_rho_*GeV2*GeV2*GeV << "\n";
output << "newdef " << name() << ":g_rho_gamma "
<<g_rho_gamma_/GeV2 << "\n";
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
LorentzVector<complex<InvEnergy> > FourPionCzyzCurrent::
baseCurrent(Energy2 Q2, tcPDPtr resonance,const int ichan,
const Lorentz5Momentum & Q , const Lorentz5Momentum & q1,
const Lorentz5Momentum & q2, const Lorentz5Momentum & q3,
const Lorentz5Momentum & q4) const {
// check the resonance
int ires0(-1);
if(resonance) {
switch(resonance->id()/1000) {
case 0:
ires0=0;
break;
case 100:
ires0=1;
break;
case 30 :
ires0=2;
break;
default:
assert(false);
}
}
// various dot products we'll need
Energy2 m12 = sqr(q1.mass()), m22 = sqr(q2.mass());
Energy2 m32 = sqr(q3.mass()), m42 = sqr(q4.mass());
Energy2 Qq1 = Q*q1, Qq2 = Q*q2, Qq3 = Q*q3, Qq4 = Q*q4;
Energy2 Qm32= Q2-2.*Qq3+m32;
Energy2 Qm42= Q2-2.*Qq4+m42;
Energy2 q1q2 = q1*q2, q1q3 = q1*q3, q1q4 = q1*q4;
Energy2 q2q3 = q2*q3, q2q4 = q2*q4, q3q4 = q3*q4;
// coefficients of the momenta
complex<InvEnergy2> c1(ZERO),c2(ZERO),c34(ZERO),cq(ZERO),c3(ZERO),c4(ZERO);
// first the a_1 terms from A.3 0804.0359 (N.B. sign due defns)
// common coefficent
if(ichan<24) {
int ires1(-1),ires2(-1),iterm(-1);
if(ichan>=0) {
ires1 = ichan/8;
ires2 = (ichan/4)%2;
iterm = ichan%4;
}
if(ires0>=0 && ires1<0) ires1=ires0;
complex<InvEnergy2> a1_fact;
if(ires1<0) {
a1_fact = -c_a1_*
Resonance::F_rho(Q2,beta_a1_,rhoMasses_Frho_,
rhoWidths_Frho_,mpip_,mpip_);
}
else {
if(ires0<0 || ires0==ires1) {
a1_fact = -c_a1_*beta_a1_[ires1]*
Resonance::BreitWignerPWave(Q2,rhoMasses_Frho_[ires1],rhoWidths_Frho_[ires1],mpip_,mpip_)
/std::accumulate(beta_a1_.begin(),beta_a1_.end(),0.);
}
else
a1_fact=ZERO;
}
// first 2 terms
if(iterm<2) {
Complex bw_a1_Qm3 = Resonance::BreitWignera1(Qm32,a1Mass_,a1Width_);
// first term
if(iterm<=0) {
Complex Brhoq1q4(0.);
if(ires2<0) {
Brhoq1q4 = bw_a1_Qm3*
Resonance::F_rho(m12+m42+2.*q1q4,beta_B_,rhoMasses_,rhoWidths_,mpip_,mpip_);
}
else {
Brhoq1q4 = bw_a1_Qm3*beta_B_[ires2]*
Resonance::BreitWignerPWave(m12+m42+2.*q1q4,rhoMasses_[ires2],rhoWidths_[ires2],mpip_,mpip_)/
std::accumulate(beta_B_.begin(),beta_B_.end(),0.);
}
double dot1 = (q1q2-q2q4)/Qm32;
double dot1B = (Qq1-Qq4)/Q2;
double dot1C = Qq3/Q2*dot1;
// coefficients of the momenta to construct the current
c1 += 0.5*a1_fact*Brhoq1q4*( 3.-dot1);
c2 += 0.5*a1_fact*Brhoq1q4*( 1.-dot1);
c34+= 0.5*a1_fact*Brhoq1q4*( 1.+dot1);
cq += 0.5*a1_fact*Brhoq1q4*(-1.+dot1-2.*dot1B-2.*dot1C);
}
// second term
if(iterm<0||iterm==1) {
Complex Brhoq2q4(0.);
if(ires2<0) {
Brhoq2q4= bw_a1_Qm3*
Resonance::F_rho(m12+m42+2.*q2q4,beta_B_,rhoMasses_,rhoWidths_,mpip_,mpip_);
}
else {
Brhoq2q4= bw_a1_Qm3*beta_B_[ires2]*
Resonance::BreitWignerPWave(m12+m42+2.*q2q4,rhoMasses_[ires2],rhoWidths_[ires2],mpip_,mpip_)/
std::accumulate(beta_B_.begin(),beta_B_.end(),0.);
}
double dot2 = (q1q2-q1q4)/Qm32;
double dot2B = (Qq2-Qq4)/Q2;
double dot2C = Qq3/Q2*dot2;
// coefficients of the momenta to construct the current
// a_1 terms
c1 += 0.5*a1_fact*Brhoq2q4*( 1.-dot2);
c2 += 0.5*a1_fact*Brhoq2q4*( 3.-dot2);
c34+= 0.5*a1_fact*Brhoq2q4*( 1.+dot2);
cq += 0.5*a1_fact*Brhoq2q4*(-1.+dot2-2.*dot2B-2.*dot2C);
}
}
// second 2 terms
if(iterm<0||iterm>=2) {
Complex bw_a1_Qm4 = Resonance::BreitWignera1(Qm42,a1Mass_,a1Width_);
// third term
if(iterm<0||iterm==2) {
Complex Brhoq1q3(0.);
if(ires2<0) {
Brhoq1q3 = bw_a1_Qm4*
Resonance::F_rho(m12+m32+2.*q1q3,beta_B_,rhoMasses_,rhoWidths_,mpip_,mpip_);
}
else {
Brhoq1q3 = bw_a1_Qm4*beta_B_[ires2]*
Resonance::BreitWignerPWave(m12+m32+2.*q1q3,rhoMasses_[ires2],rhoWidths_[ires2],mpip_,mpip_)/
std::accumulate(beta_B_.begin(),beta_B_.end(),0.);
}
double dot3 = (q1q2-q2q3)/Qm42;
double dot3B = (Qq1-Qq3)/Q2;
double dot3C = Qq4/Q2*dot3;
// coefficients of the momenta to construct the current
// a_1 terms
c1 += 0.5*a1_fact*Brhoq1q3*(-3.+dot3);
c2 += 0.5*a1_fact*Brhoq1q3*(-1.+dot3);
c34+= 0.5*a1_fact*Brhoq1q3*( 1.+dot3);
cq += 0.5*a1_fact*Brhoq1q3*( 1.-dot3+2.*dot3B+2.*dot3C);
}
// fourth term
if(iterm<0||iterm==3) {
Complex Brhoq2q3(0.);
if(ires2<0) {
Brhoq2q3 = bw_a1_Qm4*
Resonance::F_rho(m12+m32+2.*q2q3,beta_B_,rhoMasses_,rhoWidths_,mpip_,mpip_);
}
else {
Brhoq2q3 = bw_a1_Qm4*beta_B_[ires2]*
Resonance::BreitWignerPWave(m12+m32+2.*q2q3,rhoMasses_[ires2],rhoWidths_[ires2],mpip_,mpip_)/
std::accumulate(beta_B_.begin(),beta_B_.end(),0.);
}
double dot4 = (q1q2-q1q3)/Qm42;
double dot4B = (Qq2-Qq3)/Q2;
double dot4C = Qq4/Q2*dot4;
// coefficients of the momenta to construct the current
// a_1 terms
c1 += 0.5*a1_fact*Brhoq2q3*(-1.+dot4);
c2 += 0.5*a1_fact*Brhoq2q3*(-3.+dot4);
c34+= 0.5*a1_fact*Brhoq2q3*( 1.+dot4);
cq += 0.5*a1_fact*Brhoq2q3*( 1.-dot4+2.*dot4B+2.*dot4C);
}
}
}
// f_0
if(ichan<0 || (ichan>=24 && ichan<33) ) {
int ires1(-1),ires2(-2);
if(ichan>0) {
ires1 = (ichan-24)/3;
ires2 = (ichan-24)%3;
}
Complex rho1(0.);
if(ires0>=0 && ires1<0) ires1=ires0;
if(ires1<0) {
rho1 = Resonance::F_rho(Q2,beta_f0_,rhoMasses_Frho_,rhoWidths_Frho_,mpip_,mpip_);
}
else {
if(ires0<0 || ires0==ires1) {
rho1 = beta_f0_[ires1]*Resonance::BreitWignerPWave(Q2,rhoMasses_Frho_[ires1],rhoWidths_Frho_[ires1],mpip_,mpip_)/
std::accumulate(beta_f0_.begin(),beta_f0_.end(),0.);
}
else
rho1=0.;
}
Complex rho2(0.);
if(ires2<0) {
rho2 = Resonance::F_rho(m32+m42+2.*q3q4,beta_bar_,rhoMasses_,rhoWidths_,mpip_,mpip_);
}
else {
rho2 = beta_bar_[ires2]*Resonance::BreitWignerPWave(m32+m42+2.*q3q4,rhoMasses_[ires2],rhoWidths_[ires2],mpip_,mpip_)/
std::accumulate(beta_bar_.begin(),beta_bar_.end(),0.);
}
complex<InvEnergy2> f0fact = c_f0_*rho1*rho2*
Resonance::BreitWignerSWave(m12+m22+2.*q1q2,f0Mass_,f0Width_,mpip_,mpip_);
// add contribution to the coefficients
c34 -= f0fact;
cq += f0fact*(Qq3-Qq4)/Q2;
}
// omega
- if(ichan<0 || (ichan>=33&&ichan<=51) ) {
+ if(ichan<0 || (ichan>=33&&ichan<=50) ) {
int ires1(-1),iterm(-1);
if(ichan>0) {
ires1 = (ichan-33)/6;
iterm = (ichan-33)%6;
}
Complex rho1(0.);
if(ires0>=0 && ires1<0) ires1=ires0;
if(ires1<0) {
rho1 = Resonance::F_rho(Q2,beta_omega_,rhoMasses_Frho_,rhoWidths_Frho_,mpip_,mpip_);
}
else {
if(ires0<0 || ires0==ires1) {
rho1 = beta_omega_[ires1]*Resonance::BreitWignerPWave(Q2,rhoMasses_Frho_[ires1],rhoWidths_Frho_[ires1],mpip_,mpip_)/
std::accumulate(beta_omega_.begin(),beta_omega_.end(),0.);
}
else
rho1=0.;
}
complex<ThePEG::Qty<std::ratio<0,1>, std::ratio<-6,1>, std::ratio<0,1> > >
wfact = 2.*c_omega_*g_omega_pi_rho_*g_rho_pi_pi_*rho1;
if(iterm<3) {
Complex Hterm(0.);
if(iterm<0) {
Hterm = Resonance::H(rhoMasses_[0],rhoWidths_[0],m22+2.*q2q3+m32,m22+2.*q2q4+m42,
m32+2.*q3q4+m42,mpip_,mpip_);
}
else if(iterm==0)
Hterm = Resonance::BreitWignerPWave(m22+2.*q2q3+m32,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_);
else if(iterm==1)
Hterm = Resonance::BreitWignerPWave(m22+2.*q2q4+m42,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_);
else if(iterm==2)
Hterm = Resonance::BreitWignerPWave(m32+2.*q3q4+m42,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_);
else
assert(false);
complex<ThePEG::Qty<std::ratio<0,1>, std::ratio<-6,1>, std::ratio<0,1> > >
bw_omega_1 = wfact*Resonance::BreitWignerFW(m12-2.*Qq1+Q2,omegaMass_,omegaWidth_)*Hterm;
c2 -= bw_omega_1*(q1q4*Qq3-q1q3*Qq4);
c3 -= bw_omega_1*(q1q2*Qq4-q1q4*Qq2);
c4 -= bw_omega_1*(q1q3*Qq2-q1q2*Qq3);
}
if(iterm<0||iterm>=3) {
Complex Hterm(0.);
if(iterm<0) {
Hterm = Resonance::H(rhoMasses_[0],rhoWidths_[0],m12+2.*q1q3+m32,m12+2.*q1q4+m42,
m32+2.*q3q4+m42,mpip_,mpip_);
}
else if(iterm==3)
Hterm = Resonance::BreitWignerPWave(m12+2.*q1q3+m32,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_);
else if(iterm==4)
Hterm = Resonance::BreitWignerPWave(m12+2.*q1q4+m42,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_);
else if(iterm==5)
Hterm = Resonance::BreitWignerPWave(m32+2.*q3q4+m42,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_);
else
assert(false);
complex<ThePEG::Qty<std::ratio<0,1>, std::ratio<-6,1>, std::ratio<0,1> > >
bw_omega_2 = wfact*Resonance::BreitWignerFW(m22-2.*Qq2+Q2,omegaMass_,omegaWidth_)*Hterm;
c1 -= bw_omega_2*(q2q4*Qq3-q2q3*Qq4);
c3 -= bw_omega_2*(q1q2*Qq4-q2q4*Qq1);
c4 -= bw_omega_2*(q2q3*Qq1-q1q2*Qq3);
}
}
// the rho term
LorentzVector<complex<InvEnergy> > v_rho;
if(ichan<0||ichan>=51) {
int ires1(-1),ires2(-1),ires3(-1),iterm(-1);
if(ichan>0) {
- ires1 = (ichan-51)/32;
- ires2 = ((ichan-51)/16)%2;
- ires3 = ((ichan-51)/8)%2;
- iterm = (ichan-51)%8;
+ ires1 = (ichan-51)/8;
+ ires2 = ((ichan-51)/4)%2;
+ ires3 = ((ichan-51)/2)%2;
+ iterm = (ichan-51)%2;
}
// prefactor
complex<InvEnergy2> rho1;
if(ires0>=0 && ires1<0) ires1=ires0;
if(ires1<0) {
rho1 = Resonance::BreitWignerDiff(Q2,rhoMasses_[0],rhoWidths_[0],
rhoMasses_[1],rhoWidths_[1],mpip_,mpip_);
}
else {
if(ires0<0 || ires0==ires1) {
if(ires1==0)
rho1 = Resonance::BreitWignerPWave(Q2,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
else if(ires1==1)
rho1 = -Resonance::BreitWignerPWave(Q2,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
else
rho1 = ZERO;
}
else
rho1 = ZERO;
}
Complex pre_rho = c_rho_*pow(g_rho_pi_pi_,3)*g_rho_gamma_*rho1;
- complex<InvEnergy2> BW12_q1q3_A(ZERO),BW12_q1q3_B(ZERO),BW12_q1q4_A(ZERO),BW12_q1q4_B(ZERO),
- BW12_q2q3_A(ZERO),BW12_q2q3_B(ZERO),BW12_q2q4_A(ZERO),BW12_q2q4_B(ZERO);
+ complex<InvEnergy2> BW12_q1q3_A(ZERO),BW12_q1q4_A(ZERO),BW12_q2q3_B(ZERO),BW12_q2q4_B(ZERO);
if(ires2<0) {
BW12_q1q3_A = Resonance::BreitWignerDiff(m12+m32+2.*q1q3,rhoMasses_[0],rhoWidths_[0],
rhoMasses_[1],rhoWidths_[1],mpip_,mpip_);
BW12_q1q4_A = Resonance::BreitWignerDiff(m12+m42+2.*q1q4,rhoMasses_[0],rhoWidths_[0],
rhoMasses_[1],rhoWidths_[1],mpip_,mpip_);
- BW12_q2q3_A = Resonance::BreitWignerDiff(m22+m32+2.*q2q3,rhoMasses_[0],rhoWidths_[0],
+ BW12_q2q3_B = Resonance::BreitWignerDiff(m22+m32+2.*q2q3,rhoMasses_[0],rhoWidths_[0],
rhoMasses_[1],rhoWidths_[1],mpip_,mpip_);
- BW12_q2q4_A = Resonance::BreitWignerDiff(m22+m42+2.*q2q4,rhoMasses_[0],rhoWidths_[0],
+
+ BW12_q2q4_B = Resonance::BreitWignerDiff(m22+m42+2.*q2q4,rhoMasses_[0],rhoWidths_[0],
rhoMasses_[1],rhoWidths_[1],mpip_,mpip_);
+
}
else {
if(ires2==0) {
BW12_q1q3_A = Resonance::BreitWignerPWave(m12+m32+2.*q1q3,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
BW12_q1q4_A = Resonance::BreitWignerPWave(m12+m42+2.*q1q4,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
- BW12_q2q3_A = Resonance::BreitWignerPWave(m22+m32+2.*q2q3,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
- BW12_q2q4_A = Resonance::BreitWignerPWave(m22+m42+2.*q2q4,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
}
else {
BW12_q1q3_A = -Resonance::BreitWignerPWave(m12+m32+2.*q1q3,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
BW12_q1q4_A = -Resonance::BreitWignerPWave(m12+m42+2.*q1q4,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
- BW12_q2q3_A = -Resonance::BreitWignerPWave(m22+m32+2.*q2q3,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
- BW12_q2q4_A = -Resonance::BreitWignerPWave(m22+m42+2.*q2q4,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
}
}
- if(ires3<0 || ires3==ires2) {
- BW12_q1q3_B = BW12_q1q3_A;
- BW12_q1q4_B = BW12_q1q4_A;
- BW12_q2q3_B = BW12_q2q3_A;
- BW12_q2q4_B = BW12_q2q4_A;
- }
- else {
+ if(ires3>=0) {
if(ires3==0) {
- BW12_q1q3_B = Resonance::BreitWignerPWave(m12+m32+2.*q1q3,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
- BW12_q1q4_B = Resonance::BreitWignerPWave(m12+m42+2.*q1q4,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
BW12_q2q3_B = Resonance::BreitWignerPWave(m22+m32+2.*q2q3,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
BW12_q2q4_B = Resonance::BreitWignerPWave(m22+m42+2.*q2q4,rhoMasses_[0],rhoWidths_[0],mpip_,mpip_)/sqr(rhoMasses_[0]);
}
else {
- BW12_q1q3_B = -Resonance::BreitWignerPWave(m12+m32+2.*q1q3,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
- BW12_q1q4_B = -Resonance::BreitWignerPWave(m12+m42+2.*q1q4,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
BW12_q2q3_B = -Resonance::BreitWignerPWave(m22+m32+2.*q2q3,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
BW12_q2q4_B = -Resonance::BreitWignerPWave(m22+m42+2.*q2q4,rhoMasses_[1],rhoWidths_[1],mpip_,mpip_)/sqr(rhoMasses_[1]);
}
}
// now the various terms
if(iterm<=0) {
Energy2 d1 = Qq2 - Qq4 + 2.*q2q3 - 2.*q3q4;
v_rho += pre_rho*BW12_q1q3_A*(BW12_q2q4_B*d1 + 2.)*q1;
}
+ if(iterm<=0) {
+ Energy2 d5 = Qq1 - Qq3 + 2.*q1q2 - 2.*q2q3;
+ v_rho += pre_rho*BW12_q2q4_B*(BW12_q1q3_A*d5 + 2.)*q4;
+ }
if(iterm<0||iterm==1) {
- Energy2 d5 = Qq1 - Qq3 + 2.*q1q2 - 2.*q2q3;
- v_rho += pre_rho*BW12_q2q4_A*(BW12_q1q3_B*d5 + 2.)*q4;
- }
- if(iterm<0||iterm==2) {
Energy2 d2 = Qq2 - Qq3 + 2.*q2q4 - 2.*q3q4;
v_rho -= pre_rho*BW12_q1q4_A*(BW12_q2q3_B*d2 + 2.)*q1;
}
- if(iterm<0||iterm==3) {
+ if(iterm<0||iterm==1) {
Energy2 d7 = Qq1 - Qq4 + 2.*q1q2 - 2.*q2q4;
- v_rho -= pre_rho*BW12_q2q3_A*(BW12_q1q4_B*d7 + 2.)*q3;
+ v_rho -= pre_rho*BW12_q2q3_B*(BW12_q1q4_A*d7 + 2.)*q3;
}
- if(iterm<0||iterm==4) {
+ if(iterm<0||iterm==1) {
Energy2 d3 = Qq1 - Qq4 + 2.*q1q3 - 2.*q3q4;
- v_rho += pre_rho*BW12_q2q3_A*(BW12_q1q4_B*d3 + 2.)*q2;
+ v_rho += pre_rho*BW12_q2q3_B*(BW12_q1q4_A*d3 + 2.)*q2;
}
- if(iterm<0||iterm==5) {
+ if(iterm<0||iterm==1) {
Energy2 d6 = Qq2 - Qq3 + 2.*q1q2 - 2.*q1q3;
v_rho += pre_rho*BW12_q1q4_A*(BW12_q2q3_B*d6 + 2.)*q4;
}
- if(iterm<0||iterm==6) {
+ if(iterm<=0) {
Energy2 d4 = Qq1 - Qq3 + 2.*q1q4 - 2.*q3q4;
- v_rho -= pre_rho*BW12_q2q4_A*(BW12_q1q3_B*d4 + 2.)*q2;
+ v_rho -= pre_rho*BW12_q2q4_B*(BW12_q1q3_A*d4 + 2.)*q2;
}
- if(iterm<0||iterm==7) {
+ if(iterm<=0) {
Energy2 d8 = Qq2 - Qq4 + 2.*q1q2 - 2.*q1q4;
v_rho -= pre_rho*BW12_q1q3_A*(BW12_q2q4_B*d8 + 2.)*q3;
}
complex<InvEnergy2> vdot = (Q*v_rho)/Q2;
v_rho = -v_rho + vdot*Q;
}
// put everything together
return c1*q1+c2*q2+(c3+c34)*q3+(c4-c34)*q4+cq*Q+v_rho;
}
diff --git a/Decay/WeakCurrents/FourPionCzyzCurrent.h b/Decay/WeakCurrents/FourPionCzyzCurrent.h
--- a/Decay/WeakCurrents/FourPionCzyzCurrent.h
+++ b/Decay/WeakCurrents/FourPionCzyzCurrent.h
@@ -1,330 +1,335 @@
// -*- C++ -*-
#ifndef Herwig_FourPionCzyzCurrent_H
#define Herwig_FourPionCzyzCurrent_H
//
// This is the declaration of the FourPionCzyzCurrent class.
//
#include "WeakCurrent.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Decay
*
* The FourMesonCzyzCurrent class implements the currents from Phys.Rev. D77 (2008) 114005
* for 4 pions
* @see WeakCurrent.
* @see \ref FourPionCzyzCurrentInterfaces "The interfaces"
* defined for FourPionCzyzCurrent.
*
*/
class FourPionCzyzCurrent: public WeakCurrent {
public:
/**
* The default constructor.
*/
FourPionCzyzCurrent();
/** @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:
/**
* Create the channels for 1 term in the current
*/
void createChannels(unsigned int imode,
int icharge, tcPDPtr resonance,
unsigned int iloc,int ires,
tPDVector outgoing, PhaseSpaceModePtr mode,
PhaseSpaceChannel channel,
unsigned int j1, unsigned int j2,
- unsigned int j3, unsigned int j4);
+ unsigned int j3, unsigned int j4,
+ int & nchan);
/**
* Basis current in terms of which all the others can be calculated
*/
LorentzVector<complex<InvEnergy> > baseCurrent(Energy2 Q2,
tcPDPtr resonance,
const int ichan,
const Lorentz5Momentum & Q,
const Lorentz5Momentum & q1,
const Lorentz5Momentum & q2,
const Lorentz5Momentum & q3,
const Lorentz5Momentum & q4) const;
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 and
* 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.
*/
FourPionCzyzCurrent & operator=(const FourPionCzyzCurrent &);
private:
/**
* Masses and widths of the particles
*/
//@{
/**
* Rho masses (PDG for most of current)
*/
vector<Energy> rhoMasses_;
/**
* Rho widths (PDG for most of current)
*/
vector<Energy> rhoWidths_;
/**
* Rho masses for the \f$F_\rho\f$ piece
*/
vector<Energy> rhoMasses_Frho_;
/**
* Rho widths for the \f$F_\rho\f$ piece
*/
vector<Energy> rhoWidths_Frho_;
/**
* Omega mass
*/
Energy omegaMass_;
/**
* Omega widths
*/
Energy omegaWidth_;
/**
* \f$f_0\f$ mass
*/
Energy f0Mass_;
/**
* \f$f_0\f$ width
*/
Energy f0Width_;
/**
* \f$a_1\f$ mass
*/
Energy a1Mass_;
/**
* \f$a_1\f$ width
*/
Energy a1Width_;
//@}
/**
* Couplings in the model
*/
//@{
/**
* Coefficents for sum over \f$\rho\f$ resonances in \f$a_1\f$ term
*/
vector<double> beta_a1_;
/**
* Coefficents for sum over \f$\rho\f$ resonances in \f$f_0\f$ term
*/
vector<double> beta_f0_;
/**
* Coefficents for sum over \f$\rho\f$ resonances in \f$\omega\f$ term
*/
vector<double> beta_omega_;
/**
* Coefficents for sum over \f$\rho\f$ resonances in \f$B_\rho\f$ term
*/
vector<double> beta_B_;
/**
* Coefficents for sum over \f$\rho\f$ resonances in \f$T_\rho\f$ term
*/
vector<double> beta_bar_;
/**
* Coupling for the \f$a_1\f$ term
*/
InvEnergy2 c_a1_;
/**
* Coupling for the \f$f_0\f$ term
*/
InvEnergy2 c_f0_;
/**
* Coupling for the \f$\omega\f$ term
*/
InvEnergy c_omega_;
/**
* Coupling for the \f\rho\f$ term
*/
double c_rho_;
/**
* \f$g_{\rho\pi\pi}\f$
*/
double g_rho_pi_pi_;
/**
* \f$g_{\omega\pi\prho}\f$
*/
ThePEG::Qty<std::ratio<0,1>, std::ratio<-5,1>, std::ratio<0,1> > g_omega_pi_rho_;
/**
* \f$g_{\rho\gamma}\f$
*/
Energy2 g_rho_gamma_;
//@}
/**
* Pion masses
*/
Energy mpip_, mpi0_;
+ /**
+ * Map for the phase-space channels
+ */
+ vector<vector<int> > channelMap_;
};
}
#endif /* Herwig_FourPionCzyzCurrent_H */

File Metadata

Mime Type
text/x-diff
Expires
Mon, Feb 24, 6:43 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4494603
Default Alt Text
(59 KB)

Event Timeline