Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8723801
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
141 KB
Subscribers
None
View Options
diff --git a/Decay/WeakCurrents/EtaPiPiCurrent.cc b/Decay/WeakCurrents/EtaPiPiCurrent.cc
--- a/Decay/WeakCurrents/EtaPiPiCurrent.cc
+++ b/Decay/WeakCurrents/EtaPiPiCurrent.cc
@@ -1,328 +1,327 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the EtaPiPiCurrent class.
//
#include "EtaPiPiCurrent.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Helicity/epsilon.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
EtaPiPiCurrent::EtaPiPiCurrent() : fpi_(93.3*MeV) {
rhoMasses_ = {0.77549*GeV,1.470*GeV,1.76*GeV};
rhoWidths_ = {0.1494 *GeV,0.28 *GeV,0.3 *GeV};
amp_ = {1. , -0.385, -0.08};
phase_ = {0. , -5.6, -3.9};
// set up for the modes in the base class
addDecayMode(2,-1);
addDecayMode(1,-1);
addDecayMode(2,-2);
setInitialModes(3);
}
IBPtr EtaPiPiCurrent::clone() const {
return new_ptr(*this);
}
IBPtr EtaPiPiCurrent::fullclone() const {
return new_ptr(*this);
}
void EtaPiPiCurrent::doinit() {
WeakCurrent::doinit();
// check consistency of parametrers
if(rhoMasses_.size() != rhoWidths_.size())
throw InitException() << "Inconsistent parameters in EtaPiPiCurrent"
<< "::doinit()" << Exception::abortnow;
// weights for the rho channels
if(amp_.size()!=phase_.size())
throw InitException() << "The vectors containing the weights and phase for the"
<< " rho channel must be the same size in "
<< "EtaPiPiCurrent::doinit()" << Exception::runerror;
// combine mags and phase
weights_.clear();
for(unsigned int ix=0;ix<amp_.size();++ix) {
weights_.push_back(amp_[ix]*(cos(phase_[ix])+Complex(0.,1.)*sin(phase_[ix])));
}
}
void EtaPiPiCurrent::persistentOutput(PersistentOStream & os) const {
os << weights_ << amp_ << phase_ << ounit(fpi_,MeV)
<< ounit(rhoMasses_,GeV) << ounit(rhoWidths_,GeV);
}
void EtaPiPiCurrent::persistentInput(PersistentIStream & is, int) {
is >> weights_ >> amp_ >> phase_ >> iunit(fpi_,MeV)
>> iunit(rhoMasses_,GeV) >> iunit(rhoWidths_,GeV);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<EtaPiPiCurrent,WeakCurrent>
describeHerwigEtaPiPiCurrent("Herwig::EtaPiPiCurrent", "HwWeakCurrents.so");
void EtaPiPiCurrent::Init() {
static ClassDocumentation<EtaPiPiCurrent> documentation
("There is no documentation for the EtaPiPiCurrent class");
static ParVector<EtaPiPiCurrent,Energy> interfaceRhoMasses
("RhoMasses",
"The masses of the different rho resonances for the pi pi channel",
&EtaPiPiCurrent::rhoMasses_, MeV, -1, 775.8*MeV, ZERO, 10000.*MeV,
false, false, true);
static ParVector<EtaPiPiCurrent,Energy> interfaceRhoWidths
("RhoWidths",
"The widths of the different rho resonances for the pi pi channel",
&EtaPiPiCurrent::rhoWidths_, MeV, -1, 150.3*MeV, ZERO, 1000.*MeV,
false, false, true);
static ParVector<EtaPiPiCurrent,double> interfaceRhoMagnitude
("RhoMagnitude",
"Magnitude of the weight of the different resonances for the pi pi channel",
&EtaPiPiCurrent::amp_, -1, 0., 0, 0,
false, false, Interface::nolimits);
static ParVector<EtaPiPiCurrent,double> interfaceRhoPhase
("RhoPhase",
"Phase of the weight of the different resonances for the pi pi channel",
&EtaPiPiCurrent::phase_, -1, 0., 0, 0,
false, false, Interface::nolimits);
static Parameter<EtaPiPiCurrent,Energy> interfaceFPi
("FPi",
"The pion decay constant",
&EtaPiPiCurrent::fpi_, MeV, 93.3*MeV, ZERO, 200.0*MeV,
false, false, true);
}
// complete the construction of the decay mode for integration
bool EtaPiPiCurrent::createMode(int icharge, tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
unsigned int imode,PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp ) {
// check the charge
if((imode==0 && abs(icharge)!=3) ||
(imode>0 && icharge !=0)) return false;
// check the total isospin
if(Itotal!=IsoSpin::IUnknown) {
if(Itotal!=IsoSpin::IOne) return false;
}
// check I_3
if(i3!=IsoSpin::I3Unknown) {
switch(i3) {
case IsoSpin::I3Zero:
if(imode==0) return false;
break;
case IsoSpin::I3One:
if(imode==1 || icharge ==-3) return false;
break;
case IsoSpin::I3MinusOne:
if(imode==1 || icharge ==3) return false;
default:
return false;
}
}
// make sure that the decays are kinematically allowed
int iq(0),ia(0);
tPDVector part = particles(icharge,imode,iq,ia);
Energy min=ZERO;
for(tPDPtr p : part) min += p->massMin();
if(min>upp) return false;
// set up the resonances
tPDPtr res[3];
if(icharge==0) {
res[0] =getParticleData(113);
res[1] =getParticleData(100113);
res[2] =getParticleData(30113);
}
else {
res[0] =getParticleData(213);
res[1] =getParticleData(100213);
res[2] =getParticleData(30213);
if(icharge==-3) {
for(unsigned int ix=0;ix<3;++ix) {
if(res[ix]&&res[ix]->CC()) res[ix]=res[ix]->CC();
}
}
}
// create the channels
for(unsigned int ix=0;ix<3;++ix) {
if(!res[ix]) continue;
if(resonance && resonance != res[ix]) continue;
mode->addChannel((PhaseSpaceChannel(phase),ires,res[ix],ires+1,res[0],ires+1,iloc+3,
ires+2,iloc+1,ires+2,iloc+2));
}
// reset the masses in the intergrators
for(unsigned int ix=0;ix<3;++ix) {
if(ix<rhoMasses_.size()&&res[ix]) {
mode->resetIntermediate(res[ix],rhoMasses_[ix],rhoWidths_[ix]);
}
}
return true;
}
// the particles produced by the current
tPDVector EtaPiPiCurrent::particles(int icharge, unsigned int imode,
int,int) {
tPDVector output(3);
+ output[0]=getParticleData(ParticleID::piplus);
output[2]=getParticleData(ParticleID::eta);
if(imode==0) {
- output[0]=getParticleData(ParticleID::piplus);
output[1]=getParticleData(ParticleID::pi0);
if(icharge==-3) {
for(unsigned int ix=0;ix<output.size();++ix) {
if(output[ix]->CC()) output[ix]=output[ix]->CC();
}
}
}
else {
- output[0]=getParticleData(ParticleID::piplus);
output[1]=getParticleData(ParticleID::piminus);
}
return output;
}
// hadronic current
vector<LorentzPolarizationVectorE>
EtaPiPiCurrent::current(tcPDPtr resonance,
- IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
- const int imode, const int ichan,Energy & scale,
- const tPDVector & outgoing,
- const vector<Lorentz5Momentum> & momenta,
- DecayIntegrator::MEOption) const {
+ IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
+ const int imode, const int ichan,Energy & scale,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ DecayIntegrator::MEOption) const {
useMe();
// check the isospin
if(Itotal!=IsoSpin::IUnknown && Itotal!=IsoSpin::IOne)
return vector<LorentzPolarizationVectorE>();
int icharge = outgoing[0]->iCharge()+outgoing[1]->iCharge();
// check I_3
if(i3!=IsoSpin::I3Unknown) {
switch(i3) {
case IsoSpin::I3Zero:
if(imode==0) return vector<LorentzPolarizationVectorE>();
break;
case IsoSpin::I3One:
if(imode==1 || icharge ==-3) return vector<LorentzPolarizationVectorE>();
break;
case IsoSpin::I3MinusOne:
if(imode==1 || icharge ==3) return vector<LorentzPolarizationVectorE>();
default:
return vector<LorentzPolarizationVectorE>();
}
}
Lorentz5Momentum Q=momenta[0]+momenta[1]+momenta[2];
Q.rescaleMass();
Energy2 s1 = (momenta[0]+momenta[1]).m2();
Complex BW1 = Resonance::BreitWignerPWave(s1,rhoMasses_[0],rhoWidths_[0],
momenta[0].mass(),momenta[1].mass());
vector<Complex> BWs = {Resonance::BreitWignerPWave(Q.mass2(),rhoMasses_[0],rhoWidths_[0],
momenta[0].mass(),momenta[1].mass()),
Resonance::BreitWignerFW(Q.mass2(),rhoMasses_[1],
rhoWidths_[1]*pow(Q.mass()/rhoMasses_[1],3)),
Resonance::BreitWignerFW(Q.mass2(),rhoMasses_[2],
rhoWidths_[2]*pow(Q.mass()/rhoMasses_[2],3))};
Complex denom = std::accumulate(weights_.begin(),weights_.end(),Complex(0.));
unsigned int imin=0,imax=3;
if(resonance) {
switch(resonance->id()/1000) {
case 0:
imax = 1;
break;
case 100:
imin = 1;
imax = 2;
break;
case 30 :
imin = 2;
imax = 3;
break;
default:
assert(false);
}
}
if(ichan>0) {
imin = ichan;
imax = ichan+1;
}
// form factor
Complex fact(0.);
for(unsigned int ix=imin;ix<imax;++ix) {
fact += weights_[ix]*BWs[ix];
}
fact *= -0.25*Complex(0.,1.)/sqr(Constants::pi)/sqrt(3.)/denom*BW1;
if(imode==0) fact *=sqrt(2.);
scale=Q.mass();
LorentzPolarizationVectorE output = fact/pow<3,1>(fpi_)*Q.mass()*
Helicity::epsilon(momenta[0],momenta[1],momenta[2]);
return vector<LorentzPolarizationVectorE>(1,output);
}
bool EtaPiPiCurrent::accept(vector<int> id) {
// check there are only three particles
if(id.size()!=3) return false;
unsigned int npip(0),npim(0),npi0(0),neta(0);
for(unsigned int ix=0;ix<id.size();++ix) {
if(id[ix]==ParticleID::piplus) ++npip;
else if(id[ix]==ParticleID::piminus) ++npim;
else if(id[ix]==ParticleID::pi0) ++npi0;
else if(id[ix]==ParticleID::eta) ++neta;
}
if( (npip==1&&npim==1&&neta==1) ||
(npi0==1&&npim+npip==1&&neta==1))
return true;
else
return false;
}
// the decay mode
unsigned int EtaPiPiCurrent::decayMode(vector<int> idout) {
unsigned int npi(0);
for(unsigned int ix=0;ix<idout.size();++ix) {
if(abs(idout[ix])==ParticleID::piplus) ++npi;
}
if(npi==2) return 1;
else return 0;
}
// output the information for the database
void EtaPiPiCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::EtaPiPiCurrent "
<< name() << " HwWeakCurrents.so\n";
unsigned int ix;
for(ix=0;ix<rhoMasses_.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":RhoMasses " << ix << " " << rhoMasses_[ix]/MeV << "\n";
}
for(ix=0;ix<rhoWidths_.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":RhoWidths " << ix << " " << rhoWidths_[ix]/MeV << "\n";
}
for(ix=0;ix<weights_.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":RhoMagnitude " << ix << " " << amp_[ix] << "\n";
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":RhoPhase " << ix << " " << phase_[ix] << "\n";
}
output << "newdef " << name() << ":FPi " << fpi_/MeV << "\n";
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/WeakCurrents/EtaPiPiDefaultCurrent.cc b/Decay/WeakCurrents/EtaPiPiDefaultCurrent.cc
--- a/Decay/WeakCurrents/EtaPiPiDefaultCurrent.cc
+++ b/Decay/WeakCurrents/EtaPiPiDefaultCurrent.cc
@@ -1,1252 +1,364 @@
// -*- C++ -*-
//
// EtaPiPiDefaultCurrent.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the EtaPiPiDefaultCurrent class.
//
#include "EtaPiPiDefaultCurrent.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
using namespace ThePEG;
DescribeClass<EtaPiPiDefaultCurrent,WeakCurrent>
describeHerwigEtaPiPiDefaultCurrent("Herwig::EtaPiPiDefaultCurrent",
"HwWeakCurrents.so");
-HERWIG_INTERPOLATOR_CLASSDESC(EtaPiPiDefaultCurrent,Energy,Energy2)
-
EtaPiPiDefaultCurrent::EtaPiPiDefaultCurrent() {
- // the quarks for the different modes
+ // set up for the modes in the base class
addDecayMode(2,-1);
- addDecayMode(2,-1);
- addDecayMode(2,-1);
- addDecayMode(2,-1);
- addDecayMode(2,-1);
- addDecayMode(2,-3);
- addDecayMode(2,-3);
- addDecayMode(2,-3);
- addDecayMode(2,-1);
- addDecayMode(2,-1);
- addDecayMode(2,-1);
- addDecayMode(2,-1);
- setInitialModes(12);
+ addDecayMode(1,-1);
+ addDecayMode(2,-2);
+ setInitialModes(3);
// the pion decay constant
_fpi=130.7*MeV/sqrt(2.);
- _mpi=ZERO;_mK=ZERO;
+ _mpi=ZERO;
// set the initial weights for the resonances
// the rho weights
- _rhoF123wgts.push_back(1.0);_rhoF123wgts.push_back(-0.145);
- _rhoF123wgts.push_back(0.);
- _rhoF5wgts.push_back(-26.);_rhoF5wgts.push_back(6.5);
- _rhoF5wgts.push_back(1.);
- // the Kstar weights
- _kstarF123wgts.push_back(1.);
- _kstarF5wgts.push_back(1.);
- // relative rho/Kstar weights
- _rhoKstarwgt=-0.2;
- // local values of the a_1 parameters
- _a1parameters=true;_a1mass=1.251*GeV;_a1width=0.599*GeV;
- _a1opt=true;
- // local values of the K_1 parameters
- _k1parameters=true;_k1mass=1.402*GeV,_k1width=0.174*GeV;
+ _rhoF123wgts = { 1.0,-0.145,0.};
+ _rhoF5wgts = {-26., 6.5,1.};
// local values of the rho parameters
- _rhoparameters=true;
- _rhoF123masses.push_back(0.773*GeV);_rhoF123masses.push_back(1.370*GeV);
- _rhoF123masses.push_back(1.750*GeV);
- _rhoF123widths.push_back(0.145*GeV);_rhoF123widths.push_back(0.510*GeV);
- _rhoF123widths.push_back(0.120*GeV);
- _rhoF5masses.push_back(0.773*GeV);_rhoF5masses.push_back(1.500*GeV);
- _rhoF5masses.push_back(1.750*GeV);
- _rhoF5widths.push_back(0.145*GeV);_rhoF5widths.push_back(0.220*GeV);
- _rhoF5widths.push_back(0.120*GeV);
- // local values for the Kstar parameters
- _kstarparameters=true;
- _kstarF123masses.push_back(0.8921*GeV);
- _kstarF123widths.push_back(0.0513*GeV);
- _kstarF5masses.push_back(0.8921*GeV);
- _kstarF5widths.push_back(0.0513*GeV);
- // initialization of the a_1 running width
- _initializea1=false;
- double a1q2in[200]={0,15788.6,31577.3,47365.9,63154.6,78943.2,94731.9,110521,
- 126309,142098,157886,173675,189464,205252,221041,236830,
- 252618,268407,284196,299984,315773,331562,347350,363139,
- 378927,394716,410505,426293,442082,457871,473659,489448,
- 505237,521025,536814,552603,568391,584180,599969,615757,
- 631546,647334,663123,678912,694700,710489,726278,742066,
- 757855,773644,789432,805221,821010,836798,852587,868375,
- 884164,899953,915741,931530,947319,963107,978896,994685,
- 1.01047e+06,1.02626e+06,1.04205e+06,1.05784e+06,1.07363e+06,
- 1.08942e+06,1.10521e+06,1.12099e+06,1.13678e+06,1.15257e+06,
- 1.16836e+06,1.18415e+06,1.19994e+06,1.21573e+06,1.23151e+06,
- 1.2473e+06,1.26309e+06,1.27888e+06,1.29467e+06,1.31046e+06,
- 1.32625e+06,1.34203e+06,1.35782e+06,1.37361e+06,1.3894e+06,
- 1.40519e+06,1.42098e+06,1.43677e+06,1.45256e+06,1.46834e+06
- ,1.48413e+06,1.49992e+06,1.51571e+06,1.5315e+06,1.54729e+06,
- 1.56308e+06,1.57886e+06,1.59465e+06,1.61044e+06,1.62623e+06,
- 1.64202e+06,1.65781e+06,1.6736e+06,1.68939e+06,1.70517e+06,
- 1.72096e+06,1.73675e+06,1.75254e+06,1.76833e+06,1.78412e+06,
- 1.79991e+06,1.81569e+06,1.83148e+06,1.84727e+06,1.86306e+06,
- 1.87885e+06,1.89464e+06,1.91043e+06,1.92621e+06,1.942e+06,
- 1.95779e+06,1.97358e+06,1.98937e+06,2.00516e+06,2.02095e+06,
- 2.03674e+06,2.05252e+06,2.06831e+06,2.0841e+06,2.09989e+06,
- 2.11568e+06,2.13147e+06,2.14726e+06,2.16304e+06,2.17883e+06,
- 2.19462e+06,2.21041e+06,2.2262e+06,2.24199e+06,2.25778e+06,
- 2.27356e+06,2.28935e+06,2.30514e+06,2.32093e+06,2.33672e+06,
- 2.35251e+06,2.3683e+06,2.38409e+06,2.39987e+06,2.41566e+06,
- 2.43145e+06,2.44724e+06,2.46303e+06,2.47882e+06,2.49461e+06,
- 2.51039e+06,2.52618e+06,2.54197e+06,2.55776e+06,2.57355e+06,
- 2.58934e+06,2.60513e+06,2.62092e+06,2.6367e+06,2.65249e+06,
- 2.66828e+06,2.68407e+06,2.69986e+06,2.71565e+06,2.73144e+06,
- 2.74722e+06,2.76301e+06,2.7788e+06,2.79459e+06,2.81038e+06,
- 2.82617e+06,2.84196e+06,2.85774e+06,2.87353e+06,2.88932e+06,
- 2.90511e+06,2.9209e+06,2.93669e+06,2.95248e+06,2.96827e+06,
- 2.98405e+06,2.99984e+06,3.01563e+06,3.03142e+06,3.04721e+06,
- 3.063e+06,3.07879e+06,3.09457e+06,3.11036e+06,3.12615e+06,
- 3.14194e+06};
- double a1widthin[200]={0,0,0,0,0,0,0,0,0,0,0,0,0.00153933,0.0136382,0.0457614,
- 0.105567,0.199612,0.333825,0.513831,0.745192,1.0336,1.38501,
- 1.80581,2.30295,2.88403,3.5575,4.33278,5.22045,6.23243,
- 7.38223,8.68521,10.1589,11.8234,13.7018,15.8206,18.2107,
- 20.9078,23.9533,27.3954,31.2905,35.7038,40.7106,46.3984,
- 52.8654,60.2207,68.581,78.0637,88.7754,100.794,114.145,
- 128.783,144.574,161.299,178.683,196.426,214.248,231.908,
- 249.221,266.059,282.336,298.006,313.048,327.46,341.254,
- 354.448,367.066,379.133,390.677,401.726,412.304,422.439,
- 432.155,441.474,450.419,459.01,467.267,475.207,482.847,
- 490.203,497.29,504.121,510.71,517.068,523.207,529.138,
- 534.869,540.411,545.776,550.961,556.663,560.851,565.566,
- 570.137,574.569,578.869,583.041,587.091,591.023,594.843,
- 598.553,602.16,605.664,609.072,612.396,615.626,618.754,
- 621.796,624.766,627.656,630.47,633.21,635.878,638.5,
- 641.006,643.471,645.873,648.213,650.493,652.715,654.88,
- 656.99,659.047,661.052,663.007,664.963,666.771,668.6,
- 670.351,672.075,673.828,675.397,676.996,678.567,680.083,
- 681.589,683.023,684.457,685.825,687.18,688.499,689.789,
- 691.058,692.284,693.501,694.667,695.82,696.947,698.05,
- 699.129,700.186,701.221,702.234,703.226,704.198,705.158,
- 706.085,707.001,707.899,708.78,709.644,710.474,711.334,
- 712.145,712.943,713.727,714.505,715.266,716.015,716.751,
- 717.474,718.183,718.88,719.645,720.243,720.91,721.565,
- 722.211,722.851,723.473,724.094,724.697,725.296,725.886,
- 726.468,727.041,727.608,728.166,728.718,729.262,729.808,
- 730.337,730.856,731.374,731.883,732.386,732.884,733.373,
- 733.859,734.339,734.813};
-
- vector<double> tmp1(a1widthin,a1widthin+200);
- _a1runwidth.clear();
- std::transform(tmp1.begin(), tmp1.end(),
- back_inserter(_a1runwidth),
- [](double x){return x*GeV;});
-
- vector<double> tmp2(a1q2in,a1q2in+200);
- _a1runq2.clear();
- std::transform(tmp2.begin(), tmp2.end(),
- back_inserter(_a1runq2),
- [](double x){return x*GeV2;});
-
- _maxmass=ZERO;
- _maxcalc=ZERO;
+ _rhoF123masses = {0.773*GeV,1.370*GeV,1.750*GeV};
+ _rhoF123widths = {0.145*GeV,0.510*GeV,0.120*GeV};
+ _rhoF5masses = {0.773*GeV,1.500*GeV,1.750*GeV};
+ _rhoF5widths = {0.145*GeV,0.220*GeV,0.120*GeV};
}
void EtaPiPiDefaultCurrent::doinit() {
WeakCurrent::doinit();
// the particles we will use a lot
- tPDPtr a1(getParticleData(ParticleID::a_1minus)),
- k1(getParticleData(ParticleID::Kstar_1minus)),pi0(getParticleData(ParticleID::pi0)),
- piplus(getParticleData(ParticleID::piplus)),
- piminus(getParticleData(ParticleID::piminus));
+ tPDPtr piplus(getParticleData(ParticleID::piplus));
// masses for the running widths
_mpi=piplus->mass();
- _mK=getParticleData(ParticleID::Kminus)->mass();
- // the charged rho resonances
- tPDPtr rhoc[3]={getParticleData(-213),getParticleData(-100213),
- getParticleData(-30213)};
- // the charged K* resonances
- tPDPtr Kstarc[3]={getParticleData(-323),getParticleData(-100323),
- getParticleData(-30323)};
- if(!_a1parameters) {
- _a1mass=a1->mass();
- _a1width=a1->width();
- }
- // mass and width of the k_1
- if(!_k1parameters) {
- _k1mass=k1->mass();
- _k1width=k1->width();
- }
- // initialise the a_1 running width calculation
- inita1Width(-1);
- // rho parameters in the base classs
- tcPDPtr temp;
- unsigned int ix;
- if(_rhoparameters&&_rhoF123masses.size()<3) {
- ix = _rhoF123masses.size();
- _rhoF123masses.resize(3);_rhoF123widths.resize(3);
- for(;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF123masses[ix]=rhoc[ix]->mass();
- _rhoF123widths[ix]=rhoc[ix]->width();
- }
- }
- }
- else if(!_rhoparameters) {
- _rhoF123masses.resize(3);_rhoF123widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF123masses[ix]=rhoc[ix]->mass();
- _rhoF123widths[ix]=rhoc[ix]->width();
- }
- }
- }
- // K star parameters in the base class
- if(_kstarparameters&&_kstarF123masses.size()<3) {
- ix = _kstarF123masses.size();
- _kstarF123masses.resize(3);_kstarF123widths.resize(3);
- for(;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF123masses[ix]=Kstarc[ix]->mass();
- _kstarF123widths[ix]=Kstarc[ix]->width();
- }
- }
- }
- else if(!_kstarparameters) {
- _kstarF123masses.resize(3);_kstarF123widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF123masses[ix]=Kstarc[ix]->mass();
- _kstarF123widths[ix]=Kstarc[ix]->width();
- }
- }
- }
- // rho parameters here
- if(_rhoparameters&&_rhoF5masses.size()<3) {
- ix = _rhoF5masses.size();
- _rhoF5masses.resize(3);_rhoF5widths.resize(3);
- for(;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF5masses[ix]=rhoc[ix]->mass();
- _rhoF5widths[ix]=rhoc[ix]->width();
- }
- }
- }
- else if(!_rhoparameters) {
- _rhoF5masses.resize(3);_rhoF5widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF5masses[ix]=rhoc[ix]->mass();
- _rhoF5widths[ix]=rhoc[ix]->width();
- }
- }
- }
- // Kstar parameters here
- if(_kstarparameters&&_kstarF5widths.size()<3) {
- ix = _kstarF5masses.size();
- _kstarF5masses.resize(3);_kstarF5widths.resize(3);
- for(;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF5masses[ix]=Kstarc[ix]->mass();
- _kstarF5widths[ix]=Kstarc[ix]->width();
- }
- }
- }
- else if(!_kstarparameters) {
- _kstarF5masses.resize(3);_kstarF5widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF5masses[ix]=Kstarc[ix]->mass();
- _kstarF5widths[ix]=Kstarc[ix]->width();
- }
- }
- }
}
void EtaPiPiDefaultCurrent::persistentOutput(PersistentOStream & os) const {
- os << _rhoF123wgts << _kstarF123wgts << _rhoF5wgts << _kstarF5wgts
- << _rhoKstarwgt << ounit(_a1runwidth,GeV)<< ounit(_a1runq2,GeV2)
- << _initializea1
- << ounit(_a1mass,GeV)<< ounit(_a1width,GeV)<< ounit(_k1mass,GeV)
- << ounit(_k1width,GeV)<< ounit(_fpi,GeV) << ounit(_mpi,GeV)<< ounit(_mK,GeV)
- <<_rhoparameters << ounit(_rhoF123masses,GeV) << ounit(_rhoF5masses,GeV)
- << ounit(_rhoF123widths,GeV)
- << ounit(_rhoF5widths,GeV) << _kstarparameters << ounit(_kstarF123masses,GeV)
- <<ounit(_kstarF5masses,GeV)
- << ounit(_kstarF123widths,GeV) << ounit(_kstarF5widths,GeV) << _a1parameters
- << _k1parameters
- << _a1opt << ounit(_maxmass,GeV) << ounit(_maxcalc,GeV) << _a1runinter;
+ os << _rhoF123wgts << _rhoF5wgts << ounit(_fpi,GeV) << ounit(_mpi,GeV)
+ << ounit(_rhoF123masses,GeV) << ounit(_rhoF5masses,GeV)
+ << ounit(_rhoF123widths,GeV) << ounit(_rhoF5widths,GeV);
}
void EtaPiPiDefaultCurrent::persistentInput(PersistentIStream & is, int) {
- is >> _rhoF123wgts >> _kstarF123wgts >> _rhoF5wgts >> _kstarF5wgts
- >> _rhoKstarwgt >> iunit(_a1runwidth,GeV) >> iunit(_a1runq2,GeV2)
- >> _initializea1
- >> iunit(_a1mass,GeV) >> iunit(_a1width,GeV) >> iunit(_k1mass,GeV)
- >> iunit(_k1width,GeV) >> iunit(_fpi,GeV) >> iunit(_mpi,GeV) >> iunit(_mK,GeV)
- >>_rhoparameters >> iunit(_rhoF123masses,GeV) >> iunit(_rhoF5masses,GeV)
- >> iunit(_rhoF123widths,GeV)
- >> iunit(_rhoF5widths,GeV) >> _kstarparameters >> iunit(_kstarF123masses,GeV)
- >>iunit(_kstarF5masses,GeV)
- >> iunit(_kstarF123widths,GeV) >> iunit(_kstarF5widths,GeV) >> _a1parameters
- >> _k1parameters
- >> _a1opt >> iunit(_maxmass,GeV) >> iunit(_maxcalc,GeV) >> _a1runinter;
+ is >> _rhoF123wgts >> _rhoF5wgts >> iunit(_fpi,GeV) >> iunit(_mpi,GeV)
+ >> iunit(_rhoF123masses,GeV) >> iunit(_rhoF5masses,GeV)
+ >> iunit(_rhoF123widths,GeV) >> iunit(_rhoF5widths,GeV);
}
void EtaPiPiDefaultCurrent::Init() {
static ClassDocumentation<EtaPiPiDefaultCurrent> documentation
("The EtaPiPiDefaultCurrent class is designed to implement "
"the three meson decays of the tau, ie pi- pi- pi+, pi0 pi0 pi-, "
"K- pi- K+, K0 pi- Kbar0, K- pi0 K0,pi0 pi0 K-, K- pi- pi+, "
"pi- Kbar0 pi0, pi- pi0 eta. It uses the same currents as those in TAUOLA.",
"The three meson decays of the tau, ie pi- pi- pi+, pi0 pi0 pi-, "
"K- pi- K+, K0 pi- Kbar0, K- pi0 K0,pi0 pi0 K-, K- pi- pi+, "
"and pi- Kbar0 pi0, pi- pi0 eta "
"use the same currents as \\cite{Jadach:1993hs,Kuhn:1990ad,Decker:1992kj}.",
"%\\cite{Jadach:1993hs}\n"
"\\bibitem{Jadach:1993hs}\n"
" S.~Jadach, Z.~Was, R.~Decker and J.~H.~Kuhn,\n"
" %``The Tau Decay Library Tauola: Version 2.4,''\n"
" Comput.\\ Phys.\\ Commun.\\ {\\bf 76}, 361 (1993).\n"
" %%CITATION = CPHCB,76,361;%%\n"
"%\\cite{Kuhn:1990ad}\n"
"\\bibitem{Kuhn:1990ad}\n"
" J.~H.~Kuhn and A.~Santamaria,\n"
" %``Tau decays to pions,''\n"
" Z.\\ Phys.\\ C {\\bf 48}, 445 (1990).\n"
" %%CITATION = ZEPYA,C48,445;%%\n"
"%\\cite{Decker:1992kj}\n"
"\\bibitem{Decker:1992kj}\n"
" R.~Decker, E.~Mirkes, R.~Sauer and Z.~Was,\n"
" %``Tau decays into three pseudoscalar mesons,''\n"
" Z.\\ Phys.\\ C {\\bf 58}, 445 (1993).\n"
" %%CITATION = ZEPYA,C58,445;%%\n"
);
-
static ParVector<EtaPiPiDefaultCurrent,double> interfaceF123RhoWgt
("F123RhoWeight",
"The weights of the different rho resonances in the F1,2,3 form factor",
&EtaPiPiDefaultCurrent::_rhoF123wgts,
0, 0, 0, -1000, 1000, false, false, true);
- static ParVector<EtaPiPiDefaultCurrent,double> interfaceF123KstarWgt
- ("F123KstarWeight",
- "The weights of the different Kstar resonances in the F1,2,3 form factor",
- &EtaPiPiDefaultCurrent::_kstarF123wgts,
- 0, 0, 0, -1000, 1000, false, false, true);
-
static ParVector<EtaPiPiDefaultCurrent,double> interfaceF5RhoWgt
("F5RhoWeight",
"The weights of the different rho resonances in the F1,2,3 form factor",
&EtaPiPiDefaultCurrent::_rhoF5wgts,
0, 0, 0, -1000, 1000, false, false, true);
- static ParVector<EtaPiPiDefaultCurrent,double> interfaceF5KstarWgt
- ("F5KstarWeight",
- "The weights of the different Kstar resonances in the F1,2,3 form factor",
- &EtaPiPiDefaultCurrent::_kstarF5wgts,
- 0, 0, 0, -1000, 1000, false, false, true);
-
- static Parameter<EtaPiPiDefaultCurrent,double> interfaceRhoKstarWgt
- ("RhoKstarWgt",
- "The relative weights of the rho and K* in the F5 form factor",
- &EtaPiPiDefaultCurrent::_rhoKstarwgt, -0.2, -10., 10.,
- false, false, false);
-
- static Switch<EtaPiPiDefaultCurrent,bool> interfaceInitializea1
- ("Initializea1",
- "Initialise the calculation of the a_1 running width",
- &EtaPiPiDefaultCurrent::_initializea1, false, false, false);
- static SwitchOption interfaceInitializea1Initialization
- (interfaceInitializea1,
- "Yes",
- "Initialize the calculation",
- true);
- static SwitchOption interfaceInitializea1NoInitialization
- (interfaceInitializea1,
- "No",
- "Use the default values",
- false);
-
- static Switch<EtaPiPiDefaultCurrent,bool> interfaceRhoParameters
- ("RhoParameters",
- "Use local values of the rho meson masses and widths",
- &EtaPiPiDefaultCurrent::_rhoparameters, true, false, false);
- static SwitchOption interfaceRhoParameterstrue
- (interfaceRhoParameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfaceRhoParametersParticleData
- (interfaceRhoParameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<EtaPiPiDefaultCurrent,bool> interfaceKstarParameters
- ("KstarParameters",
- "Use local values of the rho meson masses and widths",
- &EtaPiPiDefaultCurrent::_kstarparameters, true, false, false);
- static SwitchOption interfaceKstarParameterstrue
- (interfaceKstarParameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfaceKstarParametersParticleData
- (interfaceKstarParameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<EtaPiPiDefaultCurrent,bool> interfacea1Parameters
- ("a1Parameters",
- "Use local values of the rho meson masses and widths",
- &EtaPiPiDefaultCurrent::_a1parameters, true, false, false);
- static SwitchOption interfacea1Parameterstrue
- (interfacea1Parameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfacea1ParametersParticleData
- (interfacea1Parameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<EtaPiPiDefaultCurrent,bool> interfaceK1Parameters
- ("K1Parameters",
- "Use local values of the rho meson masses and widths",
- &EtaPiPiDefaultCurrent::_k1parameters, true, false, false);
- static SwitchOption interfaceK1Parameterstrue
- (interfaceK1Parameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfaceK1ParametersParticleData
- (interfaceK1Parameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<EtaPiPiDefaultCurrent,bool> interfacea1WidthOption
- ("a1WidthOption",
- "Option for the treatment of the a1 width",
- &EtaPiPiDefaultCurrent::_a1opt, true, false, false);
- static SwitchOption interfacea1WidthOptionLocal
- (interfacea1WidthOption,
- "Local",
- "Use a calculation of the running width based on the parameters as"
- " interpolation table.",
- true);
- static SwitchOption interfacea1WidthOptionParam
- (interfacea1WidthOption,
- "Kuhn",
- "Use the parameterization of Kuhn and Santamaria for default parameters."
- " This should only be used for testing vs TAUOLA",
- false);
-
- static ParVector<EtaPiPiDefaultCurrent,Energy> interfacea1RunningWidth
- ("a1RunningWidth",
- "The values of the a_1 width for interpolation to giving the running width.",
- &EtaPiPiDefaultCurrent::_a1runwidth, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<EtaPiPiDefaultCurrent,Energy2> interfacea1RunningQ2
- ("a1RunningQ2",
- "The values of the q^2 for interpolation to giving the running width.",
- &EtaPiPiDefaultCurrent::_a1runq2, GeV2, -1, 1.0*GeV2, ZERO, 10.0*GeV2,
- false, false, true);
-
- static Parameter<EtaPiPiDefaultCurrent,Energy> interfaceA1Width
- ("A1Width",
- "The a_1 width if using local values.",
- &EtaPiPiDefaultCurrent::_a1width, GeV, 0.599*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
- static Parameter<EtaPiPiDefaultCurrent,Energy> interfaceA1Mass
- ("A1Mass",
- "The a_1 mass if using local values.",
- &EtaPiPiDefaultCurrent::_a1mass, GeV, 1.251*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
- static Parameter<EtaPiPiDefaultCurrent,Energy> interfaceK1Width
- ("K1Width",
- "The K_1 width if using local values.",
- &EtaPiPiDefaultCurrent::_k1width, GeV, 0.174*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
- static Parameter<EtaPiPiDefaultCurrent,Energy> interfaceK1Mass
- ("K1Mass",
- "The K_1 mass if using local values.",
- &EtaPiPiDefaultCurrent::_k1mass, GeV, 1.402*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
static ParVector<EtaPiPiDefaultCurrent,Energy> interfacerhoF123masses
("rhoF123masses",
"The masses for the rho resonances if used local values",
&EtaPiPiDefaultCurrent::_rhoF123masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
false, false, true);
static ParVector<EtaPiPiDefaultCurrent,Energy> interfacerhoF123widths
("rhoF123widths",
"The widths for the rho resonances if used local values",
&EtaPiPiDefaultCurrent::_rhoF123widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
false, false, true);
static ParVector<EtaPiPiDefaultCurrent,Energy> interfacerhoF5masses
("rhoF5masses",
"The masses for the rho resonances if used local values",
&EtaPiPiDefaultCurrent::_rhoF5masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
false, false, true);
static ParVector<EtaPiPiDefaultCurrent,Energy> interfacerhoF5widths
("rhoF5widths",
"The widths for the rho resonances if used local values",
&EtaPiPiDefaultCurrent::_rhoF5widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
false, false, true);
- static ParVector<EtaPiPiDefaultCurrent,Energy> interfaceKstarF123masses
- ("KstarF123masses",
- "The masses for the Kstar resonances if used local values",
- &EtaPiPiDefaultCurrent::_kstarF123masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<EtaPiPiDefaultCurrent,Energy> interfaceKstarF123widths
- ("KstarF123widths",
- "The widths for the Kstar resonances if used local values",
- &EtaPiPiDefaultCurrent::_kstarF123widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<EtaPiPiDefaultCurrent,Energy> interfaceKstarF5masses
- ("KstarF5masses",
- "The masses for the Kstar resonances if used local values",
- &EtaPiPiDefaultCurrent::_kstarF5masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<EtaPiPiDefaultCurrent,Energy> interfaceKstarF5widths
- ("KstarF5widths",
- "The widths for the Kstar resonances if used local values",
- &EtaPiPiDefaultCurrent::_kstarF5widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
static Parameter<EtaPiPiDefaultCurrent,Energy> interfaceFPi
("FPi",
"The pion decay constant",
&EtaPiPiDefaultCurrent::_fpi, MeV, 92.4*MeV, ZERO, 200.0*MeV,
false, false, true);
}
-
-// modes handled by this class
-bool EtaPiPiDefaultCurrent::acceptMode(int imode) const {
- return imode>=0&&imode<=8;
-}
-
-// calculate the form-factors
-EtaPiPiDefaultCurrent::FormFactors
-EtaPiPiDefaultCurrent::calculateFormFactors(const int ichan, const int imode,
- Energy2 q2, Energy2 s1,
- Energy2 s2, Energy2 s3) const {
- useMe();
- Complex F1, F2, F3, F4, F5;
- F1 = F2 = F3 = F4 = F5 = 0.0;
- // calculate the pi- pi- pi+ factor
- if(imode==0) {
- Complex a1fact(a1BreitWigner(q2)*2./3.);
- if(ichan<0) {
- F1= a1fact*BrhoF123(s1,-1);
- F2 =-a1fact*BrhoF123(s2,-1);
- }
- else if(ichan%2==0) F1 = a1fact*BrhoF123(s1, ichan/2);
- else if(ichan%2==1) F2 =-a1fact*BrhoF123(s2,(ichan-1)/2);
- }
- // calculate the pi0 pi0 pi- factor
- else if(imode==1) {
- Complex a1fact(a1BreitWigner(q2)*2./3.);
- if(ichan<0) {
- F1 = a1fact*BrhoF123(s1,-1);
- F2 =-a1fact*BrhoF123(s2,-1);
- }
- else if(ichan%2==0) F1 = a1fact*BrhoF123(s1, ichan/2);
- else if(ichan%2==1) F2 =-a1fact*BrhoF123(s2,(ichan-1)/2);
- }
- // calculate the K- pi - K+ factor
- else if(imode==2) {
- Complex a1fact(a1BreitWigner(q2)*sqrt(2.)/3.);
- if(ichan<0) {
- F1 =-a1fact*BKstarF123(s1,-1);
- F2 = a1fact*BrhoF123(s2,-1);
- F5 = BrhoF5(q2,-1)*FKrho(s1,s2,-1)*sqrt(2.);
- }
- else if(ichan%8==0) F1 =-a1fact*BKstarF123(s1,ichan/8);
- else if(ichan%8==1) F2 = a1fact*BrhoF123(s2,(ichan-1)/8);
- else if(ichan%8>=2) F5 = BrhoF5(q2,ichan/8)*FKrho(s1,s2,(ichan-2)%8)*sqrt(2.);
- }
- // calculate the K0 pi- K0bar
- else if(imode==3) {
- Complex a1fact(a1BreitWigner(q2)*sqrt(2.)/3.);
- if(ichan<0) {
- F1 =-a1fact*BKstarF123(s1,-1);
- F2 = a1fact*BrhoF123(s2,-1);
- F5 =-BrhoF5(q2,-1)*FKrho(s1,s2,-1)*sqrt(2.);
- }
- else if(ichan%8==0) F1 = -a1fact*BKstarF123(s1,ichan/8);
- else if(ichan%8==1) F2 = a1fact*BrhoF123(s2,(ichan-1)/8);
- else if(ichan%8>=2) F5 = -BrhoF5(q2,ichan/8)*FKrho(s1,s2,(ichan-2)%8)*sqrt(2.);
- }
- // calculate the K- pi0 k0
- else if(imode==4) {
- Complex a1fact(a1BreitWigner(q2));
- if(ichan<0){F2 =-a1fact*BrhoF123(s2,-1);}
- else{F2 =-a1fact*BrhoF123(s2,ichan);}
- }
- // calculate the pi0 pi0 K-
- else if(imode==5) {
- Complex K1fact(K1BreitWigner(q2)/6.);
- if(ichan<0) {
- F1 = K1fact*BKstarF123(s1,-1);
- F2 =-K1fact*BKstarF123(s2,-1);
- }
- else if(ichan%2==0) F1 = K1fact*BKstarF123(s1,ichan/2);
- else F2 =-K1fact*BKstarF123(s2,(ichan-1)/2);
- }
- // calculate the K- pi- pi+
- else if(imode==6) {
- Complex K1fact(K1BreitWigner(q2)*sqrt(2.)/3.);
- if(ichan<0) {
- F1 =-K1fact*BrhoF123(s1,-1);
- F2 = K1fact*BKstarF123(s2,-1);
- F5 =-BKstarF123(q2,-1)*FKrho(s2,s1,-1)*sqrt(2.);
- }
- else if(ichan%8==0) F1 =-K1fact*BrhoF123(s1,ichan/8);
- else if(ichan%8==1) F2 = K1fact*BKstarF123(s2,(ichan-1)/8);
- else F5 = -BKstarF123(q2,ichan/8)*FKrho(s2,s1,(ichan-2)%8)*sqrt(2.);
- }
- // calculate the pi- K0bar pi0
- else if(imode==7) {
- Complex K1fact(K1BreitWigner(q2));
- if(ichan<0) {
- F2 =-K1fact*BrhoF123(s2,-1);
- F5 =-2.*BKstarF123(q2,-1)*FKrho(s1,s2,-1);
- }
- else if(ichan%7==0) F2 =-K1fact*BrhoF123(s2,ichan/7);
- else F5 =-2.*BKstarF123(q2,ichan/7)*FKrho(s1,s2,(ichan-1)%7);
- }
- // calculate the pi- pi0 eta
- else if(imode==8) {
- if(ichan<0) F5 = BrhoF5(q2, -1)*BrhoF123(s3, -1)*sqrt(2./3.);
- else F5 = BrhoF5(q2,ichan/3)*BrhoF123(s3,ichan%3)*sqrt(2./3.);
- }
- // multiply by the prefactors
- using Constants::twopi;
- return FormFactors(F1/_fpi,
- F2/_fpi,
- F3/_fpi,
- F4/_fpi,
- -F5/sqr(twopi)/pow<3,1>(_fpi)
- );
-}
// complete the construction of the decay mode for integration
bool EtaPiPiDefaultCurrent::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 ) {
+ IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
+ unsigned int imode,PhaseSpaceModePtr mode,
+ unsigned int iloc,int ires,
+ PhaseSpaceChannel phase, Energy upp ) {
+ // check the charge
+ if((imode==0 && abs(icharge)!=3) ||
+ (imode>0 && icharge !=0)) return false;
+ // check the total isospin
+ if(Itotal!=IsoSpin::IUnknown) {
+ if(Itotal!=IsoSpin::IOne) return false;
+ }
+ // check I_3
+ if(i3!=IsoSpin::I3Unknown) {
+ switch(i3) {
+ case IsoSpin::I3Zero:
+ if(imode==0) return false;
+ break;
+ case IsoSpin::I3One:
+ if(imode==1 || icharge ==-3) return false;
+ break;
+ case IsoSpin::I3MinusOne:
+ if(imode==1 || icharge ==3) return false;
+ default:
+ return false;
+ }
+ }
+ // make sure that the decays are kinematically allowed
int iq(0),ia(0);
- if(!acceptMode(imode)) return false;
+ tPDVector part = particles(icharge,imode,iq,ia);
tPDVector extpart(particles(1,imode,iq,ia));
Energy min(ZERO);
for(unsigned int ix=0;ix<extpart.size();++ix) min+=extpart[ix]->massMin();
if(min>upp) return false;
- // the particles we will use a lot
- tPDPtr a1,k1;
- if(icharge==-3) {
- a1=getParticleData(ParticleID::a_1minus);
- k1=getParticleData(ParticleID::Kstar_1minus);
- }
- else if(icharge==3) {
- a1=getParticleData(ParticleID::a_1plus);
- k1=getParticleData(ParticleID::Kstar_1plus);
+ // set up the resonances
+ tPDPtr res[3];
+ if(icharge==0) {
+ res[0] =getParticleData(113);
+ res[1] =getParticleData(100113);
+ res[2] =getParticleData(30113);
}
else {
- return false;
- }
- _maxmass=max(_maxmass,upp);
- // the rho0 resonances
- tPDPtr rho0[3] = { getParticleData(113), getParticleData(100113), getParticleData(30113)};
- tPDPtr rhoc[3] = {getParticleData(-213),getParticleData(-100213),getParticleData(-30213)};
- tPDPtr Kstar0[3] = { getParticleData(313), getParticleData(100313), getParticleData(30313)};
- tPDPtr Kstarc[3] = {getParticleData(-323),getParticleData(-100323),getParticleData(-30323)};
- if(icharge==3) {
- for(unsigned int ix=0;ix<3;++ix) {
- rhoc [ix] = rhoc[ix]->CC();
- Kstar0[ix] = Kstar0[ix]->CC();
- Kstarc[ix] = Kstarc[ix]->CC();
- }
- }
- if(imode==0) {
- // channels for pi- pi- pi+
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rho0[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,rho0[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==1) {
- // channels for pi0 pi0 pi-
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rhoc[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,rhoc[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rhoc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==2) {
- // channels for K- pi- K+
- for(unsigned int ix=0;ix<3;++ix) {
- if(Kstar0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,Kstar0[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!rhoc[ix]) continue;
- if(Kstar0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,Kstar0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+2,ires+1,rho0[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
+ res[0] =getParticleData(213);
+ res[1] =getParticleData(100213);
+ res[2] =getParticleData(30213);
+ if(icharge==-3) {
+ for(unsigned int ix=0;ix<3;++ix) {
+ if(res[ix]&&res[ix]->CC()) res[ix]=res[ix]->CC();
}
}
}
- else if(imode==3) {
- // channels for K0 pi- K0bar
- for(unsigned int ix=0;ix<3;++ix) {
- if(Kstarc[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,Kstarc[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!rhoc[ix]) continue;
- if(Kstarc[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,Kstarc[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+2,ires+1,rho0[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
+ // channels for pi- pi0 eta
+ for(unsigned int ix=0;ix<3;++ix) {
+ if(resonance && resonance != res[ix]) continue;
+ for(unsigned int iy=0;iy<3;++iy) {
+ mode->addChannel((PhaseSpaceChannel(phase),ires,res[ix],ires+1,iloc+3,ires+1,res[iy],
+ ires+2,iloc+1,ires+2,iloc+2));
}
}
- else if(imode==4) {
- // channels for K- pi0 K0
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rhoc[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rhoc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==5) {
- // channels for pi0 pi0 K-
- for(unsigned int ix=0;ix<3;++ix) {
- if(!Kstarc[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+1,ires+1,Kstarc[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+2,ires+1,Kstarc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==6) {
- // channels for K- pi- pi+
- for(unsigned int ix=0;ix<3;++ix) {
- if(rho0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+1,ires+1,rho0[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(Kstar0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+2,ires+1,Kstar0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!Kstarc[ix]) continue;
- if(rho0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+1,ires+1,rho0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(Kstar0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+2,ires+1,Kstar0[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- }
- }
- else if(imode==7) {
- // channels for pi- kbar0 pi0
- for(unsigned int ix=0;ix<3;++ix) {
- if(rhoc[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+2,ires+1,rhoc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!Kstarc[ix]) continue;
- if(Kstar0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+1,ires+1,Kstar0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rhoc[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+2,ires+1,rhoc[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- }
- }
- else if(imode==8) {
- // channels for pi- pi0 eta
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rhoc[ix]) continue;
- for(unsigned int iy=0;iy<3;++iy) {
- if(!rho0[iy]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,rho0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- }
- }
- if(_rhoparameters) {
- if(imode!=8) {
- for(unsigned int ix=0;ix<_rhoF123masses.size();++ix) {
- if(rhoc[ix]) mode->resetIntermediate(rhoc[ix],_rhoF123masses[ix],
- _rhoF123widths[ix]);
- if(rho0[ix]) mode->resetIntermediate(rho0[ix],_rhoF123masses[ix],
- _rhoF123widths[ix]);
- }
- }
- else {
- for(unsigned int ix=0;ix<_rhoF5masses.size();++ix) {
- if(rhoc[ix]) mode->resetIntermediate(rhoc[ix],_rhoF5masses[ix],
- _rhoF5widths[ix]);
- if(rho0[ix]) mode->resetIntermediate(rho0[ix],_rhoF5masses[ix],
- _rhoF5widths[ix]);
- }
- }
- }
- // K star parameters in the base class
- if(_kstarparameters) {
- for(unsigned int ix=0;ix<_kstarF123masses.size();++ix) {
- if(Kstarc[ix]) mode->resetIntermediate(Kstarc[ix],_kstarF123masses[ix],
- _kstarF123widths[ix]);
- if(Kstar0[ix]) mode->resetIntermediate(Kstar0[ix],_kstarF123masses[ix],
- _kstarF123widths[ix]);
- }
- }
+ // reset the rho masses
+ for(unsigned int ix=0;ix<_rhoF5masses.size();++ix)
+ mode->resetIntermediate(res[ix],_rhoF5masses[ix],_rhoF5widths[ix]);
return true;
}
-// initialisation of the a_1 width
-// (iopt=-1 initialises, iopt=0 starts the interpolation)
-void EtaPiPiDefaultCurrent::inita1Width(int iopt) {
- if(iopt==-1) {
- _maxcalc=_maxmass;
- if(!_initializea1||_maxmass==ZERO) return;
- // parameters for the table of values
- Energy2 step(sqr(_maxcalc)/199.);
- // integrator to perform the integral
- vector<double> inweights;inweights.push_back(0.5);inweights.push_back(0.5);
- vector<int> intype;intype.push_back(2);intype.push_back(3);
- Energy mrho(getParticleData(ParticleID::rhoplus)->mass()),
- wrho(getParticleData(ParticleID::rhoplus)->width());
- vector<Energy> inmass(2,mrho),inwidth(2,wrho);
- vector<double> inpow(2,0.0);
- ThreeBodyAllOnCalculator<EtaPiPiDefaultCurrent>
- widthgen(inweights,intype,inmass,inwidth,inpow,*this,0,_mpi,_mpi,_mpi);
- // normalisation constant to give physical width if on shell
- double a1const(_a1width/(widthgen.partialWidth(sqr(_a1mass))));
- // loop to give the values
- _a1runq2.clear(); _a1runwidth.clear();
- for(Energy2 moff2(ZERO); moff2<=sqr(_maxcalc); moff2+=step) {
- _a1runwidth.push_back(widthgen.partialWidth(moff2)*a1const);
- _a1runq2.push_back(moff2);
- }
- }
- // set up the interpolator
- else if(iopt==0) {
- _a1runinter = make_InterpolatorPtr(_a1runwidth,_a1runq2,3);
- }
-}
-
void EtaPiPiDefaultCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header) output << "update decayers set parameters=\"";
if(create) output << "create Herwig::EtaPiPiDefaultCurrent "
<< name() << " HwWeakCurrents.so\n";
+ output << "newdef " << name() << ":FPi " << _fpi/MeV << "\n";
for(unsigned int ix=0;ix<_rhoF123wgts.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":F123RhoWeight " << ix << " " << _rhoF123wgts[ix] << "\n";
}
- for(unsigned int ix=0;ix<_kstarF123wgts.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":F123KstarWeight " << ix << " "
- << _kstarF123wgts[ix] << "\n";
- }
for(unsigned int ix=0;ix<_rhoF5wgts.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":F5RhoWeight " << ix << " " << _rhoF5wgts[ix] << "\n";
}
- for(unsigned int ix=0;ix<_kstarF5wgts.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":F5KstarWeight " << ix << " " << _kstarF5wgts[ix] << "\n";
- }
- output << "newdef " << name() << ":RhoKstarWgt " << _rhoKstarwgt << "\n";
- output << "newdef " << name() << ":Initializea1 " << _initializea1 << "\n";
- output << "newdef " << name() << ":RhoParameters " << _rhoparameters << "\n";
- output << "newdef " << name() << ":KstarParameters " << _kstarparameters << "\n";
- output << "newdef " << name() << ":a1Parameters " << _a1parameters << "\n";
- output << "newdef " << name() << ":K1Parameters " << _k1parameters << "\n";
- output << "newdef " << name() << ":a1WidthOption " << _a1opt << "\n";
- for(unsigned int ix=0;ix<_a1runwidth.size();++ix) {
- output << "newdef " << name() << ":a1RunningWidth " << ix
- << " " << _a1runwidth[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_a1runq2.size();++ix) {
- output << "newdef " << name() << ":a1RunningQ2 " << ix
- << " " << _a1runq2[ix]/GeV2 << "\n";
- }
- output << "newdef " << name() << ":A1Width " << _a1width/GeV << "\n";
- output << "newdef " << name() << ":A1Mass " << _a1mass/GeV << "\n";
- output << "newdef " << name() << ":K1Width " << _k1width/GeV << "\n";
- output << "newdef " << name() << ":K1Mass " << _k1mass/GeV << "\n";
- output << "newdef " << name() << ":FPi " << _fpi/MeV << "\n";
for(unsigned int ix=0;ix<_rhoF123masses.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":rhoF123masses " << ix
<< " " << _rhoF123masses[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<_rhoF123widths.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":rhoF123widths " << ix << " "
<< _rhoF123widths[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<_rhoF5masses.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":rhoF5masses " << ix << " "
<< _rhoF5masses[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<_rhoF5widths.size();++ix) {
if(ix<3) output << "newdef ";
else output << "insert ";
output << name() << ":rhoF5widths " << ix << " "
<< _rhoF5widths[ix]/GeV << "\n";
}
- for(unsigned int ix=0;ix<_kstarF123masses.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF123masses " << ix << " "
- << _kstarF123masses[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF123widths.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF123widths " << ix << " "
- << _kstarF123widths[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF5masses.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF5masses " << ix << " "
- << _kstarF5masses[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF5widths.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF5widths " << ix << " "
- << _kstarF5widths[ix]/GeV << "\n";
- }
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
-void EtaPiPiDefaultCurrent::doinitrun() {
- // set up the running a_1 width
- inita1Width(0);
- WeakCurrent::doinitrun();
-}
-
-void EtaPiPiDefaultCurrent::doupdate() {
- WeakCurrent::doupdate();
- // update running width if needed
- if ( !touched() ) return;
- if(_maxmass!=_maxcalc) inita1Width(-1);
-}
-
-Complex EtaPiPiDefaultCurrent::rhoKBreitWigner(Energy2 q2,unsigned int itype,
- unsigned int ires) const {
- Energy q(sqrt(q2)),mass,width,mout[2]={_mpi,_mpi};
- // get the mass and width of the requested resonance
- if(itype==0) {
- mass=_rhoF123masses[ires];
- width=_rhoF123widths[ires];
- }
- else if(itype==1) {
- mass=_rhoF5masses[ires];
- width=_rhoF5widths[ires];
- }
- else if(itype==2) {
- mass=_kstarF123masses[ires];
- width=_kstarF123widths[ires];
- }
- else if(itype==3) {
- mass=_kstarF5masses[ires];
- width=_kstarF5widths[ires];
- }
- else {
- return 0.;
- }
- // calculate the momenta for the running widths
- if(itype>1) mout[0]=_mK;
- Energy pcm0(Kinematics::pstarTwoBodyDecay(mass,mout[0],mout[1]));
- Energy pcm(ZERO);
- if(mout[0]+mout[1]<q){pcm=Kinematics::pstarTwoBodyDecay(q,mout[0],mout[1]);}
- double ratio = Math::Pow<3>(pcm/pcm0);
- Energy gamrun(width*mass*ratio/q);
- Complex ii(0.,1.);
- complex<Energy2> denom(q2-mass*mass+ii*mass*gamrun), numer(-mass*mass);
- return numer/denom;
-}
-
-double EtaPiPiDefaultCurrent::
-threeBodyMatrixElement(const int , const Energy2 q2,
- const Energy2 s3, const Energy2 s2,
- const Energy2 s1, const Energy ,
- const Energy , const Energy ) const {
- Energy2 mpi2(sqr(_mpi));
- Complex propb(BrhoF123(s1,-1)),propa(BrhoF123(s2,-1));
- // the matrix element
- Energy2 output(ZERO);
- // first resonance
- output += ((s1-4.*mpi2) + 0.25*(s3-s2)*(s3-s2)/q2) * real(propb*conj(propb));
- // second resonance
- output += ((s2-4.*mpi2) + 0.25*(s3-s1)*(s3-s1)/q2) * real(propa*conj(propa));
- // the interference term
- output += (0.5*q2-s3-0.5*mpi2+0.25*(s3-s2)*(s3-s1)/q2)*real(propa*conj(propb)+
- propb*conj(propa));
- return output/sqr(_rhoF123masses[0]);
-}
-
-
// the hadronic currents
vector<LorentzPolarizationVectorE>
EtaPiPiDefaultCurrent::current(tcPDPtr resonance,
- IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
- const int imode, const int ichan, Energy & scale,
- const tPDVector & ,
- const vector<Lorentz5Momentum> & momenta,
- DecayIntegrator::MEOption) const {
+ IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
+ const int imode, const int ichan, Energy & scale,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ DecayIntegrator::MEOption) const {
+ useMe();
+ // check the isospin
+ if(Itotal!=IsoSpin::IUnknown && Itotal!=IsoSpin::IOne)
+ return vector<LorentzPolarizationVectorE>();
+ int icharge = outgoing[0]->iCharge()+outgoing[1]->iCharge();
+ // check I_3
+ if(i3!=IsoSpin::I3Unknown) {
+ switch(i3) {
+ case IsoSpin::I3Zero:
+ if(imode==0) return vector<LorentzPolarizationVectorE>();
+ break;
+ case IsoSpin::I3One:
+ if(imode==1 || icharge ==-3) return vector<LorentzPolarizationVectorE>();
+ break;
+ case IsoSpin::I3MinusOne:
+ if(imode==1 || icharge ==3) return vector<LorentzPolarizationVectorE>();
+ default:
+ return vector<LorentzPolarizationVectorE>();
+ }
+ }
// calculate q2,s1,s2,s3
- Lorentz5Momentum q;
- for(unsigned int ix=0;ix<momenta.size();++ix)
- q+=momenta[ix];
+ Lorentz5Momentum q = momenta[0] + momenta[1] + momenta[2];
q.rescaleMass();
scale=q.mass();
Energy2 q2=q.mass2();
- Energy2 s1 = (momenta[1]+momenta[2]).m2();
- Energy2 s2 = (momenta[0]+momenta[2]).m2();
Energy2 s3 = (momenta[0]+momenta[1]).m2();
- FormFactors F = calculateFormFactors(ichan,imode,q2,s1,s2,s3);
- //if(inpart.id()==ParticleID::tauplus){F.F5=conj(F.F5);}
- // the first three form-factors
- LorentzPolarizationVector vect;
- vect = (F.F2-F.F1)*momenta[2]
- +(F.F1-F.F3)*momenta[1]
- +(F.F3-F.F2)*momenta[0];
- // multiply by the transverse projection operator
- complex<InvEnergy> dot=(vect*q)/q2;
- // scalar and parity violating terms
- vect += (F.F4-dot)*q;
- if(F.F5!=complex<InvEnergy3>())
- vect += Complex(0.,1.)*F.F5*Helicity::epsilon(momenta[0],
- momenta[1],
- momenta[2]);
+ // the form factor
+ Complex F5(0.);
+ int ires1(-1),ires2(-1);
+ if(ichan>=0) {
+ ires1 = ichan/3;
+ ires2 = ichan%3;
+ }
+ else {
+ if(resonance) {
+ switch(resonance->id()/1000) {
+ case 0:
+ ires1 = 0;
+ break;
+ case 100:
+ ires1 = 1;
+ break;
+ case 30 :
+ ires1 = 2;
+ break;
+ default:
+ assert(false);
+ }
+ }
+ }
+ F5 = BrhoF5(q2,ires1)*BrhoF123(s3,ires2)*sqrt(2./3.);
+ // constant the current
+ LorentzPolarizationVector vect = -Complex(0.,1.)*F5/sqr(Constants::twopi)/pow<3,1>(_fpi)*
+ Helicity::epsilon(momenta[0],momenta[1],momenta[2]);
// factor to get dimensions correct
return vector<LorentzPolarizationVectorE>(1,q.mass()*vect);
}
bool EtaPiPiDefaultCurrent::accept(vector<int> id) {
- int npip(0),npim(0),nkp(0),nkm(0),
- npi0(0),nk0(0),nk0bar(0),neta(0),nks(0),nkl(0);
+ // check there are only three particles
+ if(id.size()!=3) return false;
+ unsigned int npip(0),npim(0),npi0(0),neta(0);
for(unsigned int ix=0;ix<id.size();++ix) {
if(id[ix]==ParticleID::piplus) ++npip;
else if(id[ix]==ParticleID::piminus) ++npim;
- else if(id[ix]==ParticleID::Kplus) ++nkp;
- else if(id[ix]==ParticleID::Kminus) ++nkm;
else if(id[ix]==ParticleID::pi0) ++npi0;
- else if(id[ix]==ParticleID::K0) ++nk0;
- else if(id[ix]==ParticleID::Kbar0) ++nk0bar;
else if(id[ix]==ParticleID::eta) ++neta;
- else if(id[ix]==ParticleID::K_S0) ++nks;
- else if(id[ix]==ParticleID::K_L0) ++nkl;
}
- int imode(-1);
- if( (npip==2&&npim==1) || (npim==2&&npip==1) ) imode= 0;
- else if( (npip==1&&npi0==2) || (npim==1&&npi0==2) ) imode= 1;
- else if( (nkp==1&&nkm==1&&npip==1) ||
- (nkp==1&&nkm==1&&npim==1)) imode= 2;
- else if( (nk0==1&&nk0bar==1&&npip==1) ||
- (nk0==1&&nk0bar==1&&npim==1)) imode= 3;
- else if( (nkp==1&&nk0bar==1&&npi0==1) ||
- (nkm==1&&npi0==1&&nk0==1)) imode= 4;
- else if( (nkp==1&&npi0==2) || (npi0==2&&nkm==1) ) imode= 5;
- else if( (npip==1&&npim==1&&nkp==1) ||
- (nkm==1&&npim==1&&npip==1) ) imode= 6;
- else if( (nk0==1&&npip==1&&npi0==1) ||
- (npim==1&&nk0bar==1&&npi0==1)) imode= 7;
- else if( (npip==1&&npi0==1&&neta==1) ||
- (npim==1&&npi0==1&&neta==1)) imode= 8;
- else if( nks==2 && (npip==1||npim==1) ) imode= 9;
- else if( nkl==2 && (npip==1||npim==1) ) imode=10;
- else if( nks==1&&nkl==1 && (npip==1||npim==1) ) imode=11;
- return imode==-1 ? false : acceptMode(imode);
+ if( (npip==1&&npim==1&&neta==1) ||
+ (npi0==1&&npim+npip==1&&neta==1))
+ return true;
+ else
+ return false;
}
unsigned int EtaPiPiDefaultCurrent::decayMode(vector<int> id) {
- int npip(0),npim(0),nkp(0),nkm(0),
- npi0(0),nk0(0),nk0bar(0),neta(0),nks(0),nkl(0);
+ unsigned int npi(0);
for(unsigned int ix=0;ix<id.size();++ix) {
- if(id[ix]==ParticleID::piplus) ++npip;
- else if(id[ix]==ParticleID::piminus) ++npim;
- else if(id[ix]==ParticleID::Kplus) ++nkp;
- else if(id[ix]==ParticleID::Kminus) ++nkm;
- else if(id[ix]==ParticleID::pi0) ++npi0;
- else if(id[ix]==ParticleID::K0) ++nk0;
- else if(id[ix]==ParticleID::Kbar0) ++nk0bar;
- else if(id[ix]==ParticleID::eta) ++neta;
- else if(id[ix]==ParticleID::K_S0) ++nks;
- else if(id[ix]==ParticleID::K_L0) ++nkl;
+ if(abs(id[ix])==ParticleID::piplus) ++npi;
}
- int imode(-1);
- if( (npip==2&&npim==1) || (npim==2&&npip==1) ) imode= 0;
- else if( (npip==1&&npi0==2) || (npim==1&&npi0==2) ) imode= 1;
- else if( (nkp==1&&nkm==1&&npip==1) ||
- (nkp==1&&nkm==1&&npim==1)) imode= 2;
- else if( (nk0==1&&nk0bar==1&&npip==1) ||
- (nk0==1&&nk0bar==1&&npim==1)) imode= 3;
- else if( (nkp==1&&nk0bar==1&&npi0==1) ||
- (nkm==1&&npi0==1&&nk0==1)) imode= 4;
- else if( (nkp==1&&npi0==2) || (npi0==2&&nkm==1) ) imode= 5;
- else if( (npip==1&&npim==1&&nkp==1) ||
- (nkm==1&&npim==1&&npip==1) ) imode= 6;
- else if( (nk0==1&&npip==1&&npi0==1) ||
- (npim==1&&nk0bar==1&&npi0==1)) imode= 7;
- else if( (npip==1&&npi0==1&&neta==1) ||
- (npim==1&&npi0==1&&neta==1)) imode= 8;
- else if( nks==2 && (npip==1||npim==1) ) imode= 9;
- else if( nkl==2 && (npip==1||npim==1) ) imode=10;
- else if( nks==1&&nkl==1 && (npip==1||npim==1) ) imode=11;
- return imode;
+ if(npi==2) return 1;
+ else return 0;
}
tPDVector EtaPiPiDefaultCurrent::particles(int icharge, unsigned int imode,int,int) {
- tPDVector extpart(3);
+ tPDVector output(3);
+ output[0]=getParticleData(ParticleID::piplus);
+ output[2]=getParticleData(ParticleID::eta);
if(imode==0) {
- extpart[0]=getParticleData(ParticleID::piminus);
- extpart[1]=getParticleData(ParticleID::piminus);
- extpart[2]=getParticleData(ParticleID::piplus);
- }
- else if(imode==1) {
- extpart[0]=getParticleData(ParticleID::pi0);
- extpart[1]=getParticleData(ParticleID::pi0);
- extpart[2]=getParticleData(ParticleID::piminus);
- }
- else if(imode==2) {
- extpart[0]=getParticleData(ParticleID::Kminus);
- extpart[1]=getParticleData(ParticleID::piminus);
- extpart[2]=getParticleData(ParticleID::Kplus);
- }
- else if(imode==3) {
- extpart[0]=getParticleData(ParticleID::K0);
- extpart[1]=getParticleData(ParticleID::piminus);
- extpart[2]=getParticleData(ParticleID::Kbar0);
- }
- else if(imode==4) {
- extpart[0]=getParticleData(ParticleID::Kminus);
- extpart[1]=getParticleData(ParticleID::pi0);
- extpart[2]=getParticleData(ParticleID::K0);
- }
- else if(imode==5) {
- extpart[0]=getParticleData(ParticleID::pi0);
- extpart[1]=getParticleData(ParticleID::pi0);
- extpart[2]=getParticleData(ParticleID::Kminus);
- }
- else if(imode==6) {
- extpart[0]=getParticleData(ParticleID::Kminus);
- extpart[1]=getParticleData(ParticleID::piminus);
- extpart[2]=getParticleData(ParticleID::piplus);
- }
- else if(imode==7) {
- extpart[0]=getParticleData(ParticleID::piminus);
- extpart[1]=getParticleData(ParticleID::Kbar0);
- extpart[2]=getParticleData(ParticleID::pi0);
- }
- else if(imode==8) {
- extpart[0]=getParticleData(ParticleID::piminus);
- extpart[1]=getParticleData(ParticleID::pi0);
- extpart[2]=getParticleData(ParticleID::eta);
- }
- else if(imode==9) {
- extpart[0]=getParticleData(ParticleID::K_S0);
- extpart[1]=getParticleData(ParticleID::piminus);
- extpart[2]=getParticleData(ParticleID::K_S0);
- }
- else if(imode==10) {
- extpart[0]=getParticleData(ParticleID::K_L0);
- extpart[1]=getParticleData(ParticleID::piminus);
- extpart[2]=getParticleData(ParticleID::K_L0);
- }
- else if(imode==11) {
- extpart[0]=getParticleData(ParticleID::K_S0);
- extpart[1]=getParticleData(ParticleID::piminus);
- extpart[2]=getParticleData(ParticleID::K_L0);
- }
- // conjugate the particles if needed
- if(icharge==3) {
- for(unsigned int ix=0;ix<3;++ix) {
- if(extpart[ix]->CC()) extpart[ix]=extpart[ix]->CC();
+ output[1]=getParticleData(ParticleID::pi0);
+ if(icharge==-3) {
+ for(unsigned int ix=0;ix<output.size();++ix) {
+ if(output[ix]->CC()) output[ix]=output[ix]->CC();
+ }
}
}
- // return the answer
- return extpart;
+ else {
+ output[1]=getParticleData(ParticleID::piminus);
+ }
+ return output;
}
-
diff --git a/Decay/WeakCurrents/EtaPiPiDefaultCurrent.h b/Decay/WeakCurrents/EtaPiPiDefaultCurrent.h
--- a/Decay/WeakCurrents/EtaPiPiDefaultCurrent.h
+++ b/Decay/WeakCurrents/EtaPiPiDefaultCurrent.h
@@ -1,632 +1,282 @@
// -*- C++ -*-
//
// EtaPiPiDefaultCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_EtaPiPiDefaultCurrent_H
#define HERWIG_EtaPiPiDefaultCurrent_H
//
// This is the declaration of the EtaPiPiDefaultCurrent class.
//
#include "WeakCurrent.h"
#include "Herwig/Utilities/Interpolator.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "Herwig/Decay/ResonanceHelpers.h"
+#include <numeric>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Decay
*
- * The EtaPiPiDefaultCurrent class implements the currents from Z.Phys.C58:445 (1992),
- * this paper uses the form from Z.Phys.C48:445 (1990) for the \f$a_1\f$ width and
- * is the default model in TAUOLA.
+ * The EtaPiPiDefaultCurrent class implements the current from Z.Phys.C58:445 (1992),
+ * for \f$ \pi^- \pi^0 \eta \f$.
*
- * The following three meson modes are implemented.
*
- * - \f$ \pi^- \pi^- \pi^+ \f$, (imode=0)
- * - \f$ \pi^0 \pi^0 \pi^- \f$, (imode=1)
- * - \f$ K^- \pi^- K^+ \f$, (imode=2)
- * - \f$ K^0 \pi^- \bar{K}^0\f$, (imode=3)
- * - \f$ K^- \pi^0 K^0 \f$, (imode=4)
- * - \f$ \pi^0 \pi^0 K^- \f$, (imode=5)
- * - \f$ K^- \pi^- \pi^+ \f$, (imode=6)
- * - \f$ \pi^- \bar{K}^0 \pi^0 \f$, (imode=7)
- * - \f$ \pi^- \pi^0 \eta \f$, (imode=8)
- *
- * using the currents from TAUOLA
- *
- *
- * @see WeakCurrent,
- * @see WeakCurrent.
- * @see Defaulta1MatrixElement
+ * @see WeakCurrent
*
*/
class EtaPiPiDefaultCurrent: public WeakCurrent {
- /**
- * The matrix element for the running \f$a_1\f$ width is a friend to
- * keep some members private.
- */
- friend class Defaulta1MatrixElement;
-
public:
/**
* Default constructor
*/
EtaPiPiDefaultCurrent();
/**
* Hadronic current. This method is purely virtual and must be implemented in
* all classes inheriting from this one.
* @param resonance If specified only include terms with this particle
* @param Itotal If specified the total isospin of the current
* @param I3 If specified the thrid component of isospin
* @param imode The mode
* @param ichan The phase-space channel the current is needed for.
* @param scale The invariant mass of the particles in the current.
* @param outgoing The particles produced in the decay
* @param momenta The momenta of the particles produced in the decay
* @param meopt Option for the calculation of the matrix element
* @return The current.
*/
virtual vector<LorentzPolarizationVectorE>
current(tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
const int imode, const int ichan,Energy & scale,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption meopt) const;
/**
* Accept the decay. Checks the mesons against the list.
* @param id The id's of the particles in the current.
* @return Can this current have the external particles specified.
*/
virtual bool accept(vector<int> id);
/**
* Return the decay mode number for a given set of particles in the current.
* Checks the mesons against the list.
* @param id The id's of the particles in the current.
* @return The number of the mode
*/
virtual unsigned int decayMode(vector<int> id);
/**
* The particles produced by the current. This returns the mesons for the mode.
* @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);
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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
public:
/** @name Methods for the construction of the phase space integrator. */
//@{
/**
* Complete the construction of the decay mode for integration.classes inheriting
* from this one.
* This method is purely virtual and must be implemented in the classes inheriting
* from WeakCurrent.
* @param icharge The total charge of the outgoing particles in the current.
* @param resonance If specified only include terms with this particle
* @param Itotal If specified the total isospin of the current
* @param I3 If specified the thrid component of isospin
* @param imode The mode in the current being asked for.
* @param mode The phase space mode for the integration
* @param iloc The location of the of the first particle from the current in
* the list of outgoing particles.
* @param ires The location of the first intermediate for the current.
* @param phase The prototype phase space channel for the integration.
* @param upp The maximum possible mass the particles in the current are
* allowed to have.
* @return Whether the current was sucessfully constructed.
*/
virtual bool createMode(int icharge, tcPDPtr resonance,
IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
unsigned int imode,PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp );
//@}
/**
* 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;
-
- /**
- * the matrix element for the \f$a_1\f$ decay to calculate the running width
- * @param imode The mode for which the matrix element is needed.
- * @param q2 The mass of the decaying off-shell \f$a_1\f$, \f$q^2\f$.
- * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
- * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
- * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
- * @param m1 The mass of the first outgoing particle.
- * @param m2 The mass of the second outgoing particle.
- * @param m3 The mass of the third outgoing particle.
- * @return The matrix element squared summed over spins.
- */
- double threeBodyMatrixElement(const int imode, const Energy2 q2,
- const Energy2 s3, const Energy2 s2,
- const Energy2 s1, const Energy m1,
- const Energy m2, const Energy m3) const;
-
-
-protected:
-
- /**
- * Helper class for form factors
- */
- struct FormFactors {
-
- /**
- * @param F1 The \f$F_1\f$ form factor
- */
- complex<InvEnergy> F1;
-
- /**
- * @param F2 The \f$F_2\f$ form factor
- */
- complex<InvEnergy> F2;
-
- /**
- * @param F3 The \f$F_3\f$ form factor
- */
- complex<InvEnergy> F3;
-
- /**
- * @param F4 The \f$F_4\f$ form factor
- */
- complex<InvEnergy> F4;
-
- /**
- * @param F5 The \f$F_5\f$ form factor
- */
- complex<InvEnergy3> F5;
-
- /**
- * Constructor
- * @param f1 The \f$F_1\f$ form factor
- * @param f2 The \f$F_2\f$ form factor
- * @param f3 The \f$F_3\f$ form factor
- * @param f4 The \f$F_4\f$ form factor
- * @param f5 The \f$F_5\f$ form factor
- */
- FormFactors(complex<InvEnergy> f1 = InvEnergy(),
- complex<InvEnergy> f2 = InvEnergy(),
- complex<InvEnergy> f3 = InvEnergy(),
- complex<InvEnergy> f4 = InvEnergy(),
- complex<InvEnergy3> f5 = InvEnergy3())
- : F1(f1), F2(f2), F3(f3), F4(f4), F5(f5) {}
- };
-
-protected:
-
- /**
- * Can the current handle a particular set of mesons.
- * As this current includes all the allowed modes this is always true.
- */
- virtual bool acceptMode(int) const;
-
- /**
- * Calculate the form factor for the current. Implements the form factors
- * described above.
- * @param ichan The phase space channel
- * @param imode The mode
- * @param q2 The scale \f$q^2\f$ for the current.
- * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
- * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
- * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
- */
- virtual FormFactors calculateFormFactors(const int ichan, const int imode,
- Energy2 q2,
- Energy2 s1, Energy2 s2, Energy2 s3) const;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** 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 {return new_ptr(*this);}
//@}
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();
- /**
- * Initialize this object to the begining of the run phase.
- */
- virtual void doinitrun();
-
- /**
- * Check sanity of the object during the setup phase.
- */
- virtual void doupdate();
- //@}
-
private:
/**
* Private and non-existent assignment operator.
*/
EtaPiPiDefaultCurrent & operator=(const EtaPiPiDefaultCurrent &);
private:
/**
* The \f$\rho\f$ Breit-Wigner for the \f$F_{1,2,3}\f$ form factors.
* @param q2 The scale \f$q^2\f$ for the Breit-Wigner
* @param ires Which \f$\rho\f$ multiplet
* @return The Breit-Wigner
*/
Complex BrhoF123(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
- norm+=_rhoF123wgts[ix];
- }
+ Complex output(0.);
+ Complex norm = std::accumulate(_rhoF123wgts.begin(),_rhoF123wgts.end(),Complex(0.));
if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
- output+=_rhoF123wgts[ix]*rhoKBreitWigner(q2,0,ix);
+ for(unsigned int ix=0;ix<_rhoF123wgts.size();++ix) {
+ output+=_rhoF123wgts[ix]*
+ Resonance::BreitWignerPWave(q2,_rhoF123masses[ix],
+ _rhoF123widths[ix],_mpi,_mpi);
}
}
else {
- unsigned int temp(ires);
- if(temp<_rhoF123wgts.size()&&temp<3)
- output=_rhoF123wgts[temp]*rhoKBreitWigner(q2,0,temp);
- else
- output=0.;
+ assert(ires<=int(_rhoF123wgts.size()));
+ output=_rhoF123wgts[ires]*
+ Resonance::BreitWignerPWave(q2,_rhoF123masses[ires],
+ _rhoF123widths[ires],_mpi,_mpi);
}
return output/norm;
}
/**
* The \f$\rho\f$ Breit-Wigner for the \f$F_5\f$ form factors.
* @param q2 The scale \f$q^2\f$ for the Breit-Wigner
* @param ires Which \f$\rho\f$ multiplet
* @return The Breit-Wigner
*/
Complex BrhoF5(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
- norm+=_rhoF5wgts[ix];
- }
+ Complex output(0.);
+ Complex norm = std::accumulate(_rhoF5wgts.begin(),_rhoF5wgts.end(),Complex(0.0));
if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
- output+=_rhoF5wgts[ix]*rhoKBreitWigner(q2,1,ix);
+ for(unsigned int ix=0;ix<_rhoF5wgts.size();++ix) {
+ output+=_rhoF5wgts[ix]*
+ Resonance::BreitWignerPWave(q2,_rhoF5masses[ix],
+ _rhoF5widths[ix],_mpi,_mpi);
}
}
else {
- unsigned int temp(ires);
- if(temp<_rhoF5wgts.size()&&temp<3) {
- output=_rhoF5wgts[temp]*rhoKBreitWigner(q2,1,temp);
- }
+ assert(ires<=int(_rhoF123wgts.size()));
+ output=_rhoF5wgts[ires]*
+ Resonance::BreitWignerPWave(q2,_rhoF5masses[ires],
+ _rhoF5widths[ires],_mpi,_mpi);
}
return output/norm;
}
- /**
- * The \f$K^*\f$ Breit-Wigner for the \f$F_{1,2,3}\f$ form factors.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @param ires Which \f$\rho\f$ multiplet
- * @return The Breit-Wigner
- */
- Complex BKstarF123(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
- norm+=_kstarF123wgts[ix];
- }
- if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
- output+=_kstarF123wgts[ix]*rhoKBreitWigner(q2,2,ix);
- }
- }
- else {
- unsigned int temp(ires);
- if(temp<_kstarF123wgts.size()&&temp<3) {
- output=_kstarF123wgts[temp]*rhoKBreitWigner(q2,2,temp);
- }
- }
- return output/norm;
- }
-
- /**
- * The \f$K^*\f$ Breit-Wigner for the \f$F_5\f$ form factors.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @param ires Which \f$\rho\f$ multiplet
- * @return The Breit-Wigner
- */
- Complex BKstarF5(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
- norm+=_kstarF5wgts[ix];
- }
- if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
- output+=_kstarF5wgts[ix]*rhoKBreitWigner(q2,3,ix);
- }
- }
- else {
- unsigned int temp(ires);
- if(temp<_kstarF5wgts.size()&&temp<3) {
- output=_kstarF5wgts[ires]*rhoKBreitWigner(q2,3,temp);
- }
- }
- return output/norm;
- }
-
- /**
- * Mixed Breit Wigner for the \f$F_5\f$ form factor
- * @param si The scale \f$s_1\f$.
- * @param sj The scale \f$s_2\f$.
- * @param ires Which resonances to use
- * @return The mixed Breit-Wigner
- */
- Complex FKrho(Energy2 si,Energy2 sj,int ires) const {
- Complex output;
- if(ires<0){output = _rhoKstarwgt*BKstarF123(si,-1)+BrhoF123(sj,-1);}
- else if(ires%2==0){output= _rhoKstarwgt*BKstarF123(si,ires/2);}
- else if(ires%2==1){output=BrhoF123(sj,ires/2);}
- output /=(1.+_rhoKstarwgt);
- return output;
- }
-
- /**
- * \f$a_1\f$ Breit-Wigner
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @return The Breit-Wigner
- */
- Complex a1BreitWigner(Energy2 q2) const {
- if(!_a1opt)
- return Resonance::BreitWignera1(q2,_a1mass,_a1width);
- Complex ii(0.,1.);
- Energy2 m2(_a1mass*_a1mass);
- Energy q(sqrt(q2));
- Energy width = (*_a1runinter)(q2);
- return m2/(m2-q2-ii*q*width);
- }
-
- /**
- * The \f$K_1\f$ Breit-Wigner
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @return The Breit-Wigner
- */
- Complex K1BreitWigner(Energy2 q2) const {
- Energy2 m2 = sqr(_k1mass);
- Complex ii(0.,1.);
- complex<Energy2> fact(m2 - ii*_k1mass*_k1width);
- return fact/(fact-q2);
- }
-
- /**
- * Initialize the \f$a_1\f$ running width
- * @param iopt Initialization option (-1 full calculation, 0 set up the interpolation)
- */
- void inita1Width(int iopt);
-
- /**
- * Breit-Wigners for the \f$\rho\f$ and \f$K^*\f$.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner.
- * @param itype The type of Breit-Wigner, \e i.e. which masses and widths to use.x
- * @param ires Which multiplet to use.
- */
- Complex rhoKBreitWigner(Energy2 q2,unsigned int itype,unsigned int ires) const;
-
private:
/**
* Parameters for the \f$\rho\f$ Breit-Wigner in the
* \f$F_{1,2,3}\f$ form factors.
*/
vector<double> _rhoF123wgts;
-
- /**
- * Parameters for the \f$K^*\f$ Breit-Wigner in the
- * \f$F_{1,2,3}\f$ form factors.
- */
- vector<double> _kstarF123wgts;
/**
* Parameters for the \f$\rho\f$ Breit-Wigner in the
* \f$F_5\f$ form factors.
*/
vector<double> _rhoF5wgts;
/**
- * Parameters for the \f$K^*\f$ Breit-Wigner in the
- * \f$F_5\f$ form factors.
- */
- vector<double> _kstarF5wgts;
-
- /**
- * The relative weight of the \f$\rho\f$ and \f$K^*\f$ where needed.
- */
- double _rhoKstarwgt;
-
- /**
- * The \f$a_1\f$ width for the running \f$a_1\f$ width calculation.
- */
- vector<Energy> _a1runwidth;
-
- /**
- * The \f$q^2\f$ for the running \f$a_1\f$ width calculation.
- */
- vector<Energy2> _a1runq2;
-
-
- /**
- * The interpolator for the running \f$a_1\f$ width calculation.
- */
- Interpolator<Energy,Energy2>::Ptr _a1runinter;
-
- /**
- * Initialize the running \f$a_1\f$ width.
- */
- bool _initializea1;
-
- /**
- * The mass of the \f$a_1\f$ resonances.
- */
- Energy _a1mass;
-
- /**
- * The width of the \f$a_1\f$ resonances.
- */
- Energy _a1width;
-
- /**
- * The mass of the \f$aK1\f$ resonances.
- */
- Energy _k1mass;
-
- /**
- * The width of the \f$K_1\f$ resonances.
- */
- Energy _k1width;
-
- /**
* The pion decay constant, \f$f_\pi\f$.
*/
Energy _fpi;
/**
* The pion mass
*/
Energy _mpi;
/**
- * The kaon mass
- */
- Energy _mK;
-
- /**
- * use local values of the \f$\rho\f$ masses and widths
- */
- bool _rhoparameters;
-
- /**
* The \f$\rho\f$ masses for the \f$F_{1,2,3}\f$ form factors.
*/
vector<Energy> _rhoF123masses;
/**
* The \f$\rho\f$ masses for the \f$F_5\f$ form factors.
*/
vector<Energy> _rhoF5masses;
/**
* The \f$\rho\f$ widths for the \f$F_{1,2,3}\f$ form factors.
*/
vector<Energy> _rhoF123widths;
/**
* The \f$\rho\f$ widths for the \f$F_5\f$ form factors.
*/
vector<Energy> _rhoF5widths;
-
- /**
- * use local values of the \f$K^*\f$ resonances masses and widths
- */
- bool _kstarparameters;
-
- /**
- * The \f$K^*\f$ masses for the \f$F_{1,2,3}\f$ form factors.
- */
- vector<Energy> _kstarF123masses;
-
- /**
- * The \f$K^*\f$ masses for the \f$F_5\f$ form factors.
- */
- vector<Energy> _kstarF5masses;
-
- /**
- * The \f$K^*\f$ widths for the \f$F_{1,2,3}\f$ form factors.
- */
- vector<Energy> _kstarF123widths;
-
- /**
- * The \f$K^*\f$ widths for the \f$F_5\f$ form factors.
- */
- vector<Energy> _kstarF5widths;
-
- /**
- * Use local values of the \f$a_1\f$ parameters
- */
- bool _a1parameters;
-
- /**
- * Use local values of the \f$K_1\f$ parameters
- */
- bool _k1parameters;
-
- /**
- * Option for the \f$a_1\f$ width
- */
- bool _a1opt;
-
- /**
- * The maximum mass of the hadronic system
- */
- Energy _maxmass;
-
- /**
- * The maximum mass when the running width was calculated
- */
- Energy _maxcalc;
-
};
}
#endif /* HERWIG_EtaPiPiDefaultCurrent_H */
diff --git a/Decay/WeakCurrents/Makefile.am b/Decay/WeakCurrents/Makefile.am
--- a/Decay/WeakCurrents/Makefile.am
+++ b/Decay/WeakCurrents/Makefile.am
@@ -1,56 +1,54 @@
BUILT_SOURCES = WeakCurrents__all.cc
CLEANFILES = WeakCurrents__all.cc
WeakCurrents__all.cc : $(DIR_H_FILES) $(DIR_CC_FILES) Makefile
@echo "Concatenating .cc files into $@"
@$(top_srcdir)/cat_with_cpplines $(DIR_CC_FILES) > $@
EXTRA_DIST = $(ALL_H_FILES) $(ALL_CC_FILES)
DIR_H_FILES = $(addprefix $(srcdir)/,$(ALL_H_FILES))
ALL_H_FILES = \
FourPionNovosibirskCurrent.h \
ScalarMesonCurrent.h\
ThreeMesonCurrentBase.h \
-ThreeMesonDefaultCurrent.h\
ThreePionDefaultCurrent.h\
OneKaonTwoPionDefaultCurrent.h \
TwoKaonOnePionDefaultCurrent.h \
EtaPiPiDefaultCurrent.h \
ThreePionCLEOCurrent.h\
TwoPionRhoCurrent.h\
KPiKStarCurrent.h\
TwoPionPhotonCurrent.h\
VectorMesonCurrent.h\
FivePionCurrent.h \
KPiCurrent.h\
KaonThreeMesonCurrent.h\
TwoPionCzyzCurrent.h\
TwoKaonCzyzCurrent.h\
ThreePionCzyzCurrent.h\
FourPionCzyzCurrent.h\
EtaPiPiCurrent.h
DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES))
ALL_CC_FILES = \
FourPionNovosibirskCurrent.cc \
ScalarMesonCurrent.cc \
ThreeMesonCurrentBase.cc \
-ThreeMesonDefaultCurrent.cc \
ThreePionDefaultCurrent.cc \
OneKaonTwoPionDefaultCurrent.cc \
TwoKaonOnePionDefaultCurrent.cc \
EtaPiPiDefaultCurrent.cc \
ThreePionCLEOCurrent.cc \
TwoPionRhoCurrent.cc \
KPiKStarCurrent.cc \
TwoPionPhotonCurrent.cc \
VectorMesonCurrent.cc \
FivePionCurrent.cc \
KPiCurrent.cc \
KaonThreeMesonCurrent.cc \
TwoPionCzyzCurrent.cc \
TwoKaonCzyzCurrent.cc \
ThreePionCzyzCurrent.cc \
FourPionCzyzCurrent.cc\
EtaPiPiCurrent.cc
diff --git a/Decay/WeakCurrents/ThreeMesonDefaultCurrent.cc b/Decay/WeakCurrents/ThreeMesonDefaultCurrent.cc
deleted file mode 100644
--- a/Decay/WeakCurrents/ThreeMesonDefaultCurrent.cc
+++ /dev/null
@@ -1,1053 +0,0 @@
-// -*- C++ -*-
-//
-// ThreeMesonDefaultCurrent.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
-// Copyright (C) 2002-2017 The Herwig Collaboration
-//
-// Herwig is licenced under version 3 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
-//
-// This is the implementation of the non-inlined, non-templated member
-// functions of the ThreeMesonDefaultCurrent class.
-//
-
-#include "ThreeMesonDefaultCurrent.h"
-#include "ThePEG/Interface/ClassDocumentation.h"
-#include "ThePEG/Interface/Switch.h"
-#include "ThePEG/Interface/Parameter.h"
-#include "ThePEG/Interface/ParVector.h"
-#include "ThePEG/Persistency/PersistentOStream.h"
-#include "ThePEG/Persistency/PersistentIStream.h"
-#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
-#include "ThePEG/Utilities/DescribeClass.h"
-
-using namespace Herwig;
-using namespace ThePEG;
-
-DescribeClass<ThreeMesonDefaultCurrent,ThreeMesonCurrentBase>
-describeHerwigThreeMesonDefaultCurrent("Herwig::ThreeMesonDefaultCurrent",
- "HwWeakCurrents.so");
-HERWIG_INTERPOLATOR_CLASSDESC(ThreeMesonDefaultCurrent,Energy,Energy2)
-
-
-ThreeMesonDefaultCurrent::ThreeMesonDefaultCurrent() {
- // the pion decay constant
- _fpi=130.7*MeV/sqrt(2.);
- _mpi=ZERO;_mK=ZERO;
- // set the initial weights for the resonances
- // the rho weights
- _rhoF123wgts.push_back(1.0);_rhoF123wgts.push_back(-0.145);
- _rhoF123wgts.push_back(0.);
- _rhoF5wgts.push_back(-26.);_rhoF5wgts.push_back(6.5);
- _rhoF5wgts.push_back(1.);
- // the Kstar weights
- _kstarF123wgts.push_back(1.);
- _kstarF5wgts.push_back(1.);
- // relative rho/Kstar weights
- _rhoKstarwgt=-0.2;
- // local values of the a_1 parameters
- _a1parameters=true;_a1mass=1.251*GeV;_a1width=0.599*GeV;
- _a1opt=true;
- // local values of the K_1 parameters
- _k1parameters=true;_k1mass=1.402*GeV,_k1width=0.174*GeV;
- // local values of the rho parameters
- _rhoparameters=true;
- _rhoF123masses.push_back(0.773*GeV);_rhoF123masses.push_back(1.370*GeV);
- _rhoF123masses.push_back(1.750*GeV);
- _rhoF123widths.push_back(0.145*GeV);_rhoF123widths.push_back(0.510*GeV);
- _rhoF123widths.push_back(0.120*GeV);
- _rhoF5masses.push_back(0.773*GeV);_rhoF5masses.push_back(1.500*GeV);
- _rhoF5masses.push_back(1.750*GeV);
- _rhoF5widths.push_back(0.145*GeV);_rhoF5widths.push_back(0.220*GeV);
- _rhoF5widths.push_back(0.120*GeV);
- // local values for the Kstar parameters
- _kstarparameters=true;
- _kstarF123masses.push_back(0.8921*GeV);
- _kstarF123widths.push_back(0.0513*GeV);
- _kstarF5masses.push_back(0.8921*GeV);
- _kstarF5widths.push_back(0.0513*GeV);
- // initialization of the a_1 running width
- _initializea1=false;
- double a1q2in[200]={0,15788.6,31577.3,47365.9,63154.6,78943.2,94731.9,110521,
- 126309,142098,157886,173675,189464,205252,221041,236830,
- 252618,268407,284196,299984,315773,331562,347350,363139,
- 378927,394716,410505,426293,442082,457871,473659,489448,
- 505237,521025,536814,552603,568391,584180,599969,615757,
- 631546,647334,663123,678912,694700,710489,726278,742066,
- 757855,773644,789432,805221,821010,836798,852587,868375,
- 884164,899953,915741,931530,947319,963107,978896,994685,
- 1.01047e+06,1.02626e+06,1.04205e+06,1.05784e+06,1.07363e+06,
- 1.08942e+06,1.10521e+06,1.12099e+06,1.13678e+06,1.15257e+06,
- 1.16836e+06,1.18415e+06,1.19994e+06,1.21573e+06,1.23151e+06,
- 1.2473e+06,1.26309e+06,1.27888e+06,1.29467e+06,1.31046e+06,
- 1.32625e+06,1.34203e+06,1.35782e+06,1.37361e+06,1.3894e+06,
- 1.40519e+06,1.42098e+06,1.43677e+06,1.45256e+06,1.46834e+06
- ,1.48413e+06,1.49992e+06,1.51571e+06,1.5315e+06,1.54729e+06,
- 1.56308e+06,1.57886e+06,1.59465e+06,1.61044e+06,1.62623e+06,
- 1.64202e+06,1.65781e+06,1.6736e+06,1.68939e+06,1.70517e+06,
- 1.72096e+06,1.73675e+06,1.75254e+06,1.76833e+06,1.78412e+06,
- 1.79991e+06,1.81569e+06,1.83148e+06,1.84727e+06,1.86306e+06,
- 1.87885e+06,1.89464e+06,1.91043e+06,1.92621e+06,1.942e+06,
- 1.95779e+06,1.97358e+06,1.98937e+06,2.00516e+06,2.02095e+06,
- 2.03674e+06,2.05252e+06,2.06831e+06,2.0841e+06,2.09989e+06,
- 2.11568e+06,2.13147e+06,2.14726e+06,2.16304e+06,2.17883e+06,
- 2.19462e+06,2.21041e+06,2.2262e+06,2.24199e+06,2.25778e+06,
- 2.27356e+06,2.28935e+06,2.30514e+06,2.32093e+06,2.33672e+06,
- 2.35251e+06,2.3683e+06,2.38409e+06,2.39987e+06,2.41566e+06,
- 2.43145e+06,2.44724e+06,2.46303e+06,2.47882e+06,2.49461e+06,
- 2.51039e+06,2.52618e+06,2.54197e+06,2.55776e+06,2.57355e+06,
- 2.58934e+06,2.60513e+06,2.62092e+06,2.6367e+06,2.65249e+06,
- 2.66828e+06,2.68407e+06,2.69986e+06,2.71565e+06,2.73144e+06,
- 2.74722e+06,2.76301e+06,2.7788e+06,2.79459e+06,2.81038e+06,
- 2.82617e+06,2.84196e+06,2.85774e+06,2.87353e+06,2.88932e+06,
- 2.90511e+06,2.9209e+06,2.93669e+06,2.95248e+06,2.96827e+06,
- 2.98405e+06,2.99984e+06,3.01563e+06,3.03142e+06,3.04721e+06,
- 3.063e+06,3.07879e+06,3.09457e+06,3.11036e+06,3.12615e+06,
- 3.14194e+06};
- double a1widthin[200]={0,0,0,0,0,0,0,0,0,0,0,0,0.00153933,0.0136382,0.0457614,
- 0.105567,0.199612,0.333825,0.513831,0.745192,1.0336,1.38501,
- 1.80581,2.30295,2.88403,3.5575,4.33278,5.22045,6.23243,
- 7.38223,8.68521,10.1589,11.8234,13.7018,15.8206,18.2107,
- 20.9078,23.9533,27.3954,31.2905,35.7038,40.7106,46.3984,
- 52.8654,60.2207,68.581,78.0637,88.7754,100.794,114.145,
- 128.783,144.574,161.299,178.683,196.426,214.248,231.908,
- 249.221,266.059,282.336,298.006,313.048,327.46,341.254,
- 354.448,367.066,379.133,390.677,401.726,412.304,422.439,
- 432.155,441.474,450.419,459.01,467.267,475.207,482.847,
- 490.203,497.29,504.121,510.71,517.068,523.207,529.138,
- 534.869,540.411,545.776,550.961,556.663,560.851,565.566,
- 570.137,574.569,578.869,583.041,587.091,591.023,594.843,
- 598.553,602.16,605.664,609.072,612.396,615.626,618.754,
- 621.796,624.766,627.656,630.47,633.21,635.878,638.5,
- 641.006,643.471,645.873,648.213,650.493,652.715,654.88,
- 656.99,659.047,661.052,663.007,664.963,666.771,668.6,
- 670.351,672.075,673.828,675.397,676.996,678.567,680.083,
- 681.589,683.023,684.457,685.825,687.18,688.499,689.789,
- 691.058,692.284,693.501,694.667,695.82,696.947,698.05,
- 699.129,700.186,701.221,702.234,703.226,704.198,705.158,
- 706.085,707.001,707.899,708.78,709.644,710.474,711.334,
- 712.145,712.943,713.727,714.505,715.266,716.015,716.751,
- 717.474,718.183,718.88,719.645,720.243,720.91,721.565,
- 722.211,722.851,723.473,724.094,724.697,725.296,725.886,
- 726.468,727.041,727.608,728.166,728.718,729.262,729.808,
- 730.337,730.856,731.374,731.883,732.386,732.884,733.373,
- 733.859,734.339,734.813};
-
- vector<double> tmp1(a1widthin,a1widthin+200);
- _a1runwidth.clear();
- std::transform(tmp1.begin(), tmp1.end(),
- back_inserter(_a1runwidth),
- [](double x){return x*GeV;});
-
- vector<double> tmp2(a1q2in,a1q2in+200);
- _a1runq2.clear();
- std::transform(tmp2.begin(), tmp2.end(),
- back_inserter(_a1runq2),
- [](double x){return x*GeV2;});
-
- _maxmass=ZERO;
- _maxcalc=ZERO;
-}
-
-void ThreeMesonDefaultCurrent::doinit() {
- ThreeMesonCurrentBase::doinit();
- // the particles we will use a lot
- tPDPtr a1(getParticleData(ParticleID::a_1minus)),
- k1(getParticleData(ParticleID::Kstar_1minus)),pi0(getParticleData(ParticleID::pi0)),
- piplus(getParticleData(ParticleID::piplus)),
- piminus(getParticleData(ParticleID::piminus));
- // masses for the running widths
- _mpi=piplus->mass();
- _mK=getParticleData(ParticleID::Kminus)->mass();
- // the charged rho resonances
- tPDPtr rhoc[3]={getParticleData(-213),getParticleData(-100213),
- getParticleData(-30213)};
- // the charged K* resonances
- tPDPtr Kstarc[3]={getParticleData(-323),getParticleData(-100323),
- getParticleData(-30323)};
- if(!_a1parameters) {
- _a1mass=a1->mass();
- _a1width=a1->width();
- }
- // mass and width of the k_1
- if(!_k1parameters) {
- _k1mass=k1->mass();
- _k1width=k1->width();
- }
- // initialise the a_1 running width calculation
- inita1Width(-1);
- // rho parameters in the base classs
- tcPDPtr temp;
- unsigned int ix;
- if(_rhoparameters&&_rhoF123masses.size()<3) {
- ix = _rhoF123masses.size();
- _rhoF123masses.resize(3);_rhoF123widths.resize(3);
- for(;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF123masses[ix]=rhoc[ix]->mass();
- _rhoF123widths[ix]=rhoc[ix]->width();
- }
- }
- }
- else if(!_rhoparameters) {
- _rhoF123masses.resize(3);_rhoF123widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF123masses[ix]=rhoc[ix]->mass();
- _rhoF123widths[ix]=rhoc[ix]->width();
- }
- }
- }
- // K star parameters in the base class
- if(_kstarparameters&&_kstarF123masses.size()<3) {
- ix = _kstarF123masses.size();
- _kstarF123masses.resize(3);_kstarF123widths.resize(3);
- for(;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF123masses[ix]=Kstarc[ix]->mass();
- _kstarF123widths[ix]=Kstarc[ix]->width();
- }
- }
- }
- else if(!_kstarparameters) {
- _kstarF123masses.resize(3);_kstarF123widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF123masses[ix]=Kstarc[ix]->mass();
- _kstarF123widths[ix]=Kstarc[ix]->width();
- }
- }
- }
- // rho parameters here
- if(_rhoparameters&&_rhoF5masses.size()<3) {
- ix = _rhoF5masses.size();
- _rhoF5masses.resize(3);_rhoF5widths.resize(3);
- for(;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF5masses[ix]=rhoc[ix]->mass();
- _rhoF5widths[ix]=rhoc[ix]->width();
- }
- }
- }
- else if(!_rhoparameters) {
- _rhoF5masses.resize(3);_rhoF5widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(rhoc[ix]) {
- _rhoF5masses[ix]=rhoc[ix]->mass();
- _rhoF5widths[ix]=rhoc[ix]->width();
- }
- }
- }
- // Kstar parameters here
- if(_kstarparameters&&_kstarF5widths.size()<3) {
- ix = _kstarF5masses.size();
- _kstarF5masses.resize(3);_kstarF5widths.resize(3);
- for(;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF5masses[ix]=Kstarc[ix]->mass();
- _kstarF5widths[ix]=Kstarc[ix]->width();
- }
- }
- }
- else if(!_kstarparameters) {
- _kstarF5masses.resize(3);_kstarF5widths.resize(3);
- for(ix=0;ix<3;++ix) {
- if(Kstarc[ix]) {
- _kstarF5masses[ix]=Kstarc[ix]->mass();
- _kstarF5widths[ix]=Kstarc[ix]->width();
- }
- }
- }
-}
-
-void ThreeMesonDefaultCurrent::persistentOutput(PersistentOStream & os) const {
- os << _rhoF123wgts << _kstarF123wgts << _rhoF5wgts << _kstarF5wgts
- << _rhoKstarwgt << ounit(_a1runwidth,GeV)<< ounit(_a1runq2,GeV2)
- << _initializea1
- << ounit(_a1mass,GeV)<< ounit(_a1width,GeV)<< ounit(_k1mass,GeV)
- << ounit(_k1width,GeV)<< ounit(_fpi,GeV) << ounit(_mpi,GeV)<< ounit(_mK,GeV)
- <<_rhoparameters << ounit(_rhoF123masses,GeV) << ounit(_rhoF5masses,GeV)
- << ounit(_rhoF123widths,GeV)
- << ounit(_rhoF5widths,GeV) << _kstarparameters << ounit(_kstarF123masses,GeV)
- <<ounit(_kstarF5masses,GeV)
- << ounit(_kstarF123widths,GeV) << ounit(_kstarF5widths,GeV) << _a1parameters
- << _k1parameters
- << _a1opt << ounit(_maxmass,GeV) << ounit(_maxcalc,GeV) << _a1runinter;
-}
-
-void ThreeMesonDefaultCurrent::persistentInput(PersistentIStream & is, int) {
- is >> _rhoF123wgts >> _kstarF123wgts >> _rhoF5wgts >> _kstarF5wgts
- >> _rhoKstarwgt >> iunit(_a1runwidth,GeV) >> iunit(_a1runq2,GeV2)
- >> _initializea1
- >> iunit(_a1mass,GeV) >> iunit(_a1width,GeV) >> iunit(_k1mass,GeV)
- >> iunit(_k1width,GeV) >> iunit(_fpi,GeV) >> iunit(_mpi,GeV) >> iunit(_mK,GeV)
- >>_rhoparameters >> iunit(_rhoF123masses,GeV) >> iunit(_rhoF5masses,GeV)
- >> iunit(_rhoF123widths,GeV)
- >> iunit(_rhoF5widths,GeV) >> _kstarparameters >> iunit(_kstarF123masses,GeV)
- >>iunit(_kstarF5masses,GeV)
- >> iunit(_kstarF123widths,GeV) >> iunit(_kstarF5widths,GeV) >> _a1parameters
- >> _k1parameters
- >> _a1opt >> iunit(_maxmass,GeV) >> iunit(_maxcalc,GeV) >> _a1runinter;
-}
-
-void ThreeMesonDefaultCurrent::Init() {
-
- static ClassDocumentation<ThreeMesonDefaultCurrent> documentation
- ("The ThreeMesonDefaultCurrent class is designed to implement "
- "the three meson decays of the tau, ie pi- pi- pi+, pi0 pi0 pi-, "
- "K- pi- K+, K0 pi- Kbar0, K- pi0 K0,pi0 pi0 K-, K- pi- pi+, "
- "pi- Kbar0 pi0, pi- pi0 eta. It uses the same currents as those in TAUOLA.",
- "The three meson decays of the tau, ie pi- pi- pi+, pi0 pi0 pi-, "
- "K- pi- K+, K0 pi- Kbar0, K- pi0 K0,pi0 pi0 K-, K- pi- pi+, "
- "and pi- Kbar0 pi0, pi- pi0 eta "
- "use the same currents as \\cite{Jadach:1993hs,Kuhn:1990ad,Decker:1992kj}.",
- "%\\cite{Jadach:1993hs}\n"
- "\\bibitem{Jadach:1993hs}\n"
- " S.~Jadach, Z.~Was, R.~Decker and J.~H.~Kuhn,\n"
- " %``The Tau Decay Library Tauola: Version 2.4,''\n"
- " Comput.\\ Phys.\\ Commun.\\ {\\bf 76}, 361 (1993).\n"
- " %%CITATION = CPHCB,76,361;%%\n"
- "%\\cite{Kuhn:1990ad}\n"
- "\\bibitem{Kuhn:1990ad}\n"
- " J.~H.~Kuhn and A.~Santamaria,\n"
- " %``Tau decays to pions,''\n"
- " Z.\\ Phys.\\ C {\\bf 48}, 445 (1990).\n"
- " %%CITATION = ZEPYA,C48,445;%%\n"
- "%\\cite{Decker:1992kj}\n"
- "\\bibitem{Decker:1992kj}\n"
- " R.~Decker, E.~Mirkes, R.~Sauer and Z.~Was,\n"
- " %``Tau decays into three pseudoscalar mesons,''\n"
- " Z.\\ Phys.\\ C {\\bf 58}, 445 (1993).\n"
- " %%CITATION = ZEPYA,C58,445;%%\n"
- );
-
-
- static ParVector<ThreeMesonDefaultCurrent,double> interfaceF123RhoWgt
- ("F123RhoWeight",
- "The weights of the different rho resonances in the F1,2,3 form factor",
- &ThreeMesonDefaultCurrent::_rhoF123wgts,
- 0, 0, 0, -1000, 1000, false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,double> interfaceF123KstarWgt
- ("F123KstarWeight",
- "The weights of the different Kstar resonances in the F1,2,3 form factor",
- &ThreeMesonDefaultCurrent::_kstarF123wgts,
- 0, 0, 0, -1000, 1000, false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,double> interfaceF5RhoWgt
- ("F5RhoWeight",
- "The weights of the different rho resonances in the F1,2,3 form factor",
- &ThreeMesonDefaultCurrent::_rhoF5wgts,
- 0, 0, 0, -1000, 1000, false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,double> interfaceF5KstarWgt
- ("F5KstarWeight",
- "The weights of the different Kstar resonances in the F1,2,3 form factor",
- &ThreeMesonDefaultCurrent::_kstarF5wgts,
- 0, 0, 0, -1000, 1000, false, false, true);
-
- static Parameter<ThreeMesonDefaultCurrent,double> interfaceRhoKstarWgt
- ("RhoKstarWgt",
- "The relative weights of the rho and K* in the F5 form factor",
- &ThreeMesonDefaultCurrent::_rhoKstarwgt, -0.2, -10., 10.,
- false, false, false);
-
- static Switch<ThreeMesonDefaultCurrent,bool> interfaceInitializea1
- ("Initializea1",
- "Initialise the calculation of the a_1 running width",
- &ThreeMesonDefaultCurrent::_initializea1, false, false, false);
- static SwitchOption interfaceInitializea1Initialization
- (interfaceInitializea1,
- "Yes",
- "Initialize the calculation",
- true);
- static SwitchOption interfaceInitializea1NoInitialization
- (interfaceInitializea1,
- "No",
- "Use the default values",
- false);
-
- static Switch<ThreeMesonDefaultCurrent,bool> interfaceRhoParameters
- ("RhoParameters",
- "Use local values of the rho meson masses and widths",
- &ThreeMesonDefaultCurrent::_rhoparameters, true, false, false);
- static SwitchOption interfaceRhoParameterstrue
- (interfaceRhoParameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfaceRhoParametersParticleData
- (interfaceRhoParameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<ThreeMesonDefaultCurrent,bool> interfaceKstarParameters
- ("KstarParameters",
- "Use local values of the rho meson masses and widths",
- &ThreeMesonDefaultCurrent::_kstarparameters, true, false, false);
- static SwitchOption interfaceKstarParameterstrue
- (interfaceKstarParameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfaceKstarParametersParticleData
- (interfaceKstarParameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<ThreeMesonDefaultCurrent,bool> interfacea1Parameters
- ("a1Parameters",
- "Use local values of the rho meson masses and widths",
- &ThreeMesonDefaultCurrent::_a1parameters, true, false, false);
- static SwitchOption interfacea1Parameterstrue
- (interfacea1Parameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfacea1ParametersParticleData
- (interfacea1Parameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<ThreeMesonDefaultCurrent,bool> interfaceK1Parameters
- ("K1Parameters",
- "Use local values of the rho meson masses and widths",
- &ThreeMesonDefaultCurrent::_k1parameters, true, false, false);
- static SwitchOption interfaceK1Parameterstrue
- (interfaceK1Parameters,
- "Local",
- "Use local values of the parameters",
- true);
- static SwitchOption interfaceK1ParametersParticleData
- (interfaceK1Parameters,
- "ParticleData",
- "Use the masses and wdiths from the particle data objects",
- false);
-
- static Switch<ThreeMesonDefaultCurrent,bool> interfacea1WidthOption
- ("a1WidthOption",
- "Option for the treatment of the a1 width",
- &ThreeMesonDefaultCurrent::_a1opt, true, false, false);
- static SwitchOption interfacea1WidthOptionLocal
- (interfacea1WidthOption,
- "Local",
- "Use a calculation of the running width based on the parameters as"
- " interpolation table.",
- true);
- static SwitchOption interfacea1WidthOptionParam
- (interfacea1WidthOption,
- "Kuhn",
- "Use the parameterization of Kuhn and Santamaria for default parameters."
- " This should only be used for testing vs TAUOLA",
- false);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfacea1RunningWidth
- ("a1RunningWidth",
- "The values of the a_1 width for interpolation to giving the running width.",
- &ThreeMesonDefaultCurrent::_a1runwidth, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy2> interfacea1RunningQ2
- ("a1RunningQ2",
- "The values of the q^2 for interpolation to giving the running width.",
- &ThreeMesonDefaultCurrent::_a1runq2, GeV2, -1, 1.0*GeV2, ZERO, 10.0*GeV2,
- false, false, true);
-
- static Parameter<ThreeMesonDefaultCurrent,Energy> interfaceA1Width
- ("A1Width",
- "The a_1 width if using local values.",
- &ThreeMesonDefaultCurrent::_a1width, GeV, 0.599*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
- static Parameter<ThreeMesonDefaultCurrent,Energy> interfaceA1Mass
- ("A1Mass",
- "The a_1 mass if using local values.",
- &ThreeMesonDefaultCurrent::_a1mass, GeV, 1.251*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
- static Parameter<ThreeMesonDefaultCurrent,Energy> interfaceK1Width
- ("K1Width",
- "The K_1 width if using local values.",
- &ThreeMesonDefaultCurrent::_k1width, GeV, 0.174*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
- static Parameter<ThreeMesonDefaultCurrent,Energy> interfaceK1Mass
- ("K1Mass",
- "The K_1 mass if using local values.",
- &ThreeMesonDefaultCurrent::_k1mass, GeV, 1.402*GeV, ZERO, 10.0*GeV,
- false, false, false);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfacerhoF123masses
- ("rhoF123masses",
- "The masses for the rho resonances if used local values",
- &ThreeMesonDefaultCurrent::_rhoF123masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfacerhoF123widths
- ("rhoF123widths",
- "The widths for the rho resonances if used local values",
- &ThreeMesonDefaultCurrent::_rhoF123widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfacerhoF5masses
- ("rhoF5masses",
- "The masses for the rho resonances if used local values",
- &ThreeMesonDefaultCurrent::_rhoF5masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfacerhoF5widths
- ("rhoF5widths",
- "The widths for the rho resonances if used local values",
- &ThreeMesonDefaultCurrent::_rhoF5widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfaceKstarF123masses
- ("KstarF123masses",
- "The masses for the Kstar resonances if used local values",
- &ThreeMesonDefaultCurrent::_kstarF123masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfaceKstarF123widths
- ("KstarF123widths",
- "The widths for the Kstar resonances if used local values",
- &ThreeMesonDefaultCurrent::_kstarF123widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfaceKstarF5masses
- ("KstarF5masses",
- "The masses for the Kstar resonances if used local values",
- &ThreeMesonDefaultCurrent::_kstarF5masses, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static ParVector<ThreeMesonDefaultCurrent,Energy> interfaceKstarF5widths
- ("KstarF5widths",
- "The widths for the Kstar resonances if used local values",
- &ThreeMesonDefaultCurrent::_kstarF5widths, GeV, -1, 1.0*GeV, ZERO, 10.0*GeV,
- false, false, true);
-
- static Parameter<ThreeMesonDefaultCurrent,Energy> interfaceFPi
- ("FPi",
- "The pion decay constant",
- &ThreeMesonDefaultCurrent::_fpi, MeV, 92.4*MeV, ZERO, 200.0*MeV,
- false, false, true);
-}
-
-// modes handled by this class
-bool ThreeMesonDefaultCurrent::acceptMode(int imode) const {
- return imode>=0&&imode<=8;
-}
-
-// calculate the form-factors
-ThreeMesonDefaultCurrent::FormFactors
-ThreeMesonDefaultCurrent::calculateFormFactors(const int ichan, const int imode,
- Energy2 q2, Energy2 s1,
- Energy2 s2, Energy2 s3) const {
- useMe();
- Complex F1, F2, F3, F4, F5;
- F1 = F2 = F3 = F4 = F5 = 0.0;
- // calculate the pi- pi- pi+ factor
- if(imode==0) {
- Complex a1fact(a1BreitWigner(q2)*2./3.);
- if(ichan<0) {
- F1= a1fact*BrhoF123(s1,-1);
- F2 =-a1fact*BrhoF123(s2,-1);
- }
- else if(ichan%2==0) F1 = a1fact*BrhoF123(s1, ichan/2);
- else if(ichan%2==1) F2 =-a1fact*BrhoF123(s2,(ichan-1)/2);
- }
- // calculate the pi0 pi0 pi- factor
- else if(imode==1) {
- Complex a1fact(a1BreitWigner(q2)*2./3.);
- if(ichan<0) {
- F1 = a1fact*BrhoF123(s1,-1);
- F2 =-a1fact*BrhoF123(s2,-1);
- }
- else if(ichan%2==0) F1 = a1fact*BrhoF123(s1, ichan/2);
- else if(ichan%2==1) F2 =-a1fact*BrhoF123(s2,(ichan-1)/2);
- }
- // calculate the K- pi - K+ factor
- else if(imode==2) {
- Complex a1fact(a1BreitWigner(q2)*sqrt(2.)/3.);
- if(ichan<0) {
- F1 =-a1fact*BKstarF123(s1,-1);
- F2 = a1fact*BrhoF123(s2,-1);
- F5 = BrhoF5(q2,-1)*FKrho(s1,s2,-1)*sqrt(2.);
- }
- else if(ichan%8==0) F1 =-a1fact*BKstarF123(s1,ichan/8);
- else if(ichan%8==1) F2 = a1fact*BrhoF123(s2,(ichan-1)/8);
- else if(ichan%8>=2) F5 = BrhoF5(q2,ichan/8)*FKrho(s1,s2,(ichan-2)%8)*sqrt(2.);
- }
- // calculate the K0 pi- K0bar
- else if(imode==3) {
- Complex a1fact(a1BreitWigner(q2)*sqrt(2.)/3.);
- if(ichan<0) {
- F1 =-a1fact*BKstarF123(s1,-1);
- F2 = a1fact*BrhoF123(s2,-1);
- F5 =-BrhoF5(q2,-1)*FKrho(s1,s2,-1)*sqrt(2.);
- }
- else if(ichan%8==0) F1 = -a1fact*BKstarF123(s1,ichan/8);
- else if(ichan%8==1) F2 = a1fact*BrhoF123(s2,(ichan-1)/8);
- else if(ichan%8>=2) F5 = -BrhoF5(q2,ichan/8)*FKrho(s1,s2,(ichan-2)%8)*sqrt(2.);
- }
- // calculate the K- pi0 k0
- else if(imode==4) {
- Complex a1fact(a1BreitWigner(q2));
- if(ichan<0){F2 =-a1fact*BrhoF123(s2,-1);}
- else{F2 =-a1fact*BrhoF123(s2,ichan);}
- }
- // calculate the pi0 pi0 K-
- else if(imode==5) {
- Complex K1fact(K1BreitWigner(q2)/6.);
- if(ichan<0) {
- F1 = K1fact*BKstarF123(s1,-1);
- F2 =-K1fact*BKstarF123(s2,-1);
- }
- else if(ichan%2==0) F1 = K1fact*BKstarF123(s1,ichan/2);
- else F2 =-K1fact*BKstarF123(s2,(ichan-1)/2);
- }
- // calculate the K- pi- pi+
- else if(imode==6) {
- Complex K1fact(K1BreitWigner(q2)*sqrt(2.)/3.);
- if(ichan<0) {
- F1 =-K1fact*BrhoF123(s1,-1);
- F2 = K1fact*BKstarF123(s2,-1);
- F5 =-BKstarF123(q2,-1)*FKrho(s2,s1,-1)*sqrt(2.);
- }
- else if(ichan%8==0) F1 =-K1fact*BrhoF123(s1,ichan/8);
- else if(ichan%8==1) F2 = K1fact*BKstarF123(s2,(ichan-1)/8);
- else F5 = -BKstarF123(q2,ichan/8)*FKrho(s2,s1,(ichan-2)%8)*sqrt(2.);
- }
- // calculate the pi- K0bar pi0
- else if(imode==7) {
- Complex K1fact(K1BreitWigner(q2));
- if(ichan<0) {
- F2 =-K1fact*BrhoF123(s2,-1);
- F5 =-2.*BKstarF123(q2,-1)*FKrho(s1,s2,-1);
- }
- else if(ichan%7==0) F2 =-K1fact*BrhoF123(s2,ichan/7);
- else F5 =-2.*BKstarF123(q2,ichan/7)*FKrho(s1,s2,(ichan-1)%7);
- }
- // calculate the pi- pi0 eta
- else if(imode==8) {
- if(ichan<0) F5 = BrhoF5(q2, -1)*BrhoF123(s3, -1)*sqrt(2./3.);
- else F5 = BrhoF5(q2,ichan/3)*BrhoF123(s3,ichan%3)*sqrt(2./3.);
- }
- // multiply by the prefactors
- using Constants::twopi;
- return FormFactors(F1/_fpi,
- F2/_fpi,
- F3/_fpi,
- F4/_fpi,
- -F5/sqr(twopi)/pow<3,1>(_fpi)
- );
-}
-
-// complete the construction of the decay mode for integration
-bool ThreeMesonDefaultCurrent::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 ) {
- int iq(0),ia(0);
- if(!acceptMode(imode)) return false;
- tPDVector extpart(particles(1,imode,iq,ia));
- Energy min(ZERO);
- for(unsigned int ix=0;ix<extpart.size();++ix) min+=extpart[ix]->massMin();
- if(min>upp) return false;
- // the particles we will use a lot
- tPDPtr a1,k1;
- if(icharge==-3) {
- a1=getParticleData(ParticleID::a_1minus);
- k1=getParticleData(ParticleID::Kstar_1minus);
- }
- else if(icharge==3) {
- a1=getParticleData(ParticleID::a_1plus);
- k1=getParticleData(ParticleID::Kstar_1plus);
- }
- else {
- return false;
- }
- _maxmass=max(_maxmass,upp);
- // the rho0 resonances
- tPDPtr rho0[3] = { getParticleData(113), getParticleData(100113), getParticleData(30113)};
- tPDPtr rhoc[3] = {getParticleData(-213),getParticleData(-100213),getParticleData(-30213)};
- tPDPtr Kstar0[3] = { getParticleData(313), getParticleData(100313), getParticleData(30313)};
- tPDPtr Kstarc[3] = {getParticleData(-323),getParticleData(-100323),getParticleData(-30323)};
- if(icharge==3) {
- for(unsigned int ix=0;ix<3;++ix) {
- rhoc [ix] = rhoc[ix]->CC();
- Kstar0[ix] = Kstar0[ix]->CC();
- Kstarc[ix] = Kstarc[ix]->CC();
- }
- }
- if(imode==0) {
- // channels for pi- pi- pi+
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rho0[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,rho0[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==1) {
- // channels for pi0 pi0 pi-
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rhoc[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,rhoc[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rhoc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==2) {
- // channels for K- pi- K+
- for(unsigned int ix=0;ix<3;++ix) {
- if(Kstar0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,Kstar0[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!rhoc[ix]) continue;
- if(Kstar0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,Kstar0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+2,ires+1,rho0[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- }
- }
- else if(imode==3) {
- // channels for K0 pi- K0bar
- for(unsigned int ix=0;ix<3;++ix) {
- if(Kstarc[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+1,ires+1,Kstarc[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rho0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!rhoc[ix]) continue;
- if(Kstarc[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,Kstarc[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rho0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+2,ires+1,rho0[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- }
- }
- else if(imode==4) {
- // channels for K- pi0 K0
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rhoc[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,a1,ires+1,iloc+2,ires+1,rhoc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==5) {
- // channels for pi0 pi0 K-
- for(unsigned int ix=0;ix<3;++ix) {
- if(!Kstarc[ix]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+1,ires+1,Kstarc[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+2,ires+1,Kstarc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- else if(imode==6) {
- // channels for K- pi- pi+
- for(unsigned int ix=0;ix<3;++ix) {
- if(rho0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+1,ires+1,rho0[ix],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(Kstar0[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+2,ires+1,Kstar0[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!Kstarc[ix]) continue;
- if(rho0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+1,ires+1,rho0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(Kstar0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+2,ires+1,Kstar0[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- }
- }
- else if(imode==7) {
- // channels for pi- kbar0 pi0
- for(unsigned int ix=0;ix<3;++ix) {
- if(rhoc[ix]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,k1,ires+1,iloc+2,ires+1,rhoc[ix],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- for(unsigned int iy=0;iy<3;++iy) {
- if(!Kstarc[ix]) continue;
- if(Kstar0[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+1,ires+1,Kstar0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- if(rhoc[iy]) {
- mode->addChannel((PhaseSpaceChannel(phase),ires,Kstarc[ix],ires+1,iloc+2,ires+1,rhoc[iy],
- ires+2,iloc+1,ires+2,iloc+3));
- }
- }
- }
- }
- else if(imode==8) {
- // channels for pi- pi0 eta
- for(unsigned int ix=0;ix<3;++ix) {
- if(!rhoc[ix]) continue;
- for(unsigned int iy=0;iy<3;++iy) {
- if(!rho0[iy]) continue;
- mode->addChannel((PhaseSpaceChannel(phase),ires,rhoc[ix],ires+1,iloc+1,ires+1,rho0[iy],
- ires+2,iloc+2,ires+2,iloc+3));
- }
- }
- }
- if(_rhoparameters) {
- if(imode!=8) {
- for(unsigned int ix=0;ix<_rhoF123masses.size();++ix) {
- if(rhoc[ix]) mode->resetIntermediate(rhoc[ix],_rhoF123masses[ix],
- _rhoF123widths[ix]);
- if(rho0[ix]) mode->resetIntermediate(rho0[ix],_rhoF123masses[ix],
- _rhoF123widths[ix]);
- }
- }
- else {
- for(unsigned int ix=0;ix<_rhoF5masses.size();++ix) {
- if(rhoc[ix]) mode->resetIntermediate(rhoc[ix],_rhoF5masses[ix],
- _rhoF5widths[ix]);
- if(rho0[ix]) mode->resetIntermediate(rho0[ix],_rhoF5masses[ix],
- _rhoF5widths[ix]);
- }
- }
- }
- // K star parameters in the base class
- if(_kstarparameters) {
- for(unsigned int ix=0;ix<_kstarF123masses.size();++ix) {
- if(Kstarc[ix]) mode->resetIntermediate(Kstarc[ix],_kstarF123masses[ix],
- _kstarF123widths[ix]);
- if(Kstar0[ix]) mode->resetIntermediate(Kstar0[ix],_kstarF123masses[ix],
- _kstarF123widths[ix]);
- }
- }
- return true;
-}
-
-// initialisation of the a_1 width
-// (iopt=-1 initialises, iopt=0 starts the interpolation)
-void ThreeMesonDefaultCurrent::inita1Width(int iopt) {
- if(iopt==-1) {
- _maxcalc=_maxmass;
- if(!_initializea1||_maxmass==ZERO) return;
- // parameters for the table of values
- Energy2 step(sqr(_maxcalc)/199.);
- // integrator to perform the integral
- vector<double> inweights;inweights.push_back(0.5);inweights.push_back(0.5);
- vector<int> intype;intype.push_back(2);intype.push_back(3);
- Energy mrho(getParticleData(ParticleID::rhoplus)->mass()),
- wrho(getParticleData(ParticleID::rhoplus)->width());
- vector<Energy> inmass(2,mrho),inwidth(2,wrho);
- vector<double> inpow(2,0.0);
- ThreeBodyAllOnCalculator<ThreeMesonDefaultCurrent>
- widthgen(inweights,intype,inmass,inwidth,inpow,*this,0,_mpi,_mpi,_mpi);
- // normalisation constant to give physical width if on shell
- double a1const(_a1width/(widthgen.partialWidth(sqr(_a1mass))));
- // loop to give the values
- _a1runq2.clear(); _a1runwidth.clear();
- for(Energy2 moff2(ZERO); moff2<=sqr(_maxcalc); moff2+=step) {
- _a1runwidth.push_back(widthgen.partialWidth(moff2)*a1const);
- _a1runq2.push_back(moff2);
- }
- }
- // set up the interpolator
- else if(iopt==0) {
- _a1runinter = make_InterpolatorPtr(_a1runwidth,_a1runq2,3);
- }
-}
-
-void ThreeMesonDefaultCurrent::dataBaseOutput(ofstream & output,bool header,
- bool create) const {
- if(header) output << "update decayers set parameters=\"";
- if(create) output << "create Herwig::ThreeMesonDefaultCurrent "
- << name() << " HwWeakCurrents.so\n";
- for(unsigned int ix=0;ix<_rhoF123wgts.size();++ix) {
- if(ix<3) output << "newdef ";
- else output << "insert ";
- output << name() << ":F123RhoWeight " << ix << " " << _rhoF123wgts[ix] << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF123wgts.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":F123KstarWeight " << ix << " "
- << _kstarF123wgts[ix] << "\n";
- }
- for(unsigned int ix=0;ix<_rhoF5wgts.size();++ix) {
- if(ix<3) output << "newdef ";
- else output << "insert ";
- output << name() << ":F5RhoWeight " << ix << " " << _rhoF5wgts[ix] << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF5wgts.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":F5KstarWeight " << ix << " " << _kstarF5wgts[ix] << "\n";
- }
- output << "newdef " << name() << ":RhoKstarWgt " << _rhoKstarwgt << "\n";
- output << "newdef " << name() << ":Initializea1 " << _initializea1 << "\n";
- output << "newdef " << name() << ":RhoParameters " << _rhoparameters << "\n";
- output << "newdef " << name() << ":KstarParameters " << _kstarparameters << "\n";
- output << "newdef " << name() << ":a1Parameters " << _a1parameters << "\n";
- output << "newdef " << name() << ":K1Parameters " << _k1parameters << "\n";
- output << "newdef " << name() << ":a1WidthOption " << _a1opt << "\n";
- for(unsigned int ix=0;ix<_a1runwidth.size();++ix) {
- output << "newdef " << name() << ":a1RunningWidth " << ix
- << " " << _a1runwidth[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_a1runq2.size();++ix) {
- output << "newdef " << name() << ":a1RunningQ2 " << ix
- << " " << _a1runq2[ix]/GeV2 << "\n";
- }
- output << "newdef " << name() << ":A1Width " << _a1width/GeV << "\n";
- output << "newdef " << name() << ":A1Mass " << _a1mass/GeV << "\n";
- output << "newdef " << name() << ":K1Width " << _k1width/GeV << "\n";
- output << "newdef " << name() << ":K1Mass " << _k1mass/GeV << "\n";
- output << "newdef " << name() << ":FPi " << _fpi/MeV << "\n";
- for(unsigned int ix=0;ix<_rhoF123masses.size();++ix) {
- if(ix<3) output << "newdef ";
- else output << "insert ";
- output << name() << ":rhoF123masses " << ix
- << " " << _rhoF123masses[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_rhoF123widths.size();++ix) {
- if(ix<3) output << "newdef ";
- else output << "insert ";
- output << name() << ":rhoF123widths " << ix << " "
- << _rhoF123widths[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_rhoF5masses.size();++ix) {
- if(ix<3) output << "newdef ";
- else output << "insert ";
- output << name() << ":rhoF5masses " << ix << " "
- << _rhoF5masses[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_rhoF5widths.size();++ix) {
- if(ix<3) output << "newdef ";
- else output << "insert ";
- output << name() << ":rhoF5widths " << ix << " "
- << _rhoF5widths[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF123masses.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF123masses " << ix << " "
- << _kstarF123masses[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF123widths.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF123widths " << ix << " "
- << _kstarF123widths[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF5masses.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF5masses " << ix << " "
- << _kstarF5masses[ix]/GeV << "\n";
- }
- for(unsigned int ix=0;ix<_kstarF5widths.size();++ix) {
- if(ix<1) output << "newdef ";
- else output << "insert ";
- output << name() << ":KstarF5widths " << ix << " "
- << _kstarF5widths[ix]/GeV << "\n";
- }
- ThreeMesonCurrentBase::dataBaseOutput(output,false,false);
- if(header) output << "\n\" where BINARY ThePEGName=\""
- << fullName() << "\";" << endl;
-}
-
-void ThreeMesonDefaultCurrent::doinitrun() {
- // set up the running a_1 width
- inita1Width(0);
- ThreeMesonCurrentBase::doinitrun();
-}
-
-void ThreeMesonDefaultCurrent::doupdate() {
- ThreeMesonCurrentBase::doupdate();
- // update running width if needed
- if ( !touched() ) return;
- if(_maxmass!=_maxcalc) inita1Width(-1);
-}
-
-Complex ThreeMesonDefaultCurrent::rhoKBreitWigner(Energy2 q2,unsigned int itype,
- unsigned int ires) const {
- Energy q(sqrt(q2)),mass,width,mout[2]={_mpi,_mpi};
- // get the mass and width of the requested resonance
- if(itype==0) {
- mass=_rhoF123masses[ires];
- width=_rhoF123widths[ires];
- }
- else if(itype==1) {
- mass=_rhoF5masses[ires];
- width=_rhoF5widths[ires];
- }
- else if(itype==2) {
- mass=_kstarF123masses[ires];
- width=_kstarF123widths[ires];
- }
- else if(itype==3) {
- mass=_kstarF5masses[ires];
- width=_kstarF5widths[ires];
- }
- else {
- return 0.;
- }
- // calculate the momenta for the running widths
- if(itype>1) mout[0]=_mK;
- Energy pcm0(Kinematics::pstarTwoBodyDecay(mass,mout[0],mout[1]));
- Energy pcm(ZERO);
- if(mout[0]+mout[1]<q){pcm=Kinematics::pstarTwoBodyDecay(q,mout[0],mout[1]);}
- double ratio = Math::Pow<3>(pcm/pcm0);
- Energy gamrun(width*mass*ratio/q);
- Complex ii(0.,1.);
- complex<Energy2> denom(q2-mass*mass+ii*mass*gamrun), numer(-mass*mass);
- return numer/denom;
-}
-
-double ThreeMesonDefaultCurrent::
-threeBodyMatrixElement(const int , const Energy2 q2,
- const Energy2 s3, const Energy2 s2,
- const Energy2 s1, const Energy ,
- const Energy , const Energy ) const {
- Energy2 mpi2(sqr(_mpi));
- Complex propb(BrhoF123(s1,-1)),propa(BrhoF123(s2,-1));
- // the matrix element
- Energy2 output(ZERO);
- // first resonance
- output += ((s1-4.*mpi2) + 0.25*(s3-s2)*(s3-s2)/q2) * real(propb*conj(propb));
- // second resonance
- output += ((s2-4.*mpi2) + 0.25*(s3-s1)*(s3-s1)/q2) * real(propa*conj(propa));
- // the interference term
- output += (0.5*q2-s3-0.5*mpi2+0.25*(s3-s2)*(s3-s1)/q2)*real(propa*conj(propb)+
- propb*conj(propa));
- return output/sqr(_rhoF123masses[0]);
-}
diff --git a/Decay/WeakCurrents/ThreeMesonDefaultCurrent.h b/Decay/WeakCurrents/ThreeMesonDefaultCurrent.h
deleted file mode 100644
--- a/Decay/WeakCurrents/ThreeMesonDefaultCurrent.h
+++ /dev/null
@@ -1,534 +0,0 @@
-// -*- C++ -*-
-//
-// ThreeMesonDefaultCurrent.h is a part of Herwig - A multi-purpose Monte Carlo event generator
-// Copyright (C) 2002-2017 The Herwig Collaboration
-//
-// Herwig is licenced under version 3 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
-#ifndef HERWIG_ThreeMesonDefaultCurrent_H
-#define HERWIG_ThreeMesonDefaultCurrent_H
-//
-// This is the declaration of the ThreeMesonDefaultCurrent class.
-//
-#include "ThreeMesonCurrentBase.h"
-#include "Herwig/Utilities/Interpolator.h"
-#include "Herwig/Utilities/Kinematics.h"
-#include "ThePEG/StandardModel/StandardModelBase.h"
-#include "Herwig/Decay/ResonanceHelpers.h"
-
-namespace Herwig {
-using namespace ThePEG;
-
-/** \ingroup Decay
- *
- * The ThreeMesonDefaultCurrent class implements the currents from Z.Phys.C58:445 (1992),
- * this paper uses the form from Z.Phys.C48:445 (1990) for the \f$a_1\f$ width and
- * is the default model in TAUOLA.
- *
- * The following three meson modes are implemented.
- *
- * - \f$ \pi^- \pi^- \pi^+ \f$, (imode=0)
- * - \f$ \pi^0 \pi^0 \pi^- \f$, (imode=1)
- * - \f$ K^- \pi^- K^+ \f$, (imode=2)
- * - \f$ K^0 \pi^- \bar{K}^0\f$, (imode=3)
- * - \f$ K^- \pi^0 K^0 \f$, (imode=4)
- * - \f$ \pi^0 \pi^0 K^- \f$, (imode=5)
- * - \f$ K^- \pi^- \pi^+ \f$, (imode=6)
- * - \f$ \pi^- \bar{K}^0 \pi^0 \f$, (imode=7)
- * - \f$ \pi^- \pi^0 \eta \f$, (imode=8)
- *
- * using the currents from TAUOLA
- *
- *
- * @see ThreeMesonCurrentBase,
- * @see WeakCurrent.
- * @see Defaulta1MatrixElement
- *
- */
-class ThreeMesonDefaultCurrent: public ThreeMesonCurrentBase {
-
- /**
- * The matrix element for the running \f$a_1\f$ width is a friend to
- * keep some members private.
- */
- friend class Defaulta1MatrixElement;
-
-public:
-
- /**
- * Default constructor
- */
- ThreeMesonDefaultCurrent();
-
- /** @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);
- //@}
-
- /**
- * Standard Init function used to initialize the interfaces.
- */
- static void Init();
-
-public:
-
- /** @name Methods for the construction of the phase space integrator. */
- //@{
- /**
- * Complete the construction of the decay mode for integration.classes inheriting
- * from this one.
- * This method is purely virtual and must be implemented in the classes inheriting
- * from WeakCurrent.
- * @param icharge The total charge of the outgoing particles in the current.
- * @param resonance If specified only include terms with this particle
- * @param Itotal If specified the total isospin of the current
- * @param I3 If specified the thrid component of isospin
- * @param imode The mode in the current being asked for.
- * @param mode The phase space mode for the integration
- * @param iloc The location of the of the first particle from the current in
- * the list of outgoing particles.
- * @param ires The location of the first intermediate for the current.
- * @param phase The prototype phase space channel for the integration.
- * @param upp The maximum possible mass the particles in the current are
- * allowed to have.
- * @return Whether the current was sucessfully constructed.
- */
- virtual bool createMode(int icharge, tcPDPtr resonance,
- IsoSpin::IsoSpin Itotal, IsoSpin::I3 i3,
- unsigned int imode,PhaseSpaceModePtr mode,
- unsigned int iloc,int ires,
- PhaseSpaceChannel phase, Energy upp );
- //@}
-
- /**
- * 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;
-
- /**
- * the matrix element for the \f$a_1\f$ decay to calculate the running width
- * @param imode The mode for which the matrix element is needed.
- * @param q2 The mass of the decaying off-shell \f$a_1\f$, \f$q^2\f$.
- * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
- * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
- * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
- * @param m1 The mass of the first outgoing particle.
- * @param m2 The mass of the second outgoing particle.
- * @param m3 The mass of the third outgoing particle.
- * @return The matrix element squared summed over spins.
- */
- double threeBodyMatrixElement(const int imode, const Energy2 q2,
- const Energy2 s3, const Energy2 s2,
- const Energy2 s1, const Energy m1,
- const Energy m2, const Energy m3) const;
-
-protected:
-
- /**
- * Can the current handle a particular set of mesons.
- * As this current includes all the allowed modes this is always true.
- */
- virtual bool acceptMode(int) const;
-
- /**
- * Calculate the form factor for the current. Implements the form factors
- * described above.
- * @param ichan The phase space channel
- * @param imode The mode
- * @param q2 The scale \f$q^2\f$ for the current.
- * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
- * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
- * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
- */
- virtual FormFactors calculateFormFactors(const int ichan, const int imode,
- Energy2 q2,
- Energy2 s1, Energy2 s2, Energy2 s3) const;
-
-protected:
-
- /** @name Clone Methods. */
- //@{
- /**
- * Make a simple clone of this object.
- * @return a pointer to the new object.
- */
- virtual IBPtr clone() const {return new_ptr(*this);}
-
- /** 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 {return new_ptr(*this);}
- //@}
-
-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();
-
- /**
- * Initialize this object to the begining of the run phase.
- */
- virtual void doinitrun();
-
- /**
- * Check sanity of the object during the setup phase.
- */
- virtual void doupdate();
- //@}
-
-private:
-
- /**
- * Private and non-existent assignment operator.
- */
- ThreeMesonDefaultCurrent & operator=(const ThreeMesonDefaultCurrent &);
-
-private:
-
- /**
- * The \f$\rho\f$ Breit-Wigner for the \f$F_{1,2,3}\f$ form factors.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @param ires Which \f$\rho\f$ multiplet
- * @return The Breit-Wigner
- */
- Complex BrhoF123(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
- norm+=_rhoF123wgts[ix];
- }
- if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_rhoF123wgts.size()));ix<N;++ix) {
- output+=_rhoF123wgts[ix]*rhoKBreitWigner(q2,0,ix);
- }
- }
- else {
- unsigned int temp(ires);
- if(temp<_rhoF123wgts.size()&&temp<3)
- output=_rhoF123wgts[temp]*rhoKBreitWigner(q2,0,temp);
- else
- output=0.;
- }
- return output/norm;
- }
-
- /**
- * The \f$\rho\f$ Breit-Wigner for the \f$F_5\f$ form factors.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @param ires Which \f$\rho\f$ multiplet
- * @return The Breit-Wigner
- */
- Complex BrhoF5(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
- norm+=_rhoF5wgts[ix];
- }
- if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_rhoF5wgts.size()));ix<N;++ix) {
- output+=_rhoF5wgts[ix]*rhoKBreitWigner(q2,1,ix);
- }
- }
- else {
- unsigned int temp(ires);
- if(temp<_rhoF5wgts.size()&&temp<3) {
- output=_rhoF5wgts[temp]*rhoKBreitWigner(q2,1,temp);
- }
- }
- return output/norm;
- }
-
- /**
- * The \f$K^*\f$ Breit-Wigner for the \f$F_{1,2,3}\f$ form factors.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @param ires Which \f$\rho\f$ multiplet
- * @return The Breit-Wigner
- */
- Complex BKstarF123(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
- norm+=_kstarF123wgts[ix];
- }
- if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_kstarF123wgts.size()));ix<N;++ix) {
- output+=_kstarF123wgts[ix]*rhoKBreitWigner(q2,2,ix);
- }
- }
- else {
- unsigned int temp(ires);
- if(temp<_kstarF123wgts.size()&&temp<3) {
- output=_kstarF123wgts[temp]*rhoKBreitWigner(q2,2,temp);
- }
- }
- return output/norm;
- }
-
- /**
- * The \f$K^*\f$ Breit-Wigner for the \f$F_5\f$ form factors.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @param ires Which \f$\rho\f$ multiplet
- * @return The Breit-Wigner
- */
- Complex BKstarF5(Energy2 q2,int ires) const {
- Complex output(0.),norm(0.);
- for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
- norm+=_kstarF5wgts[ix];
- }
- if(ires<0) {
- for(unsigned int ix=0,N=min(3,int(_kstarF5wgts.size()));ix<N;++ix) {
- output+=_kstarF5wgts[ix]*rhoKBreitWigner(q2,3,ix);
- }
- }
- else {
- unsigned int temp(ires);
- if(temp<_kstarF5wgts.size()&&temp<3) {
- output=_kstarF5wgts[ires]*rhoKBreitWigner(q2,3,temp);
- }
- }
- return output/norm;
- }
-
- /**
- * Mixed Breit Wigner for the \f$F_5\f$ form factor
- * @param si The scale \f$s_1\f$.
- * @param sj The scale \f$s_2\f$.
- * @param ires Which resonances to use
- * @return The mixed Breit-Wigner
- */
- Complex FKrho(Energy2 si,Energy2 sj,int ires) const {
- Complex output;
- if(ires<0){output = _rhoKstarwgt*BKstarF123(si,-1)+BrhoF123(sj,-1);}
- else if(ires%2==0){output= _rhoKstarwgt*BKstarF123(si,ires/2);}
- else if(ires%2==1){output=BrhoF123(sj,ires/2);}
- output /=(1.+_rhoKstarwgt);
- return output;
- }
-
- /**
- * \f$a_1\f$ Breit-Wigner
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @return The Breit-Wigner
- */
- Complex a1BreitWigner(Energy2 q2) const {
- if(!_a1opt)
- return Resonance::BreitWignera1(q2,_a1mass,_a1width);
- Complex ii(0.,1.);
- Energy2 m2(_a1mass*_a1mass);
- Energy q(sqrt(q2));
- Energy width = (*_a1runinter)(q2);
- return m2/(m2-q2-ii*q*width);
- }
-
- /**
- * The \f$K_1\f$ Breit-Wigner
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner
- * @return The Breit-Wigner
- */
- Complex K1BreitWigner(Energy2 q2) const {
- Energy2 m2 = sqr(_k1mass);
- Complex ii(0.,1.);
- complex<Energy2> fact(m2 - ii*_k1mass*_k1width);
- return fact/(fact-q2);
- }
-
- /**
- * Initialize the \f$a_1\f$ running width
- * @param iopt Initialization option (-1 full calculation, 0 set up the interpolation)
- */
- void inita1Width(int iopt);
-
- /**
- * Breit-Wigners for the \f$\rho\f$ and \f$K^*\f$.
- * @param q2 The scale \f$q^2\f$ for the Breit-Wigner.
- * @param itype The type of Breit-Wigner, \e i.e. which masses and widths to use.x
- * @param ires Which multiplet to use.
- */
- Complex rhoKBreitWigner(Energy2 q2,unsigned int itype,unsigned int ires) const;
-
-private:
-
- /**
- * Parameters for the \f$\rho\f$ Breit-Wigner in the
- * \f$F_{1,2,3}\f$ form factors.
- */
- vector<double> _rhoF123wgts;
-
- /**
- * Parameters for the \f$K^*\f$ Breit-Wigner in the
- * \f$F_{1,2,3}\f$ form factors.
- */
- vector<double> _kstarF123wgts;
-
- /**
- * Parameters for the \f$\rho\f$ Breit-Wigner in the
- * \f$F_5\f$ form factors.
- */
- vector<double> _rhoF5wgts;
-
- /**
- * Parameters for the \f$K^*\f$ Breit-Wigner in the
- * \f$F_5\f$ form factors.
- */
- vector<double> _kstarF5wgts;
-
- /**
- * The relative weight of the \f$\rho\f$ and \f$K^*\f$ where needed.
- */
- double _rhoKstarwgt;
-
- /**
- * The \f$a_1\f$ width for the running \f$a_1\f$ width calculation.
- */
- vector<Energy> _a1runwidth;
-
- /**
- * The \f$q^2\f$ for the running \f$a_1\f$ width calculation.
- */
- vector<Energy2> _a1runq2;
-
-
- /**
- * The interpolator for the running \f$a_1\f$ width calculation.
- */
- Interpolator<Energy,Energy2>::Ptr _a1runinter;
-
- /**
- * Initialize the running \f$a_1\f$ width.
- */
- bool _initializea1;
-
- /**
- * The mass of the \f$a_1\f$ resonances.
- */
- Energy _a1mass;
-
- /**
- * The width of the \f$a_1\f$ resonances.
- */
- Energy _a1width;
-
- /**
- * The mass of the \f$aK1\f$ resonances.
- */
- Energy _k1mass;
-
- /**
- * The width of the \f$K_1\f$ resonances.
- */
- Energy _k1width;
-
- /**
- * The pion decay constant, \f$f_\pi\f$.
- */
- Energy _fpi;
-
- /**
- * The pion mass
- */
- Energy _mpi;
-
- /**
- * The kaon mass
- */
- Energy _mK;
-
- /**
- * use local values of the \f$\rho\f$ masses and widths
- */
- bool _rhoparameters;
-
- /**
- * The \f$\rho\f$ masses for the \f$F_{1,2,3}\f$ form factors.
- */
- vector<Energy> _rhoF123masses;
-
- /**
- * The \f$\rho\f$ masses for the \f$F_5\f$ form factors.
- */
- vector<Energy> _rhoF5masses;
-
- /**
- * The \f$\rho\f$ widths for the \f$F_{1,2,3}\f$ form factors.
- */
- vector<Energy> _rhoF123widths;
-
- /**
- * The \f$\rho\f$ widths for the \f$F_5\f$ form factors.
- */
- vector<Energy> _rhoF5widths;
-
- /**
- * use local values of the \f$K^*\f$ resonances masses and widths
- */
- bool _kstarparameters;
-
- /**
- * The \f$K^*\f$ masses for the \f$F_{1,2,3}\f$ form factors.
- */
- vector<Energy> _kstarF123masses;
-
- /**
- * The \f$K^*\f$ masses for the \f$F_5\f$ form factors.
- */
- vector<Energy> _kstarF5masses;
-
- /**
- * The \f$K^*\f$ widths for the \f$F_{1,2,3}\f$ form factors.
- */
- vector<Energy> _kstarF123widths;
-
- /**
- * The \f$K^*\f$ widths for the \f$F_5\f$ form factors.
- */
- vector<Energy> _kstarF5widths;
-
- /**
- * Use local values of the \f$a_1\f$ parameters
- */
- bool _a1parameters;
-
- /**
- * Use local values of the \f$K_1\f$ parameters
- */
- bool _k1parameters;
-
- /**
- * Option for the \f$a_1\f$ width
- */
- bool _a1opt;
-
- /**
- * The maximum mass of the hadronic system
- */
- Energy _maxmass;
-
- /**
- * The maximum mass when the running width was calculated
- */
- Energy _maxcalc;
-
-};
-
-}
-
-#endif /* THEPEG_ThreeMesonDefaultCurrent_H */
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:32 PM (1 d, 4 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242484
Default Alt Text
(141 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment