No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
	--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
	+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.cc
	@@ -1,384 +1,386 @@
	// -- C++ --
	//
	// MatchboxAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 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 MatchboxAmplitude class.
	//

	#include "MatchboxAmplitude.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/SpinorHelicity.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/SU2Helper.h"

	using namespace Herwig;

	MatchboxAmplitude::MatchboxAmplitude()
	: Amplitude(), theColourBasisDim(0),
	calculateTrees(true), calculateLoops(true) {}

	MatchboxAmplitude::~MatchboxAmplitude() {}

	void MatchboxAmplitude::persistentOutput(PersistentOStream & os) const {
	os << theLastXComb << theColourBasis << theColourBasisDim
	<< theCrossingMap << theColourMap << theCrossingSigns
	<< theAmplitudePartonData << theNLight;
	}

	void MatchboxAmplitude::persistentInput(PersistentIStream & is, int) {
	is >> theLastXComb >> theColourBasis >> theColourBasisDim
	>> theCrossingMap >> theColourMap >> theCrossingSigns
	>> theAmplitudePartonData >> theNLight;
	}

	void MatchboxAmplitude::cloneDependencies(const std::string&) {
	}

	void MatchboxAmplitude::dumpInfo(const string& prefix) const {
	generator()->log() << prefix << fullName()
	<< " [" << this << "]\n";
	generator()->log() << prefix << " \| XComb " << lastXCombPtr()
	<< " for ";
	if ( lastXCombPtr() ) {
	for ( cPDVector::const_iterator p = lastXComb().mePartonData().begin();
	p != lastXComb().mePartonData().end(); ++p ) {
	generator()->log() << (**p).PDGName() << " ";
	}
	}
	generator()->log() << "\n";
	}

	struct orderPartonData {

	bool operator()(const pair<tcPDPtr,int>& a,
	const pair<tcPDPtr,int>& b) const {

	if ( a.first == b.first )
	return a.second < b.second;

	int acolour = a.first->iColour();
	int bcolour = b.first->iColour();

	if ( abs(acolour) != abs(bcolour) )
	return abs(acolour) < abs(bcolour);

	if ( a.first->iSpin() != b.first->iSpin() )
	return a.first->iSpin() < b.first->iSpin();

	int acharge = a.first->iCharge();
	int bcharge = b.first->iCharge();

	if ( abs(acharge) != abs(bcharge) )
	return abs(acharge) < abs(bcharge);

	if ( abs(a.first->id()) != abs(b.first->id()) )
	return abs(a.first->id()) < abs(b.first->id());

	return a.first->id() > b.first->id();

	}

	};

	void MatchboxAmplitude::fillCrossingMap(size_t shift) {

	theLastCrossingMap = crossingMap().find(lastXCombPtr());

	if ( theLastCrossingMap != crossingMap().end() ) {
	assert(amplitudePartonData().find(lastXCombPtr()) != amplitudePartonData().end());
	assert(colourMap().find(lastXCombPtr()) != colourMap().end());
	theLastAmplitudePartonData = amplitudePartonData().find(lastXCombPtr());
	theLastColourMap = colourMap().find(lastXCombPtr());
	assert(crossingSigns().find(lastXCombPtr()) != crossingSigns().end());
	theLastCrossingSign = crossingSigns().find(lastXCombPtr())->second;
	return;
	} else {
	crossingMap()[lastXCombPtr()] = vector<int>(mePartonData().size());
	amplitudePartonData()[lastXCombPtr()] = cPDVector(mePartonData().size());
	colourMap()[lastXCombPtr()] = map<size_t,size_t>();
	theLastCrossingMap = crossingMap().find(lastXCombPtr());
	theLastAmplitudePartonData = amplitudePartonData().find(lastXCombPtr());
	theLastColourMap = colourMap().find(lastXCombPtr());
	}

	double csign = 1.;
	set<pair<tcPDPtr,int>,orderPartonData > processLegs;
	for ( unsigned int l = 0; l < mePartonData().size(); ++l ) {
	if ( l > 1 )
	processLegs.insert(make_pair(mePartonData()[l],l));
	else {
	if ( mePartonData()[l]->CC() ) {
	processLegs.insert(make_pair(mePartonData()[l]->CC(),l));
	if ( mePartonData()[l]->iSpin() == PDT::Spin1Half )
	csign *= -1.;
	} else {
	processLegs.insert(make_pair(mePartonData()[l],l));
	}
	}
	}

	lastCrossingSign(csign);
	crossingSigns()[lastXCombPtr()] = csign;

	set<pair<tcPDPtr,int> > amplitudeLegs;

	int ampCount = 0;

	// process legs are already sorted, we only need to arrange for
	// adjacent particles and anti-particles
	while ( !processLegs.empty() ) {
	set<pair<tcPDPtr,int>,orderPartonData >::iterator next
	= processLegs.begin();
	while ( next->first->id() < 0 ) {
	if ( ++next == processLegs.end() )
	break;
	}
	assert(next != processLegs.end());
	lastCrossingMap()[ampCount] = next->second - shift;
	amplitudeLegs.insert(make_pair(next->first,ampCount));
	tcPDPtr check = next->first;
	processLegs.erase(next);
	++ampCount;
	if ( check->CC() ) {
	set<pair<tcPDPtr,int>,orderPartonData>::iterator checkcc
	= processLegs.end();
	for ( set<pair<tcPDPtr,int>,orderPartonData>::iterator c = processLegs.begin();
	c != processLegs.end(); ++c ) {
	if ( c->first == check->CC() ) {
	checkcc = c; break;
	}
	}
	if ( checkcc == processLegs.end() )
	for ( set<pair<tcPDPtr,int>,orderPartonData>::iterator c = processLegs.begin();
	c != processLegs.end(); ++c ) {
	assert(SU2Helper::SU2CC(check)->CC());
	if ( c->first == SU2Helper::SU2CC(check)->CC() ) {
	checkcc = c; break;
	}
	}
	assert(checkcc != processLegs.end());
	lastCrossingMap()[ampCount] = checkcc->second - shift;
	amplitudeLegs.insert(make_pair(checkcc->first,ampCount));
	processLegs.erase(checkcc);
	++ampCount;
	}
	}

	for ( set<pair<tcPDPtr,int> >::const_iterator l = amplitudeLegs.begin();
	l != amplitudeLegs.end(); ++l )
	lastAmplitudePartonData()[l->second] = l->first;

	if ( colourBasis() ) {
	assert(colourBasis()->indexMap().find(mePartonData()) !=
	colourBasis()->indexMap().end());
	const map<size_t,size_t> colourCross =
	colourBasis()->indexMap().find(mePartonData())->second;
	for ( size_t k = 0; k < lastCrossingMap().size(); ++k ) {
	if ( colourCross.find(lastCrossingMap()[k]) !=
	colourCross.end() ) {
	lastColourMap()[k] = colourCross.find(lastCrossingMap()[k])->second;
	}
	}
	}

	}

	Lorentz5Momentum MatchboxAmplitude::amplitudeMomentum(int i) const {
	int iCrossed = lastCrossingMap()[i];
	Lorentz5Momentum res = meMomenta()[iCrossed];
	return iCrossed > 1 ? res : -res;
	}

	set<vector<int> > MatchboxAmplitude::generateHelicities() const {
	set<vector<int> > res;
	vector<int> current(lastAmplitudePartonData().size());
	doGenerateHelicities(res,current,0);
	return res;
	}

	void MatchboxAmplitude::doGenerateHelicities(set<vector<int> >& res,
	vector<int>& current,
	size_t pos) const {

	if ( pos == lastAmplitudePartonData().size() ) {
	res.insert(current);
	return;
	}

	if ( lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin0 \|\|
	( lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin1 &&
	lastAmplitudePartonData()[pos]->mass() != ZERO ) ) {
	current[pos] = 0;
	doGenerateHelicities(res,current,pos+1);
	} else if ( lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin1Half \|\|
	lastAmplitudePartonData()[pos]->iSpin() == PDT::Spin1 ) {
	current[pos] = 1;
	doGenerateHelicities(res,current,pos+1);
	current[pos] = -1;
	doGenerateHelicities(res,current,pos+1);
	}

	}

	void MatchboxAmplitude::prepareAmplitudes() {

	if ( !calculateTrees )
	return;

	if ( lastAmplitudes().empty() ) {
	set<vector<int> > helicities = generateHelicities();
	for ( set<vector<int> >::const_iterator h = helicities.begin();
	h != helicities.end(); ++h )
	lastAmplitudes().insert(make_pair(*h,CVector(colourBasisDim())));
	+ lastLargeNAmplitudes() = lastAmplitudes();
	}

	- for ( AmplitudeIterator amp = lastAmplitudes().begin();
	- amp != lastAmplitudes().end(); ++amp ) {
	+ AmplitudeIterator amp = lastAmplitudes().begin();
	+ AmplitudeIterator lamp = lastLargeNAmplitudes().begin();
	+ for ( ;amp != lastAmplitudes().end(); ++amp, ++lamp ) {
	for ( size_t k = 0; k < colourBasisDim(); ++k )
	- amp->second(k) = evaluate(k,amp->first);
	+ amp->second(k) = evaluate(k,amp->first,lamp->second(k));
	}

	calculateTrees = false;

	}

	void MatchboxAmplitude::prepareOneLoopAmplitudes() {

	if ( !calculateLoops )
	return;

	if ( lastOneLoopAmplitudes().empty() ) {
	set<vector<int> > helicities = generateHelicities();
	for ( set<vector<int> >::const_iterator h = helicities.begin();
	h != helicities.end(); ++h )
	lastOneLoopAmplitudes().insert(make_pair(*h,CVector(colourBasisDim())));
	}

	for ( AmplitudeIterator amp = lastOneLoopAmplitudes().begin();
	amp != lastOneLoopAmplitudes().end(); ++amp ) {
	for ( size_t k = 0; k < colourBasisDim(); ++k )
	amp->second(k) = evaluateOneLoop(k,amp->first);
	}

	calculateLoops = false;

	}

	Complex MatchboxAmplitude::value(const tcPDVector&,
	const vector<Lorentz5Momentum>&,
	const vector<int>&) {
	assert(false && "ThePEG::Amplitude interface is not sufficient at the moment.");
	throw Exception() << "ThePEG::Amplitude interface is not sufficient at the moment."
	<< Exception::abortnow;
	return 0.;
	}

	double MatchboxAmplitude::colourCorrelatedME2(pair<int,int> ij) const {
	double Nc = generator()->standardModel()->Nc();
	double cfac = mePartonData()[ij.first]->id() == ParticleID::g ? Nc : (sqr(Nc)-1.)/(2.*Nc);
	return colourBasis()->colourCorrelatedME2(ij,mePartonData(),lastAmplitudes())/cfac;
	}

	// compare int vectors modulo certain element
	// which needs to differe between the two
	bool equalsModulo(unsigned int i, const vector<int>& a, const vector<int>& b) {
	assert(a.size()==b.size());
	if ( a[i] == b[i] )
	return false;
	for ( unsigned int k = 0; k < a.size(); ++k ) {
	if ( k == i )
	continue;
	if ( a[k] != b[k] )
	return false;
	}
	return true;
	}

	double MatchboxAmplitude::spinColourCorrelatedME2(pair<int,int> ij,
	const SpinCorrelationTensor& c) const {

	using namespace SpinorHelicity;

	Lorentz5Momentum p = meMomenta()[ij.first];
	Lorentz5Momentum n = meMomenta()[ij.second];

	LorentzVector<complex<Energy> > num =
	PlusSpinorCurrent(PlusConjugateSpinor(n),MinusSpinor(p)).eval();

	complex<Energy> den =
	sqrt(2.)*PlusSpinorProduct(PlusConjugateSpinor(n),PlusSpinor(p)).eval();

	LorentzVector<Complex> polarization(num.x()/den,num.y()/den,num.z()/den,num.t()/den);

	Complex pFactor = (polarization*c.momentum())/sqrt(abs(c.scale()));

	double avg =
	colourCorrelatedME2(ij)(-c.diagonal()+ (c.scale() > ZERO ? 1. : -1.)norm(pFactor));

	Complex corr = 0.;

	int iCrossed = -1;
	for ( unsigned int k = 0; k < lastCrossingMap().size(); ++k )
	if ( lastCrossingMap()[k] == ij.first ) {
	iCrossed = k;
	break;
	}
	assert(iCrossed >= 0);

	set<const CVector*> done;
	for ( AmplitudeConstIterator a = theLastAmplitudes.begin();
	a != theLastAmplitudes.end(); ++a ) {
	if ( done.find(&(a->second)) != done.end() )
	continue;
	AmplitudeConstIterator b = theLastAmplitudes.begin();
	while ( !equalsModulo(iCrossed,a->first,b->first) )
	if ( ++b == theLastAmplitudes.end() )
	break;
	if ( b == theLastAmplitudes.end() \|\| done.find(&(b->second)) != done.end() )
	continue;
	done.insert(&(a->second)); done.insert(&(b->second));
	if ( a->first[iCrossed] == 1 )
	swap(a,b);
	corr += colourBasis()->interference(mePartonData(),a->second,b->second);
	}

	double Nc = generator()->standardModel()->Nc();
	double cfac = mePartonData()[ij.first]->id() == ParticleID::g ? Nc : (sqr(Nc)-1.)/(2.*Nc);

	return avg + 2.(c.scale() > ZERO ? 1. : -1.)real(corr*sqr(pFactor))/cfac;

	}

	void MatchboxAmplitude::Init() {

	static ClassDocumentation<MatchboxAmplitude> documentation
	("MatchboxAmplitude is the base class for amplitude "
	"implementations inside Matchbox.");

	static Reference<MatchboxAmplitude,ColourBasis> interfaceColourBasis
	("ColourBasis",
	"Set the colour basis implementation.",
	&MatchboxAmplitude::theColourBasis, false, false, true, true, false);

	}

	// * Attention * The following static variable is needed for the type
	// description system in ThePEG. Please check that the template arguments
	// are correct (the class and its base class), and that the constructor
	// arguments are correct (the class name and the name of the dynamically
	// loadable library where the class implementation can be found).
	DescribeAbstractClass<MatchboxAmplitude,Amplitude>
	describeMatchboxAmplitude("Herwig::MatchboxAmplitude", "HwMatchbox.so");
	diff --git a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	--- a/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	+++ b/MatrixElement/Matchbox/Base/MatchboxAmplitude.h
	@@ -1,545 +1,561 @@
	// -- C++ --
	//
	// MatchboxAmplitude.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2012 The Herwig Collaboration
	//
	// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_MatchboxAmplitude_H
	#define HERWIG_MatchboxAmplitude_H
	//
	// This is the declaration of the MatchboxAmplitude class.
	//

	#include "ThePEG/MatrixElement/Amplitude.h"
	#include "ThePEG/Handlers/LastXCombInfo.h"
	#include "Herwig++/Models/StandardModel/StandardModel.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/ColourBasis.h"
	#include "Herwig++/MatrixElement/Matchbox/Utility/SpinCorrelationTensor.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief MatchboxAmplitude is the base class for amplitude
	* implementations inside Matchbox.
	*
	* @see \ref MatchboxAmplitudeInterfaces "The interfaces"
	* defined for MatchboxAmplitude.
	*/
	class MatchboxAmplitude: public Amplitude, public LastXCombInfo<StandardXComb> {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MatchboxAmplitude();

	/**
	* The destructor.
	*/
	virtual ~MatchboxAmplitude();
	//@}

	public:

	typedef map<vector<int>,CVector> AmplitudeMap;
	typedef map<vector<int>,CVector>::iterator AmplitudeIterator;
	typedef map<vector<int>,CVector>::const_iterator AmplitudeConstIterator;

	/**
	* Return the amplitude. Needs to be implemented from
	* ThePEG::Amplitude but is actually ill-defined, as colours of the
	* external particles are not specified. To this extent, this
	* implementation just asserts.
	*/
	virtual Complex value(const tcPDVector & particles,
	const vector<Lorentz5Momentum> & momenta,
	const vector<int> & helicities);

	/** @name Subprocess information */
	//@{

	/**
	* Return true, if this amplitude can handle the given process.
	*/
	virtual bool canHandle(const PDVector&) const { return false; }

	/**
	* Return the amplitude parton data.
	*/
	const cPDVector& lastAmplitudePartonData() const { return theLastAmplitudePartonData->second; }

	/**
	* Access the amplitude parton data.
	*/
	cPDVector& lastAmplitudePartonData() { return theLastAmplitudePartonData->second; }

	/**
	* Access the amplitude parton data.
	*/
	map<tStdXCombPtr,cPDVector>& amplitudePartonData() { return theAmplitudePartonData; }

	/**
	* Return the number of light flavours
	*/
	unsigned int nLight() const { return theNLight; }

	/**
	* Set the number of light flavours
	*/
	void nLight(unsigned int n) { theNLight = n; }

	/**
	* Set the (tree-level) order in \f$g_S\f$ in which this matrix
	* element should be evaluated.
	*/
	virtual void orderInGs(unsigned int) {}

	/**
	* Return the (tree-level) order in \f$g_S\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInGs() const = 0;

	/**
	* Set the (tree-level) order in \f$g_{EM}\f$ in which this matrix
	* element should be evaluated.
	*/
	virtual void orderInGem(unsigned int) {}

	/**
	* Return the (tree-level) order in \f$g_{EM}\f$ in which this matrix
	* element is given.
	*/
	virtual unsigned int orderInGem() const = 0;

	/**
	* Return the Herwig++ StandardModel object
	*/
	Ptr<StandardModel>::tcptr standardModel() {
	if ( !theStandardModel )
	theStandardModel =
	dynamic_ptr_cast<Ptr<StandardModel>::tcptr>(HandlerBase::standardModel());
	return theStandardModel;
	}

	//@}

	/** @name Colour basis. */
	//@{

	/**
	* Return the colour basis.
	*/
	Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }

	/**
	* Set the colour basis dimensionality.
	*/
	void colourBasisDim(size_t dim) { theColourBasisDim = dim; }

	/**
	* Get the colour basis dimensionality.
	*/
	size_t colourBasisDim() const { return theColourBasisDim; }

	/**
	* Return true, if this amplitude will not require colour correlations.
	*/
	virtual bool noCorrelations() const { return !haveOneLoop(); }

	/**
	* Return true, if the colour basis is capable of assigning colour
	* flows.
	*/
	virtual bool haveColourFlows() const {
	return colourBasis() ? colourBasis()->haveColourFlows() : false;
	}

	/**
	* Return a Selector with possible colour geometries for the selected
	* diagram weighted by their relative probabilities.
	*/
	virtual Selector<const ColourLines *> colourGeometries(tcDiagPtr diag) const {
	return
	haveColourFlows() ?
	- theColourBasis->colourGeometries(diag,lastAmplitudes()) :
	+ theColourBasis->colourGeometries(diag,lastLargeNAmplitudes()) :
	Selector<const ColourLines *>();
	}

	/**
	* Return the colour crossing information as filled by the last call to
	* fillCrossingMap(...), mapping amplitude ids to colour basis ids.
	*/
	const map<size_t,size_t>& lastColourMap() const { return theLastColourMap->second; }

	/**
	* Access the colour crossing information.
	*/
	map<size_t,size_t>& lastColourMap() { return theLastColourMap->second; }

	/**
	* Access the colour crossing information.
	*/
	map<tStdXCombPtr,map<size_t,size_t> >& colourMap() { return theColourMap; }

	//@}

	/** @name Phasespace point, crossing and helicities */
	//@{

	/**
	* Set the xcomb object.
	*/
	virtual void setXComb(tStdXCombPtr xc) {
	theLastXComb = xc;
	fillCrossingMap();
	}

	/**
	* Return the momentum as crossed appropriate for this amplitude.
	*/
	Lorentz5Momentum amplitudeMomentum(int) const;

	/**
	* Perform a normal ordering of external legs and fill the
	* crossing information as. This default implementation sorts
	* lexicographically in (abs(colour)/spin/abs(charge)), putting pairs
	* of particles/anti-particles where possible.
	*/
	virtual void fillCrossingMap(size_t shift = 0);

	/**
	* Return the crossing sign.
	*/
	double lastCrossingSign() const { return theLastCrossingSign; }

	/**
	* Set the crossing sign.
	*/
	void lastCrossingSign(double s) { theLastCrossingSign = s; }

	/**
	* Return the crossing information as filled by the last call to
	* fillCrossingMap(...), mapping amplitude ids to process ids.
	*/
	const vector<int>& lastCrossingMap() const { return theLastCrossingMap->second; }

	/**
	* Access the crossing information.
	*/
	vector<int>& lastCrossingMap() { return theLastCrossingMap->second; }

	/**
	* Access the crossing information.
	*/
	map<tStdXCombPtr,vector<int> >& crossingMap() { return theCrossingMap; }

	/**
	* Access the crossing signs.
	*/
	map<tStdXCombPtr,double>& crossingSigns() { return theCrossingSigns; }

	/**
	* Generate the helicity combinations.
	*/
	virtual set<vector<int> > generateHelicities() const;

	//@}

	/** @name Tree-level amplitudes */
	//@{

	/**
	* Calculate the tree level amplitudes for the phasespace point
	* stored in lastXComb.
	*/
	virtual void prepareAmplitudes();

	/**
	* Return last evaluated helicity amplitudes.
	*/
	const AmplitudeMap& lastAmplitudes() const { return theLastAmplitudes; }

	/**
	* Access the last evaluated helicity amplitudes.
	*/
	AmplitudeMap& lastAmplitudes() { return theLastAmplitudes; }

	/**
	+ * Return last evaluated, leading colour helicity amplitudes.
	+ */
	+ const AmplitudeMap& lastLargeNAmplitudes() const { return theLastLargeNAmplitudes; }
	+
	+ /**
	+ * Access the last evaluated, leading colour helicity amplitudes.
	+ */
	+ AmplitudeMap& lastLargeNAmplitudes() { return theLastLargeNAmplitudes; }
	+
	+ /**
	* Return the matrix element squared.
	*/
	virtual double me2() const {
	return colourBasis()->me2(mePartonData(),lastAmplitudes());
	}

	/**
	* Return the colour correlated matrix element.
	*/
	virtual double colourCorrelatedME2(pair<int,int> ij) const;

	/**
	* Return the colour and spin correlated matrix element.
	*/
	virtual double spinColourCorrelatedME2(pair<int,int> emitterSpectator,
	const SpinCorrelationTensor& c) const;


	/**
	* Evaluate the amplitude for the given colour tensor id and
	* helicity assignment
	*/
	- virtual Complex evaluate(size_t, const vector<int>&) { return 0.; }
	+ virtual Complex evaluate(size_t, const vector<int>&, Complex&) { return 0.; }

	//@}

	/** @name One-loop amplitudes */
	//@{

	/**
	* Return true, if this amplitude is capable of calculating one-loop
	* (QCD) corrections.
	*/
	virtual bool haveOneLoop() const { return false; }

	/**
	* Return true, if this amplitude only provides
	* one-loop (QCD) corrections.
	*/
	virtual bool onlyOneLoop() const { return false; }

	/**
	* Return true, if one loop corrections have been calculated in
	* dimensional reduction. Otherwise conventional dimensional
	* regularization is assumed. Note that renormalization is always
	* assumed to be MSbar.
	*/
	bool isDR() const { return false; }

	/**
	* Return true, if one loop corrections are given in the conventions
	* of the integrated dipoles.
	*/
	bool isCS() const { return false; }

	/**
	* Return the value of the dimensional regularization
	* parameter. Note that renormalization scale dependence is fully
	* restored in DipoleIOperator.
	*/
	virtual Energy2 mu2() const { return 0.*GeV2; }

	/**
	* Calculate the one-loop amplitudes for the phasespace point
	* stored in lastXComb, if provided.
	*/
	virtual void prepareOneLoopAmplitudes();

	/**
	* Return last evaluated one-loop helicity amplitudes.
	*/
	const AmplitudeMap& lastOneLoopAmplitudes() const { return theLastOneLoopAmplitudes; }

	/**
	* Access the last evaluated one-loop helicity amplitudes.
	*/
	AmplitudeMap& lastOneLoopAmplitudes() { return theLastOneLoopAmplitudes; }

	/**
	* Return the one-loop/tree interference.
	*/
	virtual double oneLoopInterference() const {
	return colourBasis()->interference(mePartonData(),
	lastOneLoopAmplitudes(),lastAmplitudes());
	}

	/**
	* Evaluate the amplitude for the given colour tensor id and
	* helicity assignment
	*/
	virtual Complex evaluateOneLoop(size_t, const vector<int>&) { return 0.; }

	//@}

	/** @name Caching and helpers to setup amplitude objects. */
	//@{

	/**
	* Flush all cashes.
	*/
	virtual void flushCaches() {
	calculateTrees = true;
	calculateLoops = true;
	}

	/**
	* Clone this amplitude.
	*/
	Ptr<MatchboxAmplitude>::ptr cloneMe() const {
	return dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(clone());
	}

	/**
	* Clone the dependencies, using a given prefix.
	*/
	virtual void cloneDependencies(const std::string& prefix = "");

	//@}

	/** @name Diagnostic information */
	//@{

	/**
	* Dump xcomb hierarchies.
	*/
	void dumpInfo(const string& prefix = "") const;

	//@}

	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();

	// If needed, insert declarations of virtual function defined in the
	// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).


	private:

	/**
	* Recursively generate helicities
	*/
	void doGenerateHelicities(set<vector<int> >& res,
	vector<int>& current,
	size_t pos) const;

	/**
	* The Herwig++ StandardModel object
	*/
	Ptr<StandardModel>::tcptr theStandardModel;

	/**
	* The number of light flavours to be used.
	*/
	unsigned int theNLight;

	/**
	* The colour basis implementation to be used.
	*/
	Ptr<ColourBasis>::ptr theColourBasis;

	/**
	* The dimensionality of the colour basis for the processes covered
	* by the colour basis.
	*/
	size_t theColourBasisDim;

	/**
	* References to the amplitude values which have been contributing
	* to the last call of prepareAmplitudes.
	*/
	map<vector<int>,CVector> theLastAmplitudes;

	/**
	+ * References to the leading N amplitude values which have been
	+ * contributing to the last call of prepareAmplitudes.
	+ */
	+ map<vector<int>,CVector> theLastLargeNAmplitudes;
	+
	+ /**
	* References to the one-loop amplitude values which have been contributing
	* to the last call of prepareAmplitudes.
	*/
	map<vector<int>,CVector> theLastOneLoopAmplitudes;

	/**
	* The crossing information as filled by the last call to
	* fillCrossingMap()
	*/
	map<tStdXCombPtr,vector<int> > theCrossingMap;

	/**
	* The colour crossing information as filled by the last call to
	* fillCrossingMap()
	*/
	map<tStdXCombPtr,map<size_t,size_t> > theColourMap;

	/**
	* The crossing signs as filled by the last call to
	* fillCrossingMap()
	*/
	map<tStdXCombPtr,double> theCrossingSigns;

	/**
	* The amplitude parton data.
	*/
	map<tStdXCombPtr,cPDVector> theAmplitudePartonData;

	/**
	* The crossing information as filled by the last call to
	* fillCrossingMap()
	*/
	map<tStdXCombPtr,vector<int> >::iterator theLastCrossingMap;

	/**
	* The colour crossing information as filled by the last call to
	* fillCrossingMap()
	*/
	map<tStdXCombPtr,map<size_t,size_t> >::iterator theLastColourMap;

	/**
	* The amplitude parton data.
	*/
	map<tStdXCombPtr,cPDVector>::iterator theLastAmplitudePartonData;

	/**
	* The crossing sign.
	*/
	double theLastCrossingSign;

	/**
	* True, if tree amplitudes need to be recalculated.
	*/
	bool calculateTrees;

	/**
	* True, if loop amplitudes need to be recalculated.
	*/
	bool calculateLoops;

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

	};

	}

	#endif /* HERWIG_MatchboxAmplitude_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 8:42 PM (1 d, 3 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806129
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions