Page MenuHomeHEPForge
Files F8309972

No OneTemporary
Actions

View Options

diff --git a/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/OneOneOneEWSplitFn.cc
@@ -1,279 +1,284 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the OneOneOneEWSplitFn class.
//
#include "OneOneOneEWSplitFn.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/ParticleData.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "Herwig/Models/StandardModel/SMFFHVertex.h"
#include "ThePEG/Interface/Parameter.h"
using namespace Herwig;
IBPtr OneOneOneEWSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr OneOneOneEWSplitFn::fullclone() const {
return new_ptr(*this);
}
void OneOneOneEWSplitFn::persistentOutput(PersistentOStream & os) const {
os << gWWG_ << gWWZ_ << _theSM << _couplingValueIm << _couplingValueRe;
}
void OneOneOneEWSplitFn::persistentInput(PersistentIStream & is, int) {
is >> gWWG_ >> gWWZ_ >> _theSM >> _couplingValueIm >> _couplingValueRe;
}
// The following static variable is needed for the type description system in ThePEG.
DescribeClass<OneOneOneEWSplitFn,SplittingFunction>
describeHerwigOneOneOneEWSplitFn("Herwig::OneOneOneEWSplitFn", "HwShower.so");
void OneOneOneEWSplitFn::Init() {
static ClassDocumentation<OneOneOneEWSplitFn> documentation
("The OneOneOneEWSplitFn class implements the splitting W->WG, W->WZ and Z->ZZ");
static Parameter<OneOneOneEWSplitFn, double> interfaceCouplingValueIm
("CouplingValue.Im",
"The numerical value (imaginary part) of the splitting coupling to be imported for BSM splittings",
&OneOneOneEWSplitFn::_couplingValueIm, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
static Parameter<OneOneOneEWSplitFn, double> interfaceCouplingValueRe
("CouplingValue.Re",
"The numerical value (real part) of the splitting coupling to be imported for BSM splittings",
&OneOneOneEWSplitFn::_couplingValueRe, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
}
void OneOneOneEWSplitFn::doinit() {
SplittingFunction::doinit();
tcSMPtr sm = generator()->standardModel();
double sw2 = sm->sin2ThetaW();
// WWZ coupling
gWWZ_ = sqrt((1.-sw2)/sw2);
// WWG coupling
gWWG_ = 1.;
// to employ running masses, wherever needed
_theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
}
void OneOneOneEWSplitFn::getCouplings(Complex & gvvv, const IdList & ids) const {
if(_couplingValueIm!=0||_couplingValueRe!=0) {
double e = sqrt(4.*Constants::pi*generator()->standardModel()
->alphaEM(sqr(getParticleData(ParticleID::Z0)->mass())));
gvvv = Complex(_couplingValueRe,_couplingValueIm)/e;
}
// Z > WW
else if(ids[0]->id()==ParticleID::Z0 && abs(ids[1]->id())==ParticleID::Wplus
&& abs(ids[2]->id())==ParticleID::Wplus){
gvvv = Complex(0.,gWWZ_);
}
// W > WG
else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
&& ids[2]->id()==ParticleID::gamma){
gvvv = Complex(0.,gWWG_);
}
// W > WZ
else if(abs(ids[0]->id())==ParticleID::Wplus && abs(ids[1]->id())==ParticleID::Wplus
&& ids[2]->id()==ParticleID::Z0){
gvvv = Complex(0.,gWWZ_);
}
else
assert(false);
}
double OneOneOneEWSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
Complex gvvv(0.,0.);
getCouplings(gvvv,ids);
double abs_rho_00 = abs(rho(0,0));
double abs_rho_11 = abs(rho(1,1));
double abs_rho_22 = abs(rho(2,2));
// massless limit
double val = ((2.*sqr(1.-(1.-z)*z))/((1.-z)*z))*(abs_rho_00+abs_rho_22);
// massive limits
if(mass) {
double m0t2 = sqr(ids[0]->mass())/t;
double m1t2 = sqr(ids[1]->mass())/t;
double m2t2 = sqr(ids[2]->mass())/t;
val += 4.*m0t2*sqr(1.-z)*abs_rho_11 + (2*(m0t2*sqr(1.-(1.-z)*z)
-m2t2*(1.-sqr(1.-z)*z) -m1t2*(1.-(1.-z)*sqr(z)))*(abs_rho_00+abs_rho_22))
/((1.-z)*z);
}
return norm(gvvv)*val;
}
double OneOneOneEWSplitFn::overestimateP(const double z,
const IdList & ids) const {
Complex gvvv(0.,0.);
getCouplings(gvvv,ids);
double val = norm(gvvv)*(2./(z*(1.-z)));
return val;
}
double OneOneOneEWSplitFn::ratioP(const double z, const Energy2 t,
const IdList & ids, const bool mass,
const RhoDMatrix & rho) const {
double val(0.);
double abs_rho_00 = abs(rho(0,0));
double abs_rho_11 = abs(rho(1,1));
double abs_rho_22 = abs(rho(2,2));
// massless limit
val = sqr(1.-(1.-z)*z)*(abs_rho_00+abs_rho_22);
// massive limit
if(mass) {
double m0t2 = sqr(ids[0]->mass())/t;
double m1t2 = sqr(ids[1]->mass())/t;
double m2t2 = sqr(ids[2]->mass())/t;
val += (4.*m0t2*sqr(1.-z)*abs_rho_11 + (2*(m0t2*sqr(1.-(1.-z)*z)
-m2t2*(1.-sqr(1.-z)*z) -m1t2*(1.-(1.-z)*sqr(z)))*(abs_rho_00+abs_rho_22))
/((1.-z)*z))/(2./(z*(1.-z)));
}
return val;
}
double OneOneOneEWSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
Complex gvvv(0.);
getCouplings(gvvv,ids);
double pre = norm(gvvv);
switch (PDFfactor) {
case 0:
return 2.*pre*(log(z)-log(1.-z));
case 1:
//return -2.*pre*(1./z+log(1.-z)-log(z));
case 2:
//return 2.*pre*(2.*log(z)+(2.*z-1.)/(z*(1.-z))-2.*log(1.-z));
case 3:
//return 2.*pre*(1./(1.-z)-1./z-2.*log(1.-z)+2.*log(z));
default:
throw Exception() << "OneOneOneEWSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double OneOneOneEWSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
Complex gvvv(0.);
getCouplings(gvvv,ids);
double pre = norm(gvvv);
switch (PDFfactor) {
case 0:
return exp(0.5*r/pre)/(1.+exp(0.5*r/pre));
case 1:
case 2:
case 3:
default:
throw Exception() << "OneOneOneEWSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool OneOneOneEWSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3) return false;
- // Z > WW
- if(ids[0]->id()==ParticleID::Z0 && abs(ids[1]->id())==ParticleID::Wplus
- && ids[1]->id()==-ids[2]->id())
- return true;
-
- if(abs(ids[0]->id())==ParticleID::Wplus) {
- // W > WG
- if(ids[1]->id()==ids[0]->id() && ids[2]->id()==ParticleID::gamma)
+ if(_couplingValueIm==0&&_couplingValueRe==0) {
+ // Z > WW
+ if(ids[0]->id()==ParticleID::Z0 && abs(ids[1]->id())==ParticleID::Wplus
+ && ids[1]->id()==-ids[2]->id())
return true;
- // W > WZ
- if(ids[1]->id()==ids[0]->id() && ids[2]->id()==ParticleID::Z0)
+ if(abs(ids[0]->id())==ParticleID::Wplus) {
+ // W > WG
+ if(ids[1]->id()==ids[0]->id() && ids[2]->id()==ParticleID::gamma)
+ return true;
+ // W > WZ
+ if(ids[1]->id()==ids[0]->id() && ids[2]->id()==ParticleID::Z0)
+ return true;
+ }
+ }
+ else {
+ if(ids[0]->iSpin()==PDT::Spin1 && ids[1]->iSpin()==PDT::Spin1 && ids[2]->iSpin()==PDT::Spin1)
return true;
}
return false;
}
vector<pair<int, Complex> >
OneOneOneEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
vector<pair<int, Complex> >
OneOneOneEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
DecayMEPtr OneOneOneEWSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool) {
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin1)));
Complex gvvv(0.);
getCouplings(gvvv,ids);
// defining dummies
double m0t = ids[0]->mass()/sqrt(t);
double m1t = ids[1]->mass()/sqrt(t);
double m2t = ids[2]->mass()/sqrt(t);
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
double z1_z = z*(1.-z);
double sqrtmass = sqrt(sqr(m0t)-sqr(m1t)/z-sqr(m2t)/(1.-z)+1.);
double r2 = sqrt(2.);
// assign kernel
(*kernal)(0,0,0) = gvvv*phase*(1./sqrt(z1_z))*sqrtmass;
(*kernal)(0,0,1) = gvvv*r2*m2t*(z/(1.-z)); //2>4
(*kernal)(0,0,2) = -gvvv*cphase*sqrt(z/(1.-z))*sqrtmass;
(*kernal)(0,1,0) = -gvvv*r2*m1t*(1.-z)/z; //2>4
(*kernal)(0,1,1) = 0.;
(*kernal)(0,1,2) = 0.;
(*kernal)(0,2,0) = -gvvv*(1.-z)*cphase*sqrt((1.-z)/z)*sqrtmass;
(*kernal)(0,2,1) = 0.;
(*kernal)(0,2,2) = 0.;
(*kernal)(1,0,0) = 0.;
(*kernal)(1,0,1) = 0.; //2>4
(*kernal)(1,0,2) = -gvvv*r2*m0t*(1.-z); //2>4
(*kernal)(1,1,0) = 0.; //221>441
(*kernal)(1,1,1) = 0.; //222>444
(*kernal)(1,1,2) = 0.; //223>443
(*kernal)(1,2,0) = -gvvv*r2*m0t*(1.-z); //2>4
(*kernal)(1,2,1) = 0.; //2>4
(*kernal)(1,2,2) = 0.; //2>4
(*kernal)(2,0,0) = 0.;
(*kernal)(2,0,1) = 0.;
(*kernal)(2,0,2) = gvvv*(1.-z)*phase*sqrt((1.-z)/z)*sqrtmass;
(*kernal)(2,1,0) = 0.;
(*kernal)(2,1,1) = 0.; //2>4
(*kernal)(2,1,2) = -gvvv*r2*m1t*((1.-z)/z);//2>4
(*kernal)(2,2,0) = gvvv*phase*sqrt(z/(1.-z))*sqrtmass;
(*kernal)(2,2,1) = gvvv*r2*m2t*(z/(1.-z)); //2>4
(*kernal)(2,2,2) = -gvvv*cphase*(1./sqrt(z1_z))*sqrtmass;
// return the answer
return kernal;
}
diff --git a/Shower/QTilde/SplittingFunctions/OneOneZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/OneOneZeroEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/OneOneZeroEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/OneOneZeroEWSplitFn.cc
@@ -1,238 +1,239 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the OneOneZeroEWSplitFn class.
//
#include "OneOneZeroEWSplitFn.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/ParticleData.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "Herwig/Models/StandardModel/SMFFHVertex.h"
#include "ThePEG/Interface/Parameter.h"
using namespace Herwig;
IBPtr OneOneZeroEWSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr OneOneZeroEWSplitFn::fullclone() const {
return new_ptr(*this);
}
void OneOneZeroEWSplitFn::persistentOutput(PersistentOStream & os) const {
os << gWWH_ << gZZH_ << _theSM << _couplingValueIm << _couplingValueRe;
}
void OneOneZeroEWSplitFn::persistentInput(PersistentIStream & is, int) {
is >> gWWH_ >> gZZH_ >> _theSM >> _couplingValueIm >> _couplingValueRe;
}
// The following static variable is needed for the type description system in ThePEG.
DescribeClass<OneOneZeroEWSplitFn,SplittingFunction>
describeHerwigOneOneZeroEWSplitFn("Herwig::OneOneZeroEWSplitFn", "HwShower.so");
void OneOneZeroEWSplitFn::Init() {
static ClassDocumentation<OneOneZeroEWSplitFn> documentation
("The OneOneZeroEWSplitFn class implements the splittings W->WH and Z->ZH");
static Parameter<OneOneZeroEWSplitFn, double> interfaceCouplingValueIm
("CouplingValue.Im",
"The numerical value (imaginary part) of the splitting coupling to be imported for BSM splittings",
&OneOneZeroEWSplitFn::_couplingValueIm, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
static Parameter<OneOneZeroEWSplitFn, double> interfaceCouplingValueRe
("CouplingValue.Re",
"The numerical value (real part) of the splitting coupling to be imported for BSM splittings",
&OneOneZeroEWSplitFn::_couplingValueRe, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
}
void OneOneZeroEWSplitFn::doinit() {
SplittingFunction::doinit();
tcSMPtr sm = generator()->standardModel();
double sw2 = sm->sin2ThetaW();
// W -> W H coupling g = e/sin theta_W
gWWH_ = 1./sqrt(sw2);
// Z -> Z H coupling
gZZH_ = 1./(sqrt(sw2)*sqrt(1.-sw2));
// to employ running masses, wherever needed
_theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
}
void OneOneZeroEWSplitFn::getCouplings(Complex & g, const IdList & ids) const {
if(_couplingValueIm!=0||_couplingValueRe!=0) {
double e = sqrt(4.*Constants::pi*generator()->standardModel()
->alphaEM(sqr(getParticleData(ParticleID::Z0)->mass())));
Energy m0 = ids[0]->mass();
g = Complex(_couplingValueRe,_couplingValueIm)*GeV/e/m0;
}
else {
if(abs(ids[0]->id())==ParticleID::Wplus){
g = Complex(0.,gWWH_);
}
else if(ids[0]->id()==ParticleID::Z0){
g = Complex(0.,gZZH_);
}
else
assert(false);
}
}
double OneOneZeroEWSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
Complex gvvh(0.,0.);
getCouplings(gvvh,ids);
double rho00 = abs(rho(0,0));
double rho11 = abs(rho(1,1));
double rho22 = abs(rho(2,2));
// the splitting in the massless limit
double val = ((1.-z)*(2.*rho11+sqr(z)*(rho00+rho22)))/(4.*z);
// the massive limit
if(mass){
// get the running mass
double mBt = _theSM->mass(t,getParticleData(ids[0]->id()))/sqrt(t);
double mHt = _theSM->mass(t,getParticleData(ids[2]->id()))/sqrt(t);
val += -(sqr(mHt)*(2.*rho11+sqr(z)*(rho00+rho22)))/(4.*z)
- (sqr(mBt)*(2.*rho11+z*(-4.*rho11+z*(2.*rho11+(-1.+(-2.+z)*z)*(rho00+rho22)))))
/(4.*sqr(z));
}
return norm(gvvh)*val;
}
double OneOneZeroEWSplitFn::overestimateP(const double z,
const IdList & ids) const {
Complex gvvh(0.,0.);
getCouplings(gvvh,ids);
return norm(gvvh)/(2.*z);
}
double OneOneZeroEWSplitFn::ratioP(const double z, const Energy2 t,
const IdList & ids, const bool mass,
const RhoDMatrix & rho) const {
double rho00 = abs(rho(0,0));
double rho11 = abs(rho(1,1));
double rho22 = abs(rho(2,2));
// ratio in the massless limit
double val = ((1.-z)*(2.*rho11+sqr(z)*(rho00+rho22)))/2.;
// the massive limit
if(mass){
// get the running mass
double mBt = _theSM->mass(t,getParticleData(ids[0]->id()))/sqrt(t);
double mHt = _theSM->mass(t,getParticleData(ids[2]->id()))/sqrt(t);
val += -(sqr(mHt)*(2.*rho11+ sqr(z)*(rho00+rho22)))/2.
- (sqr(mBt)*(2.*rho11+z*(-4.*rho11+2.*z*rho11+z*(-1.+(-2.+z)*z)*(rho00+rho22))))/(2.*z);
}
return val;
}
double OneOneZeroEWSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
Complex gvvh(0.,0.);
getCouplings(gvvh,ids);
double pre = norm(gvvh);
switch (PDFfactor) {
case 0:
return pre*log(z)/2.;
case 1:
case 2:
case 3:
default:
throw Exception() << "OneOneZeroEWSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double OneOneZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
Complex gvvh(0.,0.);
getCouplings(gvvh,ids);
double pre = norm(gvvh);
switch (PDFfactor) {
case 0:
return exp(2.*r/(pre));
case 1:
case 2:
case 3:
default:
throw Exception() << "OneOneZeroEWSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool OneOneZeroEWSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3)
return false;
- if(ids[0]->id()!=ids[1]->id())
- return false;
- if(ids[0]->iSpin()==PDT::Spin1 && ids[1]->iSpin()==PDT::Spin1 && ids[2]->iSpin()==PDT::Spin0)
- return true;
- //if(abs(ids[0]->id())==ParticleID::Wplus && ids[2]->iSpin()==PDT::Spin0)
- // return true;
- //else if(ids[0]->id()==ParticleID::Z0 && ids[2]->iSpin()==PDT::Spin0)
- // return true;
- else
- return false;
+ if(_couplingValueIm==0&&_couplingValueRe==0) {
+ if(abs(ids[0]->id())==ParticleID::Wplus && ids[2]->iSpin()==PDT::Spin0)
+ return true;
+ else if(ids[0]->id()==ParticleID::Z0 && ids[2]->iSpin()==PDT::Spin0)
+ return true;
+ }
+ else {
+ if(ids[0]->iSpin()==PDT::Spin1 && ids[1]->iSpin()==PDT::Spin1 && ids[2]->iSpin()==PDT::Spin0)
+ return true;
+ }
+ return false;
}
vector<pair<int, Complex> >
OneOneZeroEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
vector<pair<int, Complex> >
OneOneZeroEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
DecayMEPtr OneOneZeroEWSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool) {
Complex gvvh(0.,0.);
getCouplings(gvvh,ids);
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin0)));
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
double r2 = sqrt(2.);
double mBt = _theSM->mass(t,getParticleData(ids[0]->id()))/sqrt(t);
double mHt = _theSM->mass(t,getParticleData(ids[2]->id()))/sqrt(t);
double sqrtmass = sqrt(-(sqr(mBt)*sqr(1.-z))-sqr(mHt)*z+(1.-z)*z);
// assign kernel
(*kernal)(0,0,0) = -gvvh*mBt/r2; // 111
(*kernal)(0,1,0) = -(cphase/2.)*sqrtmass; // 121
(*kernal)(0,2,0) = 0.; // 131
(*kernal)(1,0,0) = ( phase/2.*z)*sqrtmass; // 211
(*kernal)(1,1,0) = 0.; // 221 > 441
(*kernal)(1,2,0) = -(cphase/2.*z)*sqrtmass; // 231
(*kernal)(2,0,0) = 0.; // 311
(*kernal)(2,1,0) = ( phase/2.)*sqrtmass;; // 321
(*kernal)(2,2,0) = -gvvh*mBt/r2; // 331
// return the answer
return kernal;
}
diff --git a/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/OneZeroZeroEWSplitFn.cc
@@ -1,216 +1,213 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the OneZeroZeroEWSplitFn class.
//
#include "OneZeroZeroEWSplitFn.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/ParticleData.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "Herwig/Models/StandardModel/SMFFHVertex.h"
#include "ThePEG/Interface/Parameter.h"
using namespace Herwig;
IBPtr OneZeroZeroEWSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr OneZeroZeroEWSplitFn::fullclone() const {
return new_ptr(*this);
}
void OneZeroZeroEWSplitFn::persistentOutput(PersistentOStream & os) const {
os << _couplingValueIm << _couplingValueRe;
}
void OneZeroZeroEWSplitFn::persistentInput(PersistentIStream & is, int) {
is >> _couplingValueIm >> _couplingValueRe;
}
// The following static variable is needed for the type description system in ThePEG.
DescribeClass<OneZeroZeroEWSplitFn,SplittingFunction>
describeHerwigOneZeroZeroEWSplitFn("Herwig::OneZeroZeroEWSplitFn", "HwShower.so");
void OneZeroZeroEWSplitFn::Init() {
static ClassDocumentation<OneZeroZeroEWSplitFn> documentation
("The OneZeroZeroEWSplitFn class implements purly beyond SM electroweak splittings V->HH'");
static Parameter<OneZeroZeroEWSplitFn, double> interfaceCouplingValueIm
("CouplingValue.Im",
"The numerical value (imaginary part) of the splitting coupling to be imported for BSM splittings",
&OneZeroZeroEWSplitFn::_couplingValueIm, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
static Parameter<OneZeroZeroEWSplitFn, double> interfaceCouplingValueRe
("CouplingValue.Re",
"The numerical value (real part) of the splitting coupling to be imported for BSM splittings",
&OneZeroZeroEWSplitFn::_couplingValueRe, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
}
void OneZeroZeroEWSplitFn::doinit() {
SplittingFunction::doinit();
}
void OneZeroZeroEWSplitFn::getCouplings(Complex & g, const IdList &) const {
g = Complex(_couplingValueRe,_couplingValueIm);
}
double OneZeroZeroEWSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
Complex gvhh(0.,0.);
getCouplings(gvhh,ids);
double rho00 = abs(rho(0,0));
double rho11 = abs(rho(1,1));
double rho22 = abs(rho(2,2));
// the splitting in the massless limit
double val = (1.-z)*z*(rho00+rho22);
// the massive limit
if(mass){
// get the running mass
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
val += (m0t2*sqr(-1.+2.*z)*rho11)/2.+(-(m1t2*(1.-z))-m2t2*z+m0t2*(1.-z)*z)
*(rho00+rho22);
}
return norm(gvhh)*val;
}
double OneZeroZeroEWSplitFn::overestimateP(const double z,
const IdList & ids) const {
Complex gvhh(0.,0.);
getCouplings(gvhh,ids);
return norm(gvhh)*(1.-z)*z;
}
double OneZeroZeroEWSplitFn::ratioP(const double z, const Energy2 t,
const IdList & ids, const bool mass,
const RhoDMatrix & rho) const {
double rho00 = abs(rho(0,0));
double rho11 = abs(rho(1,1));
double rho22 = abs(rho(2,2));
// ratio in the massless limit
double val = rho00+rho22;
// the massive limit
if(mass){
// get the running mass
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
val += ((m0t2*sqr(-1.+2.*z)*rho11)/2.+(-(m1t2*(1.-z))-m2t2*z+m0t2*(1.-z)*z)
*(rho00+rho22))/((1.-z)*z);
}
return val;
}
double OneZeroZeroEWSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
Complex gvhh(0.,0.);
getCouplings(gvhh,ids);
double pre = norm(gvhh);
switch (PDFfactor) {
case 0:
return pre/6.*(3.-2.*z)*sqr(z);
case 1:
case 2:
case 3:
default:
throw Exception() << "OneZeroZeroEWSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double OneZeroZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
Complex gvhh(0.,0.);
getCouplings(gvhh,ids);
double pre = norm(gvhh);
switch (PDFfactor) {
case 0:
return (1.-pre/pow(-pow(pre,3)+12.*pow(pre,2)*r+2.*sqrt(6.)
*sqrt(-(pow(pre,4)*(pre-6.*r)*r)),1./3.)- pow(-pow(pre,3)+12.
*pow(pre,2)*r+2.*sqrt(6)*sqrt(-(pow(pre,4)*(pre-6.*r)*r)),1./3.)/pre)/2.;
case 1:
case 2:
case 3:
default:
throw Exception() << "OneZeroZeroEWSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool OneZeroZeroEWSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3)
return false;
if(_couplingValueIm==0.&&_couplingValueRe==0.)
return false;
if(ids[0]->iCharge()!=ids[1]->iCharge()+ids[2]->iCharge())
return false;
- bool isGVB = (abs(ids[0]->id())==ParticleID::Wplus
- || abs(ids[0]->id())==ParticleID::Z0
- || abs(ids[0]->id())==ParticleID::gamma);
- if(isGVB && ids[1]->iSpin()==PDT::Spin0 && ids[2]->iSpin()==PDT::Spin0)
+ if(ids[0]->iSpin()==PDT::Spin1 && ids[1]->iSpin()==PDT::Spin0 && ids[2]->iSpin()==PDT::Spin0)
return true;
else
return false;
}
vector<pair<int, Complex> >
OneZeroZeroEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
vector<pair<int, Complex> >
OneZeroZeroEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
DecayMEPtr OneZeroZeroEWSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool) {
Complex gvhh(0.,0.);
getCouplings(gvhh,ids);
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin0,PDT::Spin0)));
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
double sqrtmass = sqrt(m0t2*z*(1.-z)-m1t2*(1.-z)-m2t2*z+z*(1.-z));
// assign kernel
(*kernal)(0,0,0) = -cphase*sqrtmass; // 111
(*kernal)(1,0,0) = sqrt(m0t2)*(1.-2.*z)/sqrt(2.); // 211 -> 411
(*kernal)(2,0,0) = phase*sqrtmass; // 311
// return the answer
return kernal;
}
diff --git a/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/ZeroZeroOneEWSplitFn.cc
@@ -1,219 +1,213 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ZeroZeroOneEWSplitFn class.
//
#include "ZeroZeroOneEWSplitFn.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/ParticleData.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "Herwig/Models/StandardModel/SMFFHVertex.h"
#include "ThePEG/Interface/Parameter.h"
using namespace Herwig;
IBPtr ZeroZeroOneEWSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr ZeroZeroOneEWSplitFn::fullclone() const {
return new_ptr(*this);
}
void ZeroZeroOneEWSplitFn::persistentOutput(PersistentOStream & os) const {
os << _couplingValueIm << _couplingValueRe;
}
void ZeroZeroOneEWSplitFn::persistentInput(PersistentIStream & is, int) {
is >> _couplingValueIm >> _couplingValueRe;
}
// The following static variable is needed for the type description system in ThePEG.
DescribeClass<ZeroZeroOneEWSplitFn,SplittingFunction>
describeHerwigZeroZeroOneEWSplitFn("Herwig::ZeroZeroOneEWSplitFn", "HwShower.so");
void ZeroZeroOneEWSplitFn::Init() {
static ClassDocumentation<ZeroZeroOneEWSplitFn> documentation
("The ZeroZeroOneEWSplitFn class implements purly beyond SM electroweak splittings H->H'V");
static Parameter<ZeroZeroOneEWSplitFn, double> interfaceCouplingValueIm
("CouplingValue.Im",
"The numerical value (imaginary part) of the splitting coupling to be imported for BSM splittings",
&ZeroZeroOneEWSplitFn::_couplingValueIm, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
static Parameter<ZeroZeroOneEWSplitFn, double> interfaceCouplingValueRe
("CouplingValue.Re",
"The numerical value (real part) of the splitting coupling to be imported for BSM splittings",
&ZeroZeroOneEWSplitFn::_couplingValueRe, 0.0, -1.0E6, +1.0E6,
false, false, Interface::limited);
}
void ZeroZeroOneEWSplitFn::doinit() {
SplittingFunction::doinit();
}
void ZeroZeroOneEWSplitFn::getCouplings(Complex & g, const IdList &) const {
g = Complex(_couplingValueRe,_couplingValueIm);
}
double ZeroZeroOneEWSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
double rho00 = abs(rho(0,0));
double rho11 = abs(rho(1,1));
double rho22 = abs(rho(2,2));
// the splitting in the massless limit
double val = z*(rho00+rho22)/(1.-z);
// the massive limit
if(mass){
// get the running mass
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
val += (m2t2*sqr(1.+z)*rho11)/(2.*sqr(1.-z))+((-(m1t2*(1.-z))
-m2t2*z+m0t2*(1.-z)*z)*(rho00+rho22))/sqr(1.-z);
}
return norm(ghhv)*val;
}
double ZeroZeroOneEWSplitFn::overestimateP(const double z,
const IdList & ids) const {
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
return norm(ghhv)*z/(1.-z);
}
double ZeroZeroOneEWSplitFn::ratioP(const double z, const Energy2 t,
const IdList & ids, const bool mass,
const RhoDMatrix & rho) const {
double rho00 = abs(rho(0,0));
double rho11 = abs(rho(1,1));
double rho22 = abs(rho(2,2));
// ratio in the massless limit
double val = rho00+rho22;
// the massive limit
if(mass){
// get the running mass
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
val += ((1.-z)*(m2t2*sqr(1.+z)*rho11+2*(m1t2*(-1.+z)
-(m2t2+m0t2*(-1.+z))*z)*(rho00+rho22)))/(2.*sqr(-1.+z)*z);
}
return val;
}
double ZeroZeroOneEWSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
double pre = norm(ghhv);
switch (PDFfactor) {
case 0:
return -pre*(z+log(1.-z));
case 1:
case 2:
case 3:
default:
throw Exception() << "ZeroZeroOneEWSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double ZeroZeroOneEWSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
double pre = norm(ghhv);
switch (PDFfactor) {
case 0:
return 1.-exp(-(1+r/pre));;
case 1:
case 2:
case 3:
default:
throw Exception() << "ZeroZeroOneEWSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool ZeroZeroOneEWSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3)
return false;
if(_couplingValueIm==0.&&_couplingValueRe==0.)
return false;
if(ids[0]->iCharge()!=ids[1]->iCharge()+ids[2]->iCharge())
return false;
- bool is1GVB = (abs(ids[1]->id())==ParticleID::Wplus
- || abs(ids[1]->id())==ParticleID::Z0
- || abs(ids[1]->id())==ParticleID::gamma);
- bool is2GVB = (abs(ids[2]->id())==ParticleID::Wplus
- || abs(ids[2]->id())==ParticleID::Z0
- || abs(ids[2]->id())==ParticleID::gamma);
- if(ids[0]->iSpin()==PDT::Spin0 && is1GVB && ids[2]->iSpin()==PDT::Spin0)
+ if(ids[0]->iSpin()==PDT::Spin0 && ids[1]->iSpin()==PDT::Spin1 && ids[2]->iSpin()==PDT::Spin0)
return true;
- else if(ids[0]->iSpin()==PDT::Spin0 && ids[1]->iSpin()==PDT::Spin0 && is2GVB)
+ else if(ids[0]->iSpin()==PDT::Spin0 && ids[1]->iSpin()==PDT::Spin0 && ids[2]->iSpin()==PDT::Spin1)
return true;
else
return false;
}
vector<pair<int, Complex> >
ZeroZeroOneEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
vector<pair<int, Complex> >
ZeroZeroOneEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & ,
const RhoDMatrix &) {
// no dependence on the spin density matrix, dependence on off-diagonal terms cancels
// and rest = splitting function for Tr(rho)=1 as required by defn
return vector<pair<int, Complex> >(1,make_pair(0,1.));
}
DecayMEPtr ZeroZeroOneEWSplitFn::matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool) {
Complex ghhv(0.,0.);
getCouplings(ghhv,ids);
// calculate the kernal
DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin0,PDT::Spin0)));
double m0t2 = sqr(getParticleData(ids[0]->id())->mass())/t;
double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
Complex phase = exp(Complex(0.,1.)*phi);
Complex cphase = conj(phase);
double sqrtmass = sqrt(m0t2*z*(1.-z)-m1t2*(1.-z)-m2t2*z+z*(1.-z));
// assign kernel
(*kernal)(0,0,0) = -phase*ghhv*sqrtmass/(1.-z); // 111
(*kernal)(1,0,0) = -sqrt(m2t2)*(1.+z)/sqrt(2.*(1.-z)); // 211 -> 411
(*kernal)(2,0,0) = cphase*ghhv*sqrtmass/(1.-z); // 311
// return the answer
return kernal;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:00 PM (18 h, 35 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4020363
Default Alt Text
(34 KB)

Event Timeline