diff --git a/Shower/QTilde/SplittingFunctions/OneOneOneQEDSplitFn.cc b/Shower/QTilde/SplittingFunctions/OneOneOneQEDSplitFn.cc
--- a/Shower/QTilde/SplittingFunctions/OneOneOneQEDSplitFn.cc
+++ b/Shower/QTilde/SplittingFunctions/OneOneOneQEDSplitFn.cc
@@ -1,236 +1,237 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the OneOneOneQEDSplitFn class.
//
#include "OneOneOneQEDSplitFn.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"
using namespace Herwig;
IBPtr OneOneOneQEDSplitFn::clone() const {
return new_ptr(*this);
}
IBPtr OneOneOneQEDSplitFn::fullclone() const {
return new_ptr(*this);
}
void OneOneOneQEDSplitFn::persistentOutput(PersistentOStream & os) const {
os << gWWG_ << _theSM;
}
void OneOneOneQEDSplitFn::persistentInput(PersistentIStream & is, int) {
is >> gWWG_ >> _theSM;
}
// The following static variable is needed for the type description system in ThePEG.
DescribeClass<OneOneOneQEDSplitFn,SplittingFunction>
describeHerwigOneOneOneQEDSplitFn("Herwig::OneOneOneQEDSplitFn", "HwShower.so");
void OneOneOneQEDSplitFn::Init() {
static ClassDocumentation<OneOneOneQEDSplitFn> documentation
("The OneOneOneQEDSplitFn class implements the gamma->WW EW splitting.");
}
void OneOneOneQEDSplitFn::doinit() {
SplittingFunction::doinit();
tcSMPtr sm = generator()->standardModel();
double sw2 = sm->sin2ThetaW();
// WWG coupling
gWWG_ = 1.;
// to employ running masses, wherever needed
_theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
}
void OneOneOneQEDSplitFn::getCouplings(double & gvvv, const IdList & ids) const {
// G > WW
if(ids[0]->id()==ParticleID::gamma && abs(ids[1]->id())==ParticleID::Wplus
&& abs(ids[2]->id())==ParticleID::Wplus){
gvvv = gWWG_;
}
else
assert(false);
}
double OneOneOneQEDSplitFn::P(const double z, const Energy2 t,
const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
double gvvv(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(getParticleData(ids[0]->id())->mass())/t;
- double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
- double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
+ double m0t2 = sqr(ids[0]->mass())/t;
+ double m1t2 = sqr(ids[1]->mass())/t;
+ double m2t2 = sqr(ids[2]->mass())/t;
val += (-2.*(m2t2*(1.-sqr(1.-z)*z)+m1t2*(1.-(1.-z)*sqr(z)))*(abs_rho_00+abs_rho_22))/((1.-z)*z)
+ (2.*m0t2*(2.*pow(1.-z,3)*z*abs_rho_11+sqr(1.-(1.-z)*z)*(abs_rho_00+abs_rho_22)))/((1.-z)*z);
}
return sqr(gvvv)*val;
}
double OneOneOneQEDSplitFn::overestimateP(const double z,
const IdList & ids) const {
double gvvv(0.);
getCouplings(gvvv,ids);
return sqr(gvvv)*(2./(z*(1.-z)));
}
double OneOneOneQEDSplitFn::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(getParticleData(ids[0]->id())->mass())/t;
- double m1t2 = sqr(getParticleData(ids[1]->id())->mass())/t;
- double m2t2 = sqr(getParticleData(ids[2]->id())->mass())/t;
+ cerr<<getParticleData(ids[1]->id())->mass()/GeV<<" "<<ids[1]->mass()/GeV<<"\n";
+ double m0t2 = sqr(ids[0]->mass())/t;
+ double m1t2 = sqr(ids[1]->mass())/t;
+ double m2t2 = sqr(ids[2]->mass())/t;
val += -(m2t2*(1.-sqr(1.-z)*z) + m1t2*(1.-(1.-z)*sqr(z)))*(abs_rho_00+abs_rho_22)
+ m0t2*(2.*pow(1.-z,3)*z*abs_rho_11+sqr(1.-(1.-z)*z)*(abs_rho_00+abs_rho_22));
}
return val;
}
double OneOneOneQEDSplitFn::integOverP(const double z,
const IdList & ids,
unsigned int PDFfactor) const {
double gvvv(0.);
getCouplings(gvvv,ids);
double pre = sqr(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() << "OneOneOneQEDSplitFn::integOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
double OneOneOneQEDSplitFn::invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor) const {
double gvvv(0.);
getCouplings(gvvv,ids);
double pre = sqr(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() << "OneOneOneQEDSplitFn::invIntegOverP() invalid PDFfactor = "
<< PDFfactor << Exception::runerror;
}
}
bool OneOneOneQEDSplitFn::accept(const IdList &ids) const {
if(ids.size()!=3) return false;
if(ids[0]->id()==ParticleID::gamma && abs(ids[1]->id())==ParticleID::Wplus
&& ids[1]->id()==-ids[2]->id())
return true;
return false;
}
vector<pair<int, Complex> >
OneOneOneQEDSplitFn::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> >
OneOneOneQEDSplitFn::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 OneOneOneQEDSplitFn::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)));
double 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;
}