No OneTemporary
Actions

Size

29 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Lepton/MEee2Higgs2SM.cc b/MatrixElement/Lepton/MEee2Higgs2SM.cc
	--- a/MatrixElement/Lepton/MEee2Higgs2SM.cc
	+++ b/MatrixElement/Lepton/MEee2Higgs2SM.cc
	@@ -1,354 +1,422 @@
	// -- C++ --
	//
	// MEee2Higgs2SM.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 MEee2Higgs2SM class.
	//

	#include "MEee2Higgs2SM.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Handlers/StandardXComb.h"
	#include "Herwig/MatrixElement/HardVertex.h"
	#include "ThePEG/PDT/EnumParticles.h"

	using namespace Herwig;

	void MEee2Higgs2SM::doinit() {
	ME2to2Base::doinit();
	h0_ = getParticleData(ThePEG::ParticleID::h0);
	tcHwSMPtr hwsm=ThePEG::dynamic_ptr_cast<tcHwSMPtr>(standardModel());
	// do the initialisation
	if(hwsm) {
	FFHVertex_ = hwsm->vertexFFH();
	HGGVertex_ = hwsm->vertexHGG();
	+ HWWVertex_ = hwsm->vertexWWH();
	}
	else {
	throw InitException() << "Must have Herwig StandardModel object in "
	<< "MEee2Higgs2SM::doinit()" << Exception::runerror;
	}
	}

	Energy2 MEee2Higgs2SM::scale() const {
	return sHat();
	}

	unsigned int MEee2Higgs2SM::orderInAlphaS() const {
	return 0;
	}

	unsigned int MEee2Higgs2SM::orderInAlphaEW() const {
	return 2;
	}

	void MEee2Higgs2SM::getDiagrams() const {
	// specific the diagrams
	tcPDPtr ep = getParticleData(ParticleID::eplus );
	tcPDPtr em = getParticleData(ParticleID::eminus);
	// outgoing quarks
	for(int i=1;i<=5;++i) {
	if(allowed_==0 \|\| allowed_==1 \|\| allowed_==i+2) {
	tcPDPtr lm = getParticleData(i);
	tcPDPtr lp = lm->CC();
	add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, lm, 3, lp, -1)));
	}
	}
	// outgoing leptons
	for( int i =11;i<=16;i+=2) {
	if(allowed_==0 \|\| allowed_==2 \|\| allowed_==(i+7)/2) {
	tcPDPtr lm = getParticleData(i);
	tcPDPtr lp = lm->CC();
	add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, lm, 3, lp, -1)));
	}
	}
	if(allowed_==0 \|\| allowed_==12) {
	tcPDPtr g = getParticleData(ParticleID::g);
	add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, g, 3, g, -1)));
	}
	+ if(allowed_==0 \|\| allowed_==13) {
	+ tcPDPtr Wp = getParticleData(ParticleID::Wplus);
	+ add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, Wp, 3, Wp->CC(), -1)));
	+ }
	+ if(allowed_==0 \|\| allowed_==14) {
	+ tcPDPtr Z0 = getParticleData(ParticleID::Z0);
	+ add(new_ptr((Tree2toNDiagram(2), em, ep, 1, h0_, 3, Z0, 3, Z0, -1)));
	+ }
	}

	double MEee2Higgs2SM::me2() const {
	double aver=0.;
	// get the order right
	int ielectron(0),ipositron(1),ilp(2),ilm(3);
	if(mePartonData()[0]->id()!=11) swap(ielectron,ipositron);
	if(mePartonData()[2]->id()>mePartonData()[3]->id()) swap(ilp,ilm);
	// the arrays for the wavefunction to be passed to the matrix element
	vector<SpinorWaveFunction> fin;
	vector<SpinorBarWaveFunction> ain;
	for(unsigned int ihel=0;ihel<2;++ihel) {
	fin .push_back(SpinorWaveFunction(meMomenta()[ielectron],
	mePartonData()[ielectron],ihel,incoming));
	ain .push_back(SpinorBarWaveFunction(meMomenta()[ipositron],
	mePartonData()[ipositron],ihel,incoming));
	}
	// H -> f fbar
	- if(mePartonData()[2]->id()!=ParticleID::g) {
	+ if(abs(mePartonData()[2]->id())<20) {
	vector<SpinorWaveFunction> aout;
	vector<SpinorBarWaveFunction> fout;
	for(unsigned int ihel=0;ihel<2;++ihel) {
	fout.push_back(SpinorBarWaveFunction(meMomenta()[ilm],
	mePartonData()[ilm],ihel,outgoing));
	aout.push_back(SpinorWaveFunction(meMomenta()[ilp],
	mePartonData()[ilp],ihel,outgoing));
	}
	HelicityME(fin,ain,fout,aout,aver);
	if(mePartonData()[ilm]->id()<=6) aver*=3.;
	}
	- else {
	+ else if(mePartonData()[2]->id()==21) {
	vector<VectorWaveFunction> g1,g2;
	for(unsigned int ihel=0;ihel<2;++ihel) {
	g1.push_back(VectorWaveFunction(meMomenta()[2],mePartonData()[2],
	2*ihel,outgoing));
	g2.push_back(VectorWaveFunction(meMomenta()[3],mePartonData()[3],
	2*ihel,outgoing));
	}
	ggME(fin,ain,g1,g2,aver);
	aver *= 8.;
	}
	+ else if(abs(mePartonData()[2]->id())==24 \|\|mePartonData()[2]->id()==23) {
	+ vector<VectorWaveFunction> g1,g2;
	+ for(unsigned int ihel=0;ihel<3;++ihel) {
	+ g1.push_back(VectorWaveFunction(meMomenta()[2],mePartonData()[2],
	+ ihel,outgoing));
	+ g2.push_back(VectorWaveFunction(meMomenta()[3],mePartonData()[3],
	+ ihel,outgoing));
	+ }
	+ WWME(fin,ain,g1,g2,aver);
	+ }
	+ else {
	+ assert(false);
	+ }
	return aver;
	}

	// the helicity amplitude matrix element
	ProductionMatrixElement MEee2Higgs2SM::HelicityME(vector<SpinorWaveFunction> fin,
	vector<SpinorBarWaveFunction> ain,
	vector<SpinorBarWaveFunction> fout,
	vector<SpinorWaveFunction> aout,
	double & aver) const {
	// the particles should be in the order
	// for the incoming
	// 0 incoming fermion (u spinor)
	// 1 incoming antifermion (vbar spinor)
	// for the outgoing
	// 0 outgoing fermion (ubar spinor)
	// 1 outgoing antifermion (v spinor)
	// me to be returned
	ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half,
	PDT::Spin1Half,PDT::Spin1Half);
	// wavefunctions for the intermediate particles
	ScalarWaveFunction interh;
	// temporary storage of the different diagrams
	Complex diag;
	aver=0.;
	// sum over helicities to get the matrix element
	unsigned int inhel1,inhel2,outhel1,outhel2;
	for(inhel1=0;inhel1<2;++inhel1) {
	for(inhel2=0;inhel2<2;++inhel2) {
	interh = FFHVertex_->evaluate(sHat(),1,h0_,fin[inhel1],ain[inhel2]);
	for(outhel1=0;outhel1<2;++outhel1) {
	for(outhel2=0;outhel2<2;++outhel2) {
	diag = FFHVertex_->evaluate(sHat(),aout[outhel2],
	fout[outhel1],interh);
	output(inhel1,inhel2,outhel1,outhel2)=diag;
	aver +=real(diag*conj(diag));
	}
	}
	}
	}
	return output;
	}

	Selector<MEBase::DiagramIndex>
	MEee2Higgs2SM::diagrams(const DiagramVector &) const {
	Selector<DiagramIndex> sel;sel.insert(1.0, 0);
	return sel;
	}
	Selector<const ColourLines *>
	MEee2Higgs2SM::colourGeometries(tcDiagPtr diag) const {
	static ColourLines neutral ( " " );
	static ColourLines quarks ( "-5 4");
	static ColourLines gluons ( "-5 4, 5 -4");
	Selector<const ColourLines *> sel;
	int id = abs((diag->partons()[2])->id());
	if (id<=6 )
	sel.insert(1.0, &quarks);
	else if (id==21)
	sel.insert(1.0, &gluons);
	else
	sel.insert(1.0, &neutral);
	return sel;
	}

	void MEee2Higgs2SM::persistentOutput(PersistentOStream & os) const {
	os << FFHVertex_ << HGGVertex_ << h0_ << allowed_;
	}

	void MEee2Higgs2SM::persistentInput(PersistentIStream & is, int) {
	is >> FFHVertex_ >> HGGVertex_ >> h0_ >> allowed_;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<MEee2Higgs2SM,ME2to2Base>
	describeHerwigMEee2Higgs2SM("Herwig::MEee2Higgs2SM", "HwMELepton.so");

	void MEee2Higgs2SM::Init() {

	static ClassDocumentation<MEee2Higgs2SM> documentation
	("The MEee2Higgs2SM class implements the matrix element for e+e- to"
	" SM particle via Higgs exchnage and is designed to be a process to test"
	" things involving scalar particles. ");

	static Switch<MEee2Higgs2SM,int> interfaceallowed
	("Allowed",
	"Allowed outgoing particles",
	&MEee2Higgs2SM::allowed_, 0, false, false);
	static SwitchOption interfaceallowedAll
	(interfaceallowed,
	"All",
	"All SM particles allowed",
	0);
	static SwitchOption interfaceallowed1
	(interfaceallowed,
	"Quarks",
	"Only the quarks allowed",
	1);
	static SwitchOption interfaceallowed2
	(interfaceallowed,
	"Leptons",
	"Only the leptons allowed",
	2);
	static SwitchOption interfacealloweddown
	(interfaceallowed,
	"Down",
	"Only d dbar allowed",
	3);
	static SwitchOption interfaceallowedup
	(interfaceallowed,
	"Up",
	"Only u ubar allowed",
	4);
	static SwitchOption interfaceallowedstrange
	(interfaceallowed,
	"Strange",
	"Only s sbar allowed",
	5);
	static SwitchOption interfaceallowedcharm
	(interfaceallowed,
	"Charm",
	"Only c cbar allowed",
	6);
	static SwitchOption interfaceallowedbottom
	(interfaceallowed,
	"Bottom",
	"Only b bbar allowed",
	7);
	static SwitchOption interfaceallowedtop
	(interfaceallowed,
	"Top",
	"Only t tbar allowed",
	8);
	static SwitchOption interfaceallowedelectron
	(interfaceallowed,
	"Electron",
	"Only e+e- allowed",
	9);
	static SwitchOption interfaceallowedMuon
	(interfaceallowed,
	"Muon",
	"Only mu+mu- allowed",
	10);
	static SwitchOption interfaceallowedTau
	(interfaceallowed,
	"Tau",
	"Only tau+tau- allowed",
	11);
	static SwitchOption interfaceallowedGluon
	(interfaceallowed,
	"Gluon",
	"Only gg allowed",
	12);
	+ static SwitchOption interfaceallowedW
	+ (interfaceallowed,
	+ "W",
	+ "Only W+W- allowed",
	+ 13);
	+ static SwitchOption interfaceallowedZ
	+ (interfaceallowed,
	+ "Z",
	+ "Only ZZ allowed",
	+ 14);
	}

	void MEee2Higgs2SM::constructVertex(tSubProPtr sub) {
	// extract the particles in the hard process
	ParticleVector hard;
	hard.push_back(sub->incoming().first);hard.push_back(sub->incoming().second);
	hard.push_back(sub->outgoing()[0]);hard.push_back(sub->outgoing()[1]);
	if(hard[0]->id()<hard[1]->id()) swap(hard[0],hard[1]);
	if(hard[2]->id()<hard[3]->id()) swap(hard[2],hard[3]);
	vector<SpinorWaveFunction> fin;
	vector<SpinorBarWaveFunction> ain;
	SpinorWaveFunction( fin ,hard[0],incoming,false,true);
	SpinorBarWaveFunction(ain ,hard[1],incoming,false,true);
	double me;
	ProductionMatrixElement prodme;
	- if(hard[2]->id()!=ParticleID::g) {
	+ if(abs(hard[2]->id())<20) {
	vector<SpinorWaveFunction> aout;
	vector<SpinorBarWaveFunction> fout;
	SpinorBarWaveFunction(fout,hard[2],outgoing,true ,true);
	SpinorWaveFunction( aout,hard[3],outgoing,true ,true);
	// calculate the matrix element
	prodme=HelicityME(fin,ain,fout,aout,me);
	}
	- else {
	+ else if(abs(hard[2]->id())==21) {
	vector<VectorWaveFunction> g1,g2;
	VectorWaveFunction(g1,hard[2],outgoing,true,true);
	VectorWaveFunction(g2,hard[3],outgoing,true,true);
	g1[1]=g1[2];
	g2[1]=g2[2];
	prodme=ggME(fin,ain,g1,g2,me);
	}
	+ else {
	+ vector<VectorWaveFunction> g1,g2;
	+ VectorWaveFunction(g1,hard[2],outgoing,true,false);
	+ VectorWaveFunction(g2,hard[3],outgoing,true,false);
	+ prodme=WWME(fin,ain,g1,g2,me);
	+ }
	// construct the vertex
	HardVertexPtr hardvertex=new_ptr(HardVertex());
	// set the matrix element for the vertex
	hardvertex->ME(prodme);
	// set the pointers and to and from the vertex
	for(unsigned int ix=0;ix<4;++ix) {
	(hard[ix]->spinInfo())->productionVertex(hardvertex);
	}
	}

	void MEee2Higgs2SM::rebind(const TranslationMap & trans) {
	// dummy = trans.translate(dummy);
	FFHVertex_ = trans.translate(FFHVertex_);
	HGGVertex_ = trans.translate(HGGVertex_);
	h0_ = trans.translate(h0_);
	ME2to2Base::rebind(trans);
	}

	IVector MEee2Higgs2SM::getReferences() {
	IVector ret = ME2to2Base::getReferences();
	ret.push_back(FFHVertex_);
	ret.push_back(HGGVertex_);
	ret.push_back(h0_);
	return ret;
	}
	// the helicity amplitude matrix element
	ProductionMatrixElement MEee2Higgs2SM::ggME(vector<SpinorWaveFunction> fin,
	vector<SpinorBarWaveFunction> ain,
	vector<VectorWaveFunction> g1,
	vector<VectorWaveFunction> g2,
	double & aver) const {
	ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half,
	PDT::Spin1,PDT::Spin1);
	// wavefunctions for the intermediate particles
	ScalarWaveFunction interh;
	// temporary storage of the different diagrams
	Complex diag;
	aver=0.;
	// sum over helicities to get the matrix element
	unsigned int inhel1,inhel2,outhel1,outhel2;
	for(inhel1=0;inhel1<2;++inhel1) {
	for(inhel2=0;inhel2<2;++inhel2) {
	interh = FFHVertex_->evaluate(sHat(),1,h0_,fin[inhel1],ain[inhel2]);
	for(outhel1=0;outhel1<2;++outhel1) {
	for(outhel2=0;outhel2<2;++outhel2) {
	diag = HGGVertex_->evaluate(sHat(),g1[outhel1],g2[outhel2],interh);
	output(inhel1,inhel2,2outhel1,2outhel2)=diag;
	- aver +=real(diag*conj(diag));
	+ aver += norm(diag);
	}
	}
	}
	}
	return output;
	}
	+
	+// the helicity amplitude matrix element
	+ProductionMatrixElement MEee2Higgs2SM::WWME(vector<SpinorWaveFunction> fin,
	+ vector<SpinorBarWaveFunction> ain,
	+ vector<VectorWaveFunction> g1,
	+ vector<VectorWaveFunction> g2,
	+ double & aver) const {
	+ ProductionMatrixElement output(PDT::Spin1Half,PDT::Spin1Half,
	+ PDT::Spin1,PDT::Spin1);
	+ // wavefunctions for the intermediate particles
	+ ScalarWaveFunction interh;
	+ // temporary storage of the different diagrams
	+ Complex diag;
	+ aver=0.;
	+ // sum over helicities to get the matrix element
	+ unsigned int inhel1,inhel2,outhel1,outhel2;
	+ for(inhel1=0;inhel1<2;++inhel1) {
	+ for(inhel2=0;inhel2<2;++inhel2) {
	+ interh = FFHVertex_->evaluate(sHat(),1,h0_,fin[inhel1],ain[inhel2]);
	+ for(outhel1=0;outhel1<3;++outhel1) {
	+ for(outhel2=0;outhel2<3;++outhel2) {
	+ diag = HWWVertex_->evaluate(sHat(),g1[outhel1],g2[outhel2],interh);
	+ output(inhel1,inhel2,outhel1,outhel2)=diag;
	+ aver += norm(diag);
	+ }
	+ }
	+ }
	+ }
	+ return output;
	+}
	diff --git a/MatrixElement/Lepton/MEee2Higgs2SM.h b/MatrixElement/Lepton/MEee2Higgs2SM.h
	--- a/MatrixElement/Lepton/MEee2Higgs2SM.h
	+++ b/MatrixElement/Lepton/MEee2Higgs2SM.h
	@@ -1,236 +1,254 @@
	// -- C++ --
	//
	// MEee2Higgs2SM.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_MEee2Higgs2SM_H
	#define HERWIG_MEee2Higgs2SM_H
	//
	// This is the declaration of the MEee2Higgs2SM class.
	//

	#include "ThePEG/MatrixElement/ME2to2Base.h"
	#include "Herwig/Models/StandardModel/StandardModel.h"
	#include "Herwig/MatrixElement/ProductionMatrixElement.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* The MEee2Higgs2SM class implements the production of an \f$s\f$-channel
	* Higgs in \f$e^+e^-\f$ collisions in order to allow easy tests of Higgs
	* decays. It should not be used for physics studies.
	*
	* @see \ref MEee2Higgs2SMInterfaces "The interfaces"
	* defined for MEee2Higgs2SM.
	*/
	class MEee2Higgs2SM: public ME2to2Base {

	public:

	/**
	* The default constructor.
	*/
	inline MEee2Higgs2SM() : allowed_(0) {}

	/** @name Virtual functions required by the MEBase class. */
	//@{
	/**
	* Return the order in \f$\alpha_S\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInAlphaS() const;

	/**
	* Return the order in \f$\alpha_{EW}\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInAlphaEW() const;

	/**
	* The matrix element for the kinematical configuration
	* previously provided by the last call to setKinematics(), suitably
	* scaled by sHat() to give a dimension-less number.
	* @return the matrix element scaled with sHat() to give a
	* dimensionless number.
	*/
	virtual double me2() const;

	/**
	* Return the scale associated with the last set phase space point.
	*/
	virtual Energy2 scale() const;

	/**
	* Add all possible diagrams with the add() function.
	*/
	virtual void getDiagrams() const;

	/**
	* Get diagram selector. With the information previously supplied with the
	* setKinematics method, a derived class may optionally
	* override this method to weight the given diagrams with their
	* (although certainly not physical) relative probabilities.
	* @param dv the diagrams to be weighted.
	* @return a Selector relating the given diagrams to their weights.
	*/
	virtual Selector<DiagramIndex> diagrams(const DiagramVector & dv) const;

	/**
	* Return a Selector with possible colour geometries for the selected
	* diagram weighted by their relative probabilities.
	* @param diag the diagram chosen.
	* @return the possible colour geometries weighted by their
	* relative probabilities.
	*/
	virtual Selector<const ColourLines *>
	colourGeometries(tcDiagPtr diag) const;

	/**
	* set up the spin correlations
	*/
	virtual void constructVertex(tSubProPtr sub);
	//@}

	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);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	inline 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.
	*/
	inline virtual IBPtr fullclone() const {return new_ptr(*this);}
	//@}

	protected:

	/** @name Standard Interfaced functions. */
	//@{
	/**
	* Initialize this object after the setup phase before saving an
	* EventGenerator to disk.
	* @throws InitException if object could not be initialized properly.
	*/
	virtual void doinit();

	/**
	* Rebind pointer to other Interfaced objects. Called in the setup phase
	* after all objects used in an EventGenerator has been cloned so that
	* the pointers will refer to the cloned objects afterwards.
	* @param trans a TranslationMap relating the original objects to
	* their respective clones.
	* @throws RebindException if no cloned object was found for a given
	* pointer.
	*/
	virtual void rebind(const TranslationMap & trans)
	;

	/**
	* Return a vector of all pointers to Interfaced objects used in this
	* object.
	* @return a vector of pointers.
	*/
	virtual IVector getReferences();
	//@}

	private:

	/**
	* The matrix element
	* @param fin The incoming spinor wavefunction
	* @param ain The incoming spinorbar wavefunction
	* @param fout The outgoing spinor bar wavefunction
	* @param aout The outgoing spinor wavefunction
	* @param me The spin averaged matrix element
	*/
	ProductionMatrixElement HelicityME(vector<SpinorWaveFunction> fin,
	vector<SpinorBarWaveFunction> ain,
	vector<SpinorBarWaveFunction> fout,
	vector<SpinorWaveFunction> aout,double& me) const;
	/**
	* \f$H\to gg\f$ matrix element
	* @param fin The incoming spinor wavefunction
	* @param ain The incoming spinorbar wavefunction
	* @param g1 Outgoing gluon wavefunction
	* @param g2 Outgoing gluon wavefunction
	* @param me The spin averaged matrix element
	*/
	ProductionMatrixElement ggME(vector<SpinorWaveFunction> fin,
	vector<SpinorBarWaveFunction> ain,
	vector<VectorWaveFunction> g1,
	vector<VectorWaveFunction> g2,
	double & me) const;
	+ /**
	+ * \f$H\to WW\f$ matrix element
	+ * @param fin The incoming spinor wavefunction
	+ * @param ain The incoming spinorbar wavefunction
	+ * @param g1 Outgoing W wavefunction
	+ * @param g2 Outgoing W wavefunction
	+ * @param me The spin averaged matrix element
	+ */
	+ ProductionMatrixElement WWME(vector<SpinorWaveFunction> fin,
	+ vector<SpinorBarWaveFunction> ain,
	+ vector<VectorWaveFunction> g1,
	+ vector<VectorWaveFunction> g2,
	+ double & me) const;

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	MEee2Higgs2SM & operator=(const MEee2Higgs2SM &) = delete;

	private:

	/**
	* Pointer to the Higgs fermion-antifermion vertex
	*/
	AbstractFFSVertexPtr FFHVertex_;

	/**
	* Pointer to Higgs-gluon-gluon vertex
	*/
	AbstractVVSVertexPtr HGGVertex_;

	/**
	+ * Pointer to Higgs-WW vertex
	+ */
	+ AbstractVVSVertexPtr HWWVertex_;
	+
	+ /**
	* Allowed outgoing particles
	*/
	int allowed_;

	/**
	* Pointer to the Higgs ParticleData object
	*/
	PDPtr h0_;
	};

	}

	#endif /* HERWIG_MEee2Higgs2SM_H */
	diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
	--- a/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
	+++ b/Shower/QTilde/SplittingFunctions/HalfHalfZeroEWSplitFn.cc
	@@ -1,242 +1,248 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the HalfHalfZeroEWSplitFn class.
	//

	#include "HalfHalfZeroEWSplitFn.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 HalfHalfZeroEWSplitFn::clone() const {
	return new_ptr(*this);
	}

	IBPtr HalfHalfZeroEWSplitFn::fullclone() const {
	return new_ptr(*this);
	}

	void HalfHalfZeroEWSplitFn::persistentOutput(PersistentOStream & os) const {
	os << ghqq_ << _theSM;
	}

	void HalfHalfZeroEWSplitFn::persistentInput(PersistentIStream & is, int) {
	is >> ghqq_ >> _theSM;
	}

	// The following static variable is needed for the type description system in ThePEG.
	DescribeClass<HalfHalfZeroEWSplitFn,SplittingFunction>
	describeHerwigHalfHalfZeroEWSplitFn("Herwig::HalfHalfZeroEWSplitFn", "HwShower.so");

	void HalfHalfZeroEWSplitFn::Init() {

	static ClassDocumentation<HalfHalfZeroEWSplitFn> documentation
	("The HalfHalfZeroEWSplitFn class implements the splitting q->qH");

	}

	void HalfHalfZeroEWSplitFn::doinit() {
	SplittingFunction::doinit();
	tcSMPtr sm = generator()->standardModel();
	double sw2 = sm->sin2ThetaW();
	ghqq_ = 1./sqrt(4.*sw2);

	_theSM = dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel());
	}

	void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids) const {
	if(abs(ids[2]->id())==ParticleID::h0) {
	//get quark masses
	Energy mq;
	if(abs(ids[0]->id())==ParticleID::c)
	mq = getParticleData(ParticleID::c)->mass();
	else if(abs(ids[0]->id())==ParticleID::b)
	mq = getParticleData(ParticleID::b)->mass();
	else if(abs(ids[0]->id())==ParticleID::t)
	mq = getParticleData(ParticleID::t)->mass();
	Energy mW = getParticleData(ParticleID::Wplus)->mass();
	gH = ghqq_*(mq/mW);
	}
	else
	assert(false);
	}

	void HalfHalfZeroEWSplitFn::getCouplings(double & gH, const IdList & ids, const Energy2 t) const {
	if(abs(ids[2]->id())==ParticleID::h0) {
	//get quark masses
	Energy mq;
	if(abs(ids[0]->id())==ParticleID::c)
	mq = _theSM->mass(t,getParticleData(ParticleID::c));
	else if(abs(ids[0]->id())==ParticleID::b)
	mq = _theSM->mass(t,getParticleData(ParticleID::b));
	else if(abs(ids[0]->id())==ParticleID::t)
	mq = _theSM->mass(t,getParticleData(ParticleID::t));
	Energy mW = getParticleData(ParticleID::Wplus)->mass();
	//Energy mW = _theSM->mass(t,getParticleData(ParticleID::Wplus));
	gH = ghqq_*(mq/mW);
	}
	else
	assert(false);
	}

	double HalfHalfZeroEWSplitFn::P(const double z, const Energy2 t,
	const IdList &ids, const bool mass, const RhoDMatrix & rho) const {
	double gH(0.);
	getCouplings(gH,ids,t);
	double val = (1.-z);
	Energy mq, mH;
	//get masses
	if(mass) {
	mq = ids[0]->mass();
	mH = ids[2]->mass();
	}
	else { // to assure the particle mass in non-zero
	if(abs(ids[0]->id())==ParticleID::c)
	mq = getParticleData(ParticleID::c)->mass();
	else if(abs(ids[0]->id())==ParticleID::b)
	mq = getParticleData(ParticleID::b)->mass();
	else if(abs(ids[0]->id())==ParticleID::t)
	mq = getParticleData(ParticleID::t)->mass();
	mH = getParticleData(ParticleID::h0)->mass();
	}
	val += (4.sqr(mq) - sqr(mH))/(t(1. - z)*z);
	val *= sqr(gH);
	return colourFactor(ids)*val;
	}


	double HalfHalfZeroEWSplitFn::overestimateP(const double z,
	const IdList & ids) const {
	double gH(0.);
	getCouplings(gH,ids);
	return sqr(gH)colourFactor(ids)(1.-z);
	}

	double HalfHalfZeroEWSplitFn::ratioP(const double z, const Energy2 t,
	const IdList & ids, const bool mass,
	const RhoDMatrix & rho) const {
	double gH(0.);
	getCouplings(gH,ids,t);
	double val = 1.;
	Energy mq, mH;
	if(mass) {
	mq = ids[0]->mass();
	mH = ids[2]->mass();
	}
	else { // to assure the particle mass in non-zero
	if(abs(ids[0]->id())==ParticleID::c)
	mq = getParticleData(ParticleID::c)->mass();
	else if(abs(ids[0]->id())==ParticleID::b)
	mq = getParticleData(ParticleID::b)->mass();
	else if(abs(ids[0]->id())==ParticleID::t)
	mq = getParticleData(ParticleID::t)->mass();
	mH = getParticleData(ParticleID::h0)->mass();
	}
	val += (4.sqr(mq) - sqr(mH))/(t(1. - z)*z);
	return val;
	}

	double HalfHalfZeroEWSplitFn::integOverP(const double z,
	const IdList & ids,
	unsigned int PDFfactor) const {
	double gH(0.);
	getCouplings(gH,ids);
	double pre = colourFactor(ids)*sqr(gH);
	switch (PDFfactor) {
	case 0: //OverP
	return pre*(z-sqr(z)/2.);
	case 1: //OverP/z
	return pre*(log(z)-z);
	case 2: //OverP/(1-z)
	return pre*z;
	case 3: //OverP/[z(1-z)]
	return pre*log(z);
	default:
	throw Exception() << "HalfHalfZeroEWSplitFn::integOverP() invalid PDFfactor = "
	<< PDFfactor << Exception::runerror;
	}
	}

	double HalfHalfZeroEWSplitFn::invIntegOverP(const double r, const IdList & ids,
	unsigned int PDFfactor) const {
	double gH(0.);
	getCouplings(gH,ids);
	double pre = colourFactor(ids)*sqr(gH);
	switch (PDFfactor) {
	case 0:
	return 1. - sqrt(1. - 2.*r/pre);
	case 1: //OverP/z
	case 2: //OverP/(1-z)
	return r/pre;
	case 3: //OverP/[z(1-z)]
	return exp(r/pre);
	default:
	throw Exception() << "HalfHalfZeroEWSplitFn::invIntegOverP() invalid PDFfactor = "
	<< PDFfactor << Exception::runerror;
	}
	}

	bool HalfHalfZeroEWSplitFn::accept(const IdList &ids) const {
	if(ids.size()!=3) return false;
	if(ids[2]->id()==ParticleID::h0) {
	if(ids[0]->id()==ids[1]->id() && (ids[0]->id()==4 \|\| ids[0]->id()==5 \|\| ids[0]->id()==6))
	return true;
	}
	return false;
	}

	vector<pair<int, Complex> >
	HalfHalfZeroEWSplitFn::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> >
	HalfHalfZeroEWSplitFn::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 HalfHalfZeroEWSplitFn::matrixElement(const double z, const Energy2 t,
	const IdList & ids, const double phi,
	bool) {
	// calculate the kernal
	DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0)));
	//get masses
	Energy mq, mH;
	if(abs(ids[0]->id())==ParticleID::c)
	mq = getParticleData(ParticleID::c)->mass();
	else if(abs(ids[0]->id())==ParticleID::b)
	mq = getParticleData(ParticleID::b)->mass();
	else if(abs(ids[0]->id())==ParticleID::t)
	mq = getParticleData(ParticleID::t)->mass();
	mH = getParticleData(ParticleID::h0)->mass();
	double gH(0.);
	getCouplings(gH,ids,t);
	double mqt = mq/sqrt(t);
	double mHt = mH/sqrt(t);
	double num1 = gH(1.+z)mqt;
	double num2 = gHsqrt(-sqr(mqt)(1.-z) - sqr(mHt)z + z(1.-z)(sqr(mqt)+z(1.-z))); //watch this
	double dnum = sqrt(2.)sqrt((1.-z)sqr(z));
	Complex phase = exp(Complex(0.,1.)*phi);
	Complex cphase = conj(phase);
	(*kernal)(0,0,0) = num1/dnum;
	(kernal)(0,1,0) = cphasenum2/dnum;
	(kernal)(1,0,0) = -phasenum2/dnum;
	(*kernal)(1,1,0) = num1/dnum;
	+ cerr << "testing kernal";
	+ for(unsigned int ix=0;ix<2;++ix) {
	+ for(unsigned int iy=0;iy<2;++iy)
	+ cerr << (*kernal)(ix,iy,0) << " ";
	+ cerr << "\n";
	+ }
	// return the answer
	return kernal;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 4:09 PM (1 d, 20 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023388
Default Alt Text: (29 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions