Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7879786
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
36 KB
Subscribers
None
View Options
diff --git a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbar.cc b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbar.cc
--- a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbar.cc
+++ b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbar.cc
@@ -1,184 +1,192 @@
// -*- C++ -*-
//
// MatchboxAmplitudelnuqqbar.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 MatchboxAmplitudelnuqqbar class.
//
#include "Herwig/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbar.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
using namespace Herwig;
MatchboxAmplitudelnuqqbar::MatchboxAmplitudelnuqqbar()
: theDiagonal(false) {}
MatchboxAmplitudelnuqqbar::~MatchboxAmplitudelnuqqbar() {}
void MatchboxAmplitudelnuqqbar::doinit() {
MatchboxAmplitude::doinit();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
CA = SM().Nc();
CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
theCKM = standardCKM(SM())->getUnsquaredMatrix(6);
nPoints(4);
}
void MatchboxAmplitudelnuqqbar::doinitrun() {
MatchboxAmplitude::doinitrun();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
+ MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
+ GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
+ CA = SM().Nc();
+ CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
nPoints(4);
}
IBPtr MatchboxAmplitudelnuqqbar::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxAmplitudelnuqqbar::fullclone() const {
return new_ptr(*this);
}
bool MatchboxAmplitudelnuqqbar::canHandle(const PDVector& proc) const {
PDVector xproc = proc;
if ( xproc[0]->CC() )
xproc[0] = xproc[0]->CC();
if ( xproc[1]->CC() )
xproc[1] = xproc[1]->CC();
PDVector::iterator elektron = xproc.begin();
for ( ; elektron != xproc.end(); ++elektron )
if ( abs((*elektron)->id()) >= 11 && abs((*elektron)->id()) <= 16 && abs((*elektron)->id()) % 2 == 1 ) break;
if ( elektron == xproc.end() ) return false;
PDPtr e = *elektron;
xproc.erase(elektron);
PDVector::iterator neutrino = xproc.begin();
for ( ; neutrino != xproc.end(); ++neutrino )
if ( abs((*neutrino)->id()) >= 11 && abs((*neutrino)->id()) <= 16 && abs((*neutrino)->id()) % 2 == 0 ) break;
if ( neutrino == xproc.end() ) return false;
PDPtr n = *neutrino;
xproc.erase(neutrino);
PDVector::iterator quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 && abs((*quark)->id()) % 2 == 1 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr d = *quark;
xproc.erase(quark);
quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 && abs((*quark)->id()) % 2 == 0 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr u = *quark;
xproc.erase(quark);
if ( u->iCharge() + d->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(u) != SU2Helper::family(d) ) return false;
if ( SU2Helper::family(e) != SU2Helper::family(n) ) return false;
return xproc.empty();
}
void MatchboxAmplitudelnuqqbar::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
amplitudeScale(sqrt(lastSHat()));
setupLeptons(0,amplitudeMomentum(0),1,amplitudeMomentum(1));
momentum(2,amplitudeMomentum(2));
momentum(3,amplitudeMomentum(3));
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex MatchboxAmplitudelnuqqbar::evaluate(size_t, const vector<int>& hel, Complex& largeN) {
if ( abs(hel[2]+hel[3]) != 2 ) {
largeN = 0.;
return 0.;
}
Complex ckmelement = 1.;
if ( !theDiagonal ) {
bool wPlus = ( abs(amplitudePartonData()[0]->id()) % 2 == 0 ) ?
amplitudePartonData()[1]->id() < 0:
amplitudePartonData()[0]->id() < 0;
pair<int,int> tmp(
SU2Helper::family(amplitudePartonData()[2])-1,
SU2Helper::family(amplitudePartonData()[3])-1);
if ( amplitudePartonData()[3]->id() < 0 ) swap(tmp.first,tmp.second);
ckmelement = theCKM[tmp.first][tmp.second];
if ( !wPlus ) ckmelement = conj(ckmelement);
}
Complex wPropergator =
1./Complex(((amplitudeMomentum(0)+amplitudeMomentum(1)).m2()-sqr(MW))/lastSHat(),MW*GW/lastSHat());
Complex wVertices =
2.*SM().alphaEMMZ()*Constants::pi/SM().sin2ThetaW()*ckmelement;
const LorentzVector<Complex>& leptonCurrent = llbarLeftCurrent(0,hel[0],1,hel[1]);
const LorentzVector<Complex>& quarkCurrent = qqbarLeftCurrent(2,hel[2],3,hel[3]);
Complex current = hel[2] == 1 ? Complex(0.,-1)*leptonCurrent.dot(quarkCurrent): 0.;
Complex res = current*wVertices*wPropergator;
largeN = res;
return res;
}
Complex MatchboxAmplitudelnuqqbar::evaluateOneLoop(size_t, const vector<int>& hel) {
if ( abs(hel[2]+hel[3]) != 2 ) return 0.;
Complex ckmelement = 1.;
if ( !theDiagonal ) {
bool wPlus = ( abs(amplitudePartonData()[0]->id()) % 2 == 0 ) ?
amplitudePartonData()[1]->id() < 0:
amplitudePartonData()[0]->id() < 0;
pair<int,int> tmp(
SU2Helper::family(amplitudePartonData()[2])-1,
SU2Helper::family(amplitudePartonData()[3])-1);
if ( amplitudePartonData()[3]->id() < 0 ) swap(tmp.first,tmp.second);
ckmelement = theCKM[tmp.first][tmp.second];
if ( !wPlus ) ckmelement = conj(ckmelement);
}
Complex wPropergator =
1./Complex(((amplitudeMomentum(0)+amplitudeMomentum(1)).m2()-sqr(MW))/lastSHat(),MW*GW/lastSHat());
Complex wVertices =
2.*SM().alphaEMMZ()*Constants::pi/SM().sin2ThetaW()*ckmelement;
const LorentzVector<Complex>& leptonCurrent = llbarLeftCurrent(0,hel[0],1,hel[1]);
const LorentzVector<Complex>& quarkCurrent = qqbarLeftOneLoopCurrent(2,hel[2],3,hel[3]);
Complex current = hel[2] == 1 ? Complex(0.,-1)*leptonCurrent.dot(quarkCurrent): 0.;
Complex res = (SM().alphaS()/(2.*Constants::pi))*current*wVertices*wPropergator;
return res;
}
void MatchboxAmplitudelnuqqbar::persistentOutput(PersistentOStream & os) const {
- os << theDiagonal << theCKM << ounit(MW,GeV) << ounit(GW,GeV) << CA << CF;
+ os << theDiagonal << theCKM ;
}
void MatchboxAmplitudelnuqqbar::persistentInput(PersistentIStream & is, int) {
- is >> theDiagonal >> theCKM >> iunit(MW,GeV) >> iunit(GW,GeV) >> CA >> CF;
+ is >> theDiagonal >> theCKM ;
}
DescribeClass<MatchboxAmplitudelnuqqbar,MatchboxAmplitude>
describeHerwigMatchboxAmplitudelnuqqbar("Herwig::MatchboxAmplitudelnuqqbar", "HwMatchboxBuiltin.so");
void MatchboxAmplitudelnuqqbar::Init() {
static ClassDocumentation<MatchboxAmplitudelnuqqbar> documentation
("MatchboxAmplitudelnuqqbar");
static Switch<MatchboxAmplitudelnuqqbar,bool> interfaceDiagonal
("Diagonal",
"Use a diagonal CKM matrix (ignoring the CKM object of the StandardModel).",
&MatchboxAmplitudelnuqqbar::theDiagonal, false, false, false);
static SwitchOption interfaceDiagonalYes
(interfaceDiagonal,
"Yes",
"Use a diagonal CKM matrix.",
true);
static SwitchOption interfaceDiagonalNo
(interfaceDiagonal,
"No",
"Use the CKM object as used by the StandardModel.",
false);
}
diff --git a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarg.cc b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarg.cc
--- a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarg.cc
+++ b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarg.cc
@@ -1,195 +1,203 @@
// -*- C++ -*-
//
// MatchboxAmplitudelnuqqbarg.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 MatchboxAmplitudelnuqqbarg class.
//
#include "Herwig/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarg.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
using namespace Herwig;
MatchboxAmplitudelnuqqbarg::MatchboxAmplitudelnuqqbarg()
: theDiagonal(false) {}
MatchboxAmplitudelnuqqbarg::~MatchboxAmplitudelnuqqbarg() {}
void MatchboxAmplitudelnuqqbarg::doinit() {
MatchboxAmplitude::doinit();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
CA = SM().Nc();
CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
theCKM = standardCKM(SM())->getUnsquaredMatrix(6);
nPoints(5);
}
void MatchboxAmplitudelnuqqbarg::doinitrun() {
MatchboxAmplitude::doinitrun();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
+ MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
+ GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
+ CA = SM().Nc();
+ CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
nPoints(5);
}
IBPtr MatchboxAmplitudelnuqqbarg::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxAmplitudelnuqqbarg::fullclone() const {
return new_ptr(*this);
}
bool MatchboxAmplitudelnuqqbarg::canHandle(const PDVector& proc) const {
PDVector xproc = proc;
if ( xproc[0]->CC() )
xproc[0] = xproc[0]->CC();
if ( xproc[1]->CC() )
xproc[1] = xproc[1]->CC();
PDVector::iterator elektron = xproc.begin();
for ( ; elektron != xproc.end(); ++elektron )
if ( abs((*elektron)->id()) >= 11 && abs((*elektron)->id()) <= 16 && abs((*elektron)->id()) % 2 == 1 ) break;
if ( elektron == xproc.end() ) return false;
PDPtr e = *elektron;
xproc.erase(elektron);
PDVector::iterator neutrino = xproc.begin();
for ( ; neutrino != xproc.end(); ++neutrino )
if ( abs((*neutrino)->id()) >= 11 && abs((*neutrino)->id()) <= 16 && abs((*neutrino)->id()) % 2 == 0 ) break;
if ( neutrino == xproc.end() ) return false;
PDPtr n = *neutrino;
xproc.erase(neutrino);
PDVector::iterator quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 && abs((*quark)->id()) % 2 == 1 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr d = *quark;
xproc.erase(quark);
quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 && abs((*quark)->id()) % 2 == 0 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr u = *quark;
xproc.erase(quark);
PDVector::iterator gluon = xproc.begin();
for ( ; gluon != xproc.end(); ++gluon )
if ( (*gluon)->id() == 21 ) break;
if ( gluon == xproc.end() ) return false;
xproc.erase(gluon);
if ( SU2Helper::family(e) != SU2Helper::family(n) ) return false;
if ( u->iCharge() + d->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(u) != SU2Helper::family(d) ) return false;
return xproc.empty();
}
void MatchboxAmplitudelnuqqbarg::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
amplitudeScale(sqrt(lastSHat()));
setupLeptons(0,amplitudeMomentum(0),1,amplitudeMomentum(1));
momentum(2,amplitudeMomentum(2));
momentum(3,amplitudeMomentum(3));
momentum(4,amplitudeMomentum(4));
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex MatchboxAmplitudelnuqqbarg::evaluate(size_t, const vector<int>& hel, Complex& largeN) {
if ( abs(hel[2]+hel[3]) != 2 ) {
largeN = 0.;
return 0.;
}
Complex ckmelement = 1.;
if ( !theDiagonal ) {
bool wPlus = ( abs(amplitudePartonData()[0]->id()) % 2 == 0 ) ?
amplitudePartonData()[1]->id() < 0:
amplitudePartonData()[0]->id() < 0;
pair<int,int> tmp(
SU2Helper::family(amplitudePartonData()[2])-1,
SU2Helper::family(amplitudePartonData()[3])-1);
if ( amplitudePartonData()[3]->id() < 0 ) swap(tmp.first,tmp.second);
ckmelement = theCKM[tmp.first][tmp.second];
if ( !wPlus ) ckmelement = conj(ckmelement);
}
Complex wPropergator =
1./Complex(((amplitudeMomentum(0)+amplitudeMomentum(1)).m2()-sqr(MW))/lastSHat(),MW*GW/lastSHat());
Complex wVertices =
2.*SM().alphaEMMZ()*Constants::pi/SM().sin2ThetaW()*ckmelement;
Complex sVertex =
sqrt(4.*Constants::pi*SM().alphaS());
const LorentzVector<Complex>& leptonCurrent = llbarLeftCurrent(0,hel[0],1,hel[1]);
const LorentzVector<Complex>& quarkCurrent = qqbargLeftCurrent(2,hel[2],3,hel[3],4,hel[4]);
Complex current = hel[2] == 1 ? Complex(0.,-1)*leptonCurrent.dot(quarkCurrent): 0.;
Complex res = current*wVertices*wPropergator*sVertex;
largeN = res;
return res;
}
Complex MatchboxAmplitudelnuqqbarg::evaluateOneLoop(size_t, const vector<int>& hel) {
if ( abs(hel[2]+hel[3]) != 2 ) return 0.;
Complex ckmelement = 1.;
if ( !theDiagonal ) {
bool wPlus = ( abs(amplitudePartonData()[0]->id()) % 2 == 0 ) ?
amplitudePartonData()[1]->id() < 0:
amplitudePartonData()[0]->id() < 0;
pair<int,int> tmp(
SU2Helper::family(amplitudePartonData()[2])-1,
SU2Helper::family(amplitudePartonData()[3])-1);
if ( amplitudePartonData()[3]->id() < 0 ) swap(tmp.first,tmp.second);
ckmelement = theCKM[tmp.first][tmp.second];
if ( !wPlus ) ckmelement = conj(ckmelement);
}
Complex wPropergator =
1./Complex(((amplitudeMomentum(0)+amplitudeMomentum(1)).m2()-sqr(MW))/lastSHat(),MW*GW/lastSHat());
Complex wVertices =
2.*SM().alphaEMMZ()*Constants::pi/SM().sin2ThetaW()*ckmelement;
Complex sVertex =
sqrt(4.*Constants::pi*SM().alphaS());
const LorentzVector<Complex>& leptonCurrent = llbarLeftCurrent(0,hel[0],1,hel[1]);
const LorentzVector<Complex>& quarkCurrent = qqbargLeftOneLoopCurrent(2,hel[2],3,hel[3],4,hel[4]);
Complex current = hel[2] == 1 ? Complex(0.,-1)*leptonCurrent.dot(quarkCurrent): 0.;
Complex res = (SM().alphaS()/(2.*Constants::pi))*current*wVertices*wPropergator*sVertex;
return res;
}
void MatchboxAmplitudelnuqqbarg::persistentOutput(PersistentOStream & os) const {
- os << theDiagonal << theCKM << ounit(MW,GeV) << ounit(GW,GeV) << CA << CF;
+ os << theDiagonal << theCKM ;
}
void MatchboxAmplitudelnuqqbarg::persistentInput(PersistentIStream & is, int) {
- is >> theDiagonal >> theCKM >> iunit(MW,GeV) >> iunit(GW,GeV) >> CA >> CF;
+ is >> theDiagonal >> theCKM ;
}
DescribeClass<MatchboxAmplitudelnuqqbarg,MatchboxAmplitude>
describeHerwigMatchboxAmplitudelnuqqbarg("Herwig::MatchboxAmplitudelnuqqbarg", "HwMatchboxBuiltin.so");
void MatchboxAmplitudelnuqqbarg::Init() {
static ClassDocumentation<MatchboxAmplitudelnuqqbarg> documentation
("MatchboxAmplitudelnuqqbarg");
static Switch<MatchboxAmplitudelnuqqbarg,bool> interfaceDiagonal
("Diagonal",
"Use a diagonal CKM matrix (ignoring the CKM object of the StandardModel).",
&MatchboxAmplitudelnuqqbarg::theDiagonal, false, false, false);
static SwitchOption interfaceDiagonalYes
(interfaceDiagonal,
"Yes",
"Use a diagonal CKM matrix.",
true);
static SwitchOption interfaceDiagonalNo
(interfaceDiagonal,
"No",
"Use the CKM object as used by the StandardModel.",
false);
}
diff --git a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbargg.cc b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbargg.cc
--- a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbargg.cc
+++ b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbargg.cc
@@ -1,193 +1,201 @@
// -*- C++ -*-
//
// MatchboxAmplitudelnuqqbargg.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 MatchboxAmplitudelnuqqbargg class.
//
#include "Herwig/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbargg.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
using namespace Herwig;
MatchboxAmplitudelnuqqbargg::MatchboxAmplitudelnuqqbargg()
: theDiagonal(false) {}
MatchboxAmplitudelnuqqbargg::~MatchboxAmplitudelnuqqbargg() {}
void MatchboxAmplitudelnuqqbargg::doinit() {
MatchboxAmplitude::doinit();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
CA = SM().Nc();
CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
theCKM = standardCKM(SM())->getUnsquaredMatrix(6);
nPoints(6);
}
void MatchboxAmplitudelnuqqbargg::doinitrun() {
MatchboxAmplitude::doinitrun();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
+ MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
+ GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
+ CA = SM().Nc();
+ CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
nPoints(6);
}
IBPtr MatchboxAmplitudelnuqqbargg::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxAmplitudelnuqqbargg::fullclone() const {
return new_ptr(*this);
}
bool MatchboxAmplitudelnuqqbargg::canHandle(const PDVector& proc) const {
PDVector xproc = proc;
if ( xproc[0]->CC() )
xproc[0] = xproc[0]->CC();
if ( xproc[1]->CC() )
xproc[1] = xproc[1]->CC();
PDVector::iterator elektron = xproc.begin();
for ( ; elektron != xproc.end(); ++elektron )
if ( abs((*elektron)->id()) >= 11 && abs((*elektron)->id()) <= 16 && abs((*elektron)->id()) % 2 == 1 ) break;
if ( elektron == xproc.end() ) return false;
PDPtr e = *elektron;
xproc.erase(elektron);
PDVector::iterator neutrino = xproc.begin();
for ( ; neutrino != xproc.end(); ++neutrino )
if ( abs((*neutrino)->id()) >= 11 && abs((*neutrino)->id()) <= 16 && abs((*neutrino)->id()) % 2 == 0 ) break;
if ( neutrino == xproc.end() ) return false;
PDPtr n = *neutrino;
xproc.erase(neutrino);
PDVector::iterator quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 && abs((*quark)->id()) % 2 == 1 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr d = *quark;
xproc.erase(quark);
quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 && abs((*quark)->id()) % 2 == 0 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr u = *quark;
xproc.erase(quark);
PDVector::iterator gluon = xproc.begin();
for ( ; gluon != xproc.end(); ++gluon )
if ( (*gluon)->id() == 21 ) break;
if ( gluon == xproc.end() ) return false;
xproc.erase(gluon);
gluon = xproc.begin();
for ( ; gluon != xproc.end(); ++gluon )
if ( (*gluon)->id() == 21 ) break;
if ( gluon == xproc.end() ) return false;
xproc.erase(gluon);
if ( SU2Helper::family(e) != SU2Helper::family(n) ) return false;
if ( u->iCharge() + d->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(u) != SU2Helper::family(d) ) return false;
return xproc.empty();
}
void MatchboxAmplitudelnuqqbargg::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
amplitudeScale(sqrt(lastSHat()));
setupLeptons(0,amplitudeMomentum(0),1,amplitudeMomentum(1));
momentum(2,amplitudeMomentum(2));
momentum(3,amplitudeMomentum(3));
momentum(4,amplitudeMomentum(4));
momentum(5,amplitudeMomentum(5));
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex MatchboxAmplitudelnuqqbargg::evaluate(size_t a, const vector<int>& hel, Complex& largeN) {
if ( abs(hel[2]+hel[3]) != 2 ) {
largeN = 0.;
return 0.;
}
assert ( amplitudeToColourMap()[2] == 0 && amplitudeToColourMap()[3] == 1 );
int g1,hg1,g2,hg2;
if ( amplitudeToColourMap()[4] == 2 && amplitudeToColourMap()[5] == 3 ) {
if ( a == 0 ) {
g1 = 4; hg1 = hel[4];
g2 = 5; hg2 = hel[5];
} else if ( a == 1 ) {
g1 = 5; hg1 = hel[5];
g2 = 4; hg2 = hel[4];
} else assert ( false );
} else if ( amplitudeToColourMap()[4] == 3 && amplitudeToColourMap()[5] == 2 ) {
if ( a == 0 ) {
g1 = 5; hg1 = hel[5];
g2 = 4; hg2 = hel[4];
} else if ( a == 1 ) {
g1 = 4; hg1 = hel[4];
g2 = 5; hg2 = hel[5];
} else assert ( false );
} else assert ( false );
Complex ckmelement = 1.;
if ( !theDiagonal ) {
bool wPlus = ( abs(amplitudePartonData()[0]->id()) % 2 == 0 ) ?
amplitudePartonData()[1]->id() < 0:
amplitudePartonData()[0]->id() < 0;
pair<int,int> tmp(
SU2Helper::family(amplitudePartonData()[2]),
SU2Helper::family(amplitudePartonData()[3]));
if ( amplitudePartonData()[3]->id() < 0 ) swap(tmp.first,tmp.second);
ckmelement = theCKM[tmp.first][tmp.second];
if ( !wPlus ) ckmelement = conj(ckmelement);
}
Complex wPropergator =
1./Complex(((amplitudeMomentum(0)+amplitudeMomentum(1)).m2()-sqr(MW))/lastSHat(),MW*GW/lastSHat());
Complex wVertices =
2.*SM().alphaEMMZ()*Constants::pi/SM().sin2ThetaW()*ckmelement;
Complex sVertices =
4.*Constants::pi*SM().alphaS();
const LorentzVector<Complex>& leptonCurrent = llbarLeftCurrent(0,hel[0],1,hel[1]);
const LorentzVector<Complex>& quarkCurrent = qqbarggLeftCurrent(2,hel[2],3,hel[3],g1,hg1,g2,hg2);
Complex current = hel[2] == 1 ? Complex(0.,-1)*leptonCurrent.dot(quarkCurrent): 0.;
Complex res = current*wVertices*wPropergator*sVertices;
largeN = res;
return res;
}
void MatchboxAmplitudelnuqqbargg::persistentOutput(PersistentOStream & os) const {
- os << theDiagonal << theCKM << ounit(MW,GeV) << ounit(GW,GeV) << CA << CF;
+ os << theDiagonal << theCKM ;
}
void MatchboxAmplitudelnuqqbargg::persistentInput(PersistentIStream & is, int) {
- is >> theDiagonal >> theCKM >> iunit(MW,GeV) >> iunit(GW,GeV) >> CA >> CF;
+ is >> theDiagonal >> theCKM ;
}
DescribeClass<MatchboxAmplitudelnuqqbargg,MatchboxAmplitude>
describeHerwigMatchboxAmplitudelnuqqbargg("Herwig::MatchboxAmplitudelnuqqbargg", "HwMatchboxBuiltin.so");
void MatchboxAmplitudelnuqqbargg::Init() {
static ClassDocumentation<MatchboxAmplitudelnuqqbargg> documentation
("MatchboxAmplitudelnuqqbargg");
static Switch<MatchboxAmplitudelnuqqbargg,bool> interfaceDiagonal
("Diagonal",
"Use a diagonal CKM matrix (ignoring the CKM object of the StandardModel).",
&MatchboxAmplitudelnuqqbargg::theDiagonal, false, false, false);
static SwitchOption interfaceDiagonalYes
(interfaceDiagonal,
"Yes",
"Use a diagonal CKM matrix.",
true);
static SwitchOption interfaceDiagonalNo
(interfaceDiagonal,
"No",
"Use the CKM object as used by the StandardModel.",
false);
}
diff --git a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarqqbar.cc b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarqqbar.cc
--- a/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarqqbar.cc
+++ b/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarqqbar.cc
@@ -1,267 +1,275 @@
// -*- C++ -*-
//
// MatchboxAmplitudelnuqqbarqqbar.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 MatchboxAmplitudelnuqqbarqqbar class.
//
#include "Herwig/MatrixElement/Matchbox/Builtin/Amplitudes/MatchboxAmplitudelnuqqbarqqbar.h"
using namespace Herwig;
MatchboxAmplitudelnuqqbarqqbar::MatchboxAmplitudelnuqqbarqqbar()
: theDiagonal(false) {}
MatchboxAmplitudelnuqqbarqqbar::~MatchboxAmplitudelnuqqbarqqbar() {}
void MatchboxAmplitudelnuqqbarqqbar::doinit() {
MatchboxAmplitude::doinit();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
CA = SM().Nc();
CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
theCKM = standardCKM(SM())->getUnsquaredMatrix(6);
nPoints(6);
}
void MatchboxAmplitudelnuqqbarqqbar::doinitrun() {
MatchboxAmplitude::doinitrun();
+ MZ = getParticleData(ParticleID::Z0)->hardProcessMass();
+ GZ = getParticleData(ParticleID::Z0)->hardProcessWidth();
+ MW = getParticleData(ParticleID::Wplus)->hardProcessMass();
+ GW = getParticleData(ParticleID::Wplus)->hardProcessWidth();
+ CA = SM().Nc();
+ CF = (SM().Nc()*SM().Nc()-1.)/(2.*SM().Nc());
nPoints(6);
}
IBPtr MatchboxAmplitudelnuqqbarqqbar::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxAmplitudelnuqqbarqqbar::fullclone() const {
return new_ptr(*this);
}
bool MatchboxAmplitudelnuqqbarqqbar::canHandle(const PDVector& proc) const {
PDVector xproc = proc;
if ( xproc[0]->CC() )
xproc[0] = xproc[0]->CC();
if ( xproc[1]->CC() )
xproc[1] = xproc[1]->CC();
// Charge charge(ZERO);
PDVector::iterator elektron = xproc.begin();
for ( ; elektron != xproc.end(); ++elektron )
if ( abs((*elektron)->id()) >= 11 && abs((*elektron)->id()) <= 16 && abs((*elektron)->id()) % 2 == 1 ) break;
if ( elektron == xproc.end() ) return false;
PDPtr e = *elektron;
xproc.erase(elektron);
PDVector::iterator neutrino = xproc.begin();
for ( ; neutrino != xproc.end(); ++neutrino )
if ( abs((*neutrino)->id()) >= 11 && abs((*neutrino)->id()) <= 16 && abs((*neutrino)->id()) % 2 == 0 ) break;
if ( neutrino == xproc.end() ) return false;
PDPtr n = *neutrino;
xproc.erase(neutrino);
PDVector::iterator quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr q1 = *quark;
xproc.erase(quark);
quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr q2 = *quark;
xproc.erase(quark);
quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr q3 = *quark;
xproc.erase(quark);
quark = xproc.begin();
for ( ; quark != xproc.end(); ++quark )
if ( abs((*quark)->id()) >= 1 && abs((*quark)->id()) <= 6 ) {
assert( (*quark)->hardProcessMass() == ZERO );
break;
}
if ( quark == xproc.end() ) return false;
PDPtr q4 = *quark;
xproc.erase(quark);
if ( q1->id() == -q2->id() ) {
if ( q3->iCharge() + q4->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(q3) != SU2Helper::family(q4) ) return false;
} else if ( q1->id() == -q3->id() ) {
if ( q2->iCharge() + q4->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(q2) != SU2Helper::family(q4) ) return false;
} else if ( q1->id() == -q4->id() ) {
if ( q2->iCharge() + q3->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(q2) != SU2Helper::family(q3) ) return false;
} else if ( q2->id() == -q3->id() ) {
if ( q1->iCharge() + q4->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(q1) != SU2Helper::family(q4) ) return false;
} else if ( q2->id() == -q4->id() ) {
if ( q1->iCharge() + q3->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(q1) != SU2Helper::family(q3) ) return false;
} else if ( q3->id() == -q4->id() ) {
if ( q1->iCharge() + q2->iCharge() + e->iCharge() != PDT::ChargeNeutral ) return false;
if ( theDiagonal && SU2Helper::family(q1) != SU2Helper::family(q2) ) return false;
} else return false;
if ( SU2Helper::family(e) != SU2Helper::family(n) ) return false;
return xproc.empty();
}
void MatchboxAmplitudelnuqqbarqqbar::prepareAmplitudes(Ptr<MatchboxMEBase>::tcptr me) {
if ( !calculateTreeAmplitudes() ) {
MatchboxAmplitude::prepareAmplitudes(me);
return;
}
amplitudeScale(sqrt(lastSHat()));
setupLeptons(0,amplitudeMomentum(0),1,amplitudeMomentum(1));
momentum(2,amplitudeMomentum(2));
momentum(3,amplitudeMomentum(3));
momentum(4,amplitudeMomentum(4));
momentum(5,amplitudeMomentum(5));
MatchboxAmplitude::prepareAmplitudes(me);
}
Complex MatchboxAmplitudelnuqqbarqqbar::evaluate(size_t a, const vector<int>& hel, Complex& largeN) {
const LorentzVector<Complex>& leptonCurrent = llbarLeftCurrent(0,hel[0],1,hel[1]);
Complex Current2345 =
hel[2] == 1 && hel[3] == 1 && abs(hel[4]+hel[5]) == 2 &&
abs(amplitudePartonData()[4]->id()) == abs(amplitudePartonData()[5]->id()) ?
leptonCurrent.dot(qqbarqqbarLeftCurrent(2,hel[2],3,hel[3],4,hel[4],5,hel[5])) : 0.;
Complex Current4523 =
hel[4] == 1 && hel[5] == 1 && abs(hel[2]+hel[3]) == 2 &&
abs(amplitudePartonData()[2]->id()) == abs(amplitudePartonData()[3]->id()) ?
leptonCurrent.dot(qqbarqqbarLeftCurrent(4,hel[4],5,hel[5],2,hel[2],3,hel[3])) : 0.;
Complex Current2543 =
hel[2] == 1 && hel[5] == 1 && abs(hel[4]+hel[3]) == 2 &&
abs(amplitudePartonData()[4]->id()) == abs(amplitudePartonData()[3]->id()) ?
-leptonCurrent.dot(qqbarqqbarLeftCurrent(2,hel[2],5,hel[5],4,hel[4],3,hel[3])) : 0.;
Complex Current4325 =
hel[4] == 1 && hel[3] == 1 && abs(hel[2]+hel[5]) == 2 &&
abs(amplitudePartonData()[2]->id()) == abs(amplitudePartonData()[5]->id()) ?
-leptonCurrent.dot(qqbarqqbarLeftCurrent(4,hel[4],3,hel[3],2,hel[2],5,hel[5])) : 0.;
Complex ckmelement23 = 1.;
Complex ckmelement45 = 1.;
if ( !theDiagonal ) {
bool wPlus = ( abs(amplitudePartonData()[0]->id()) % 2 == 0 ) ?
amplitudePartonData()[1]->id() < 0:
amplitudePartonData()[0]->id() < 0;
pair<int,int> tmp23(
SU2Helper::family(amplitudePartonData()[2])-1,
SU2Helper::family(amplitudePartonData()[3])-1);
pair<int,int> tmp45(
SU2Helper::family(amplitudePartonData()[4])-1,
SU2Helper::family(amplitudePartonData()[5])-1);
if ( amplitudePartonData()[3]->id() < 0 ) swap(tmp23.first,tmp23.second);
if ( amplitudePartonData()[5]->id() < 0 ) swap(tmp45.first,tmp45.second);
ckmelement23 = theCKM[tmp23.first][tmp23.second];
ckmelement45 = theCKM[tmp45.first][tmp45.second];
if ( !wPlus ) {
ckmelement23 = conj(ckmelement23);
ckmelement45 = conj(ckmelement45);
}
}
Complex wPropergator =
1./Complex(((amplitudeMomentum(0)+amplitudeMomentum(1)).m2()-sqr(MW))/lastSHat(),MW*GW/lastSHat());
Complex wVertices23 =
2.*SM().alphaEMMZ()*Constants::pi/SM().sin2ThetaW()*ckmelement23;
Complex wVertices45 =
2.*SM().alphaEMMZ()*Constants::pi/SM().sin2ThetaW()*ckmelement45;
Complex sVertices =
4.*Constants::pi*SM().alphaS();
Complex res2345 =
Complex(0.,-1.)*wPropergator*sVertices*Current2345*wVertices23;
Complex res2543 =
Complex(0.,-1.)*wPropergator*sVertices*Current2543*wVertices23;
Complex res4523 =
Complex(0.,-1.)*wPropergator*sVertices*Current4523*wVertices45;
Complex res4325 =
Complex(0.,-1.)*wPropergator*sVertices*Current4325*wVertices45;
double Nc = SM().Nc();
Complex resLeading = 0.;
Complex resSubLeading = 0.;
if ( amplitudeToColourMap()[2] == 0 && amplitudeToColourMap()[3] == 1 &&
amplitudeToColourMap()[4] == 2 && amplitudeToColourMap()[5] == 3 ) {
if ( a == 0 ) { //(23)(45)
resLeading = res2543 + res4325;
resSubLeading = res2345 + res4523;
} else if ( a == 1 ) { //(25)(43)
resLeading = res2345 + res4523;
resSubLeading = res2543 + res4325;
} else assert(false);
} else if ( amplitudeToColourMap()[2] == 0 && amplitudeToColourMap()[3] == 3 &&
amplitudeToColourMap()[4] == 2 && amplitudeToColourMap()[5] == 1 ) {
if ( a == 0 ) { // (25)(43)
resLeading = res2345 + res4523;
resSubLeading = res2543 + res4325;
} else if ( a == 1 ) { // (23)(45)
resLeading = res2543 + res4325;
resSubLeading = res2345 + res4523;
} else assert(false);
} else if ( amplitudeToColourMap()[2] == 2 && amplitudeToColourMap()[3] == 3 &&
amplitudeToColourMap()[4] == 0 && amplitudeToColourMap()[5] == 1 ) {
if ( a == 0 ) { //(23)(45)
resLeading = res2543 + res4325;
resSubLeading = res2345 + res4523;
} else if ( a == 1 ) { //(25)(43)
resLeading = res2345 + res4523;
resSubLeading = res2543 + res4325;
} else assert(false);
} else if ( amplitudeToColourMap()[2] == 2 && amplitudeToColourMap()[3] == 1 &&
amplitudeToColourMap()[4] == 0 && amplitudeToColourMap()[5] == 3 ) {
if ( a == 0 ) { //(25)(43)
resLeading = res2345 + res4523;
resSubLeading = res2543 + res4325;
} else if ( a == 1 ) { //(23)(45)
resLeading = res2543 + res4325;
resSubLeading = res2345 + res4523;
} else assert(false);
} else assert(false);
resSubLeading *= -1./Nc;
largeN = resLeading/2.;
return (resLeading + resSubLeading)/2.;
}
void MatchboxAmplitudelnuqqbarqqbar::persistentOutput(PersistentOStream & os) const {
- os << theDiagonal << theCKM << ounit(MW,GeV) << ounit(GW,GeV) << CA << CF;
+ os << theDiagonal << theCKM ;
}
void MatchboxAmplitudelnuqqbarqqbar::persistentInput(PersistentIStream & is, int) {
- is >> theDiagonal >> theCKM >> iunit(MW,GeV) >> iunit(GW,GeV) >> CA >> CF;
+ is >> theDiagonal >> theCKM ;
}
DescribeClass<MatchboxAmplitudelnuqqbarqqbar,MatchboxAmplitude>
describeHerwigMatchboxAmplitudelnuqqbarqqbar("Herwig::MatchboxAmplitudelnuqqbarqqbar", "HwMatchboxBuiltin.so");
void MatchboxAmplitudelnuqqbarqqbar::Init() {
static ClassDocumentation<MatchboxAmplitudelnuqqbarqqbar> documentation
("MatchboxAmplitudelnuqqbarqqbar");
static Switch<MatchboxAmplitudelnuqqbarqqbar,bool> interfaceDiagonal
("Diagonal",
"Use a diagonal CKM matrix (ignoring the CKM object of the StandardModel).",
&MatchboxAmplitudelnuqqbarqqbar::theDiagonal, false, false, false);
static SwitchOption interfaceDiagonalYes
(interfaceDiagonal,
"Yes",
"Use a diagonal CKM matrix.",
true);
static SwitchOption interfaceDiagonalNo
(interfaceDiagonal,
"No",
"Use the CKM object as used by the StandardModel.",
false);
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 9:00 PM (1 d, 40 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3801875
Default Alt Text
(36 KB)
Attached To
R563 testingHerwigHG
Event Timeline
Log In to Comment