diff --git a/Decay/WeakCurrents/FourPionCzyzCurrent.cc b/Decay/WeakCurrents/FourPionCzyzCurrent.cc
--- a/Decay/WeakCurrents/FourPionCzyzCurrent.cc
+++ b/Decay/WeakCurrents/FourPionCzyzCurrent.cc
@@ -1,1120 +1,1120 @@
// -*- 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), 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};
+ rhoMasses_Frho_ = {0.7755*GeV,1.437*GeV ,1.738*GeV ,2.12*GeV};
+ rhoWidths_Frho_ = {0.1494*GeV,0.6784258847826656*GeV,0.8049153117715863*GeV,0.20924569673342333*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_a1_ ={1.,-0.051864694215520334, -0.0416090742847935,-0.0018940213981381284 };
+ beta_f0_ ={1., 73673.55747104406, -26116.259838644528, 332.84652898870786};
+ beta_omega_ ={1., -0.3668543468911129,0.03624673983924276 , -0.0047186283455050064 };
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;
+ beta_bar_ ={1.,0.08,-0.0075};
+ // couplings for the various terms
+ c_a1_ = -201.7929015686851/GeV2;
+ c_f0_ = 124.09577796839454/GeV2;
+ c_omega_ = -1.5792052826500913/GeV;
+ c_rho_ = -2.308949096570198;
// 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;
+ 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) << 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) >> 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,
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) {
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;
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];
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];
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];
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) {
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,
FlavourInfo flavour,
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(flavour.I!=IsoSpin::IUnknown) {
if(flavour.I!=IsoSpin::IOne) return false;
}
// check I_3
if(flavour.I3!=IsoSpin::I3Unknown) {
switch(flavour.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;
break;
default:
return false;
}
}
// and other flavour
if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return false;
if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return false;
if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) 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,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,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,nchan);
}
else {
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,
FlavourInfo flavour,
const int imode, const int ichan,Energy & scale,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption) const {
int icharge(0);
for(tPDPtr out : outgoing) icharge+=out->iCharge();
// check the total isospin
if(flavour.I!=IsoSpin::IUnknown) {
if(flavour.I!=IsoSpin::IOne) return vector<LorentzPolarizationVectorE>();
}
// check I_3
if(flavour.I3!=IsoSpin::I3Unknown) {
switch(flavour.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>();
break;
default:
return vector<LorentzPolarizationVectorE>();
}
}
if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return vector<LorentzPolarizationVectorE>();
if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return vector<LorentzPolarizationVectorE>();
if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) 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(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(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,ichannelB,Q,q3,q4,q1,q2);
}
else if(imode==4||imode==5) {
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<=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)/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_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_B = Resonance::BreitWignerDiff(m22+m32+2.*q2q3,rhoMasses_[0],rhoWidths_[0],
rhoMasses_[1],rhoWidths_[1],mpip_,mpip_);
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]);
}
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]);
}
}
if(ires3>=0) {
if(ires3==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_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 d2 = Qq2 - Qq3 + 2.*q2q4 - 2.*q3q4;
v_rho -= pre_rho*BW12_q1q4_A*(BW12_q2q3_B*d2 + 2.)*q1;
}
if(iterm<0||iterm==1) {
Energy2 d7 = Qq1 - Qq4 + 2.*q1q2 - 2.*q2q4;
v_rho -= pre_rho*BW12_q2q3_B*(BW12_q1q4_A*d7 + 2.)*q3;
}
if(iterm<0||iterm==1) {
Energy2 d3 = Qq1 - Qq4 + 2.*q1q3 - 2.*q3q4;
v_rho += pre_rho*BW12_q2q3_B*(BW12_q1q4_A*d3 + 2.)*q2;
}
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) {
Energy2 d4 = Qq1 - Qq3 + 2.*q1q4 - 2.*q3q4;
v_rho -= pre_rho*BW12_q2q4_B*(BW12_q1q3_A*d4 + 2.)*q2;
}
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;
}