No OneTemporary
Actions

Size

32 KB

Subscribers

None

View Options

	diff --git a/MatrixElement/Lepton/MEee2Mesons.cc b/MatrixElement/Lepton/MEee2Mesons.cc
	--- a/MatrixElement/Lepton/MEee2Mesons.cc
	+++ b/MatrixElement/Lepton/MEee2Mesons.cc
	@@ -1,400 +1,183 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the MEee2Mesons class.
	//

	#include "MEee2Mesons.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 "Herwig/Decay/DecayIntegrator2.fh"
	+#include "Herwig/Decay/PhaseSpaceMode.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/EnumParticles.h"
	#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
	-#include "ThePEG/Utilities/SimplePhaseSpace.h"
	-#include "ThePEG/Cuts/Cuts.h"
	+#include "ThePEG/StandardModel/StandardModelBase.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
	#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
	-#include "ThePEG/StandardModel/StandardModelBase.h"
	-#include "Herwig/Decay/DecayIntegrator.h"
	#include "ThePEG/PDF/PolarizedBeamParticleData.h"

	using namespace Herwig;

	typedef LorentzVector<complex<InvEnergy> > LorentzPolarizationVectorInvE;

	-void MEee2Mesons::getDiagrams() const {
	- // make sure the current got initialised
	- current_->init();
	- tPDPtr gamma = getParticleData(ParticleID::gamma);
	- Energy Emax = generator()->maximumCMEnergy();
	- tcPDPtr em = getParticleData(ParticleID::eminus);
	- tcPDPtr ep = getParticleData(ParticleID::eplus);
	- // loop at the modes
	- for(unsigned int imode=0;imode<current_->numberOfModes();++imode) {
	- // get the external particles for this mode
	- int iq(0),ia(0);
	- tPDVector extpart = {gamma};
	- tPDVector ptemp = current_->particles(0,imode,iq,ia);
	- extpart.insert(std::end(extpart), std::begin(ptemp), std::end(ptemp));
	- // create the mode
	- DecayPhaseSpaceModePtr mode = new_ptr(DecayPhaseSpaceMode(extpart,DecayIntegratorPtr()));
	- // create the first piece of the channel
	- DecayPhaseSpaceChannelPtr channel = new_ptr(DecayPhaseSpaceChannel(mode));
	- channel->addIntermediate(extpart[0],0,0.0,-1,1);
	- if(!current_->createMode(0,imode,mode,1,1,channel,Emax)) continue;
	- nmult_ = int(extpart.size())-1;
	- int ndiag=0;
	- for(unsigned int imode=0;imode<mode->numberChannels();++imode) {
	- tcDecayPhaseSpaceChannelPtr iChannel = mode->channel(imode);
	- // extract the intermediates
	- vector<pair<int,int> > children;
	- vector<tcPDPtr> intermediate;
	- vector<unsigned int> jacType;
	- vector<Energy> intMass;
	- vector<Energy> intWidth;
	- vector<double> intPower;
	- for(unsigned int inter=1;inter<iChannel->numberOfIntermediates();++inter) {
	- children.push_back(make_pair(0,0));
	- intermediate.push_back(tcPDPtr());
	- jacType.push_back(0);
	- intMass.push_back(ZERO);
	- intWidth.push_back(ZERO);
	- intPower.push_back(0.);
	- iChannel->intermediateInfo(inter,intermediate.back(),jacType.back(),
	- intMass.back(),intWidth.back(),intPower.back(),
	- children.back().first,children.back().second);
	- }
	- unsigned int isize=2;
	- map<unsigned,unsigned int> ires;
	- // create the diagram
	- ThePEG::Ptr<ThePEG::Tree2toNDiagram>::pointer diag;
	- for(unsigned int ix=0;ix<intermediate.size();++ix) {
	- if(ix==0) {
	- diag = new_ptr((Tree2toNDiagram(2), em, ep,1,intermediate[ix]));
	- ndiag+=1;
	- isize+=1;
	- ires[ix] = isize;
	- }
	- if(children[ix].first>0) {
	- diag = new_ptr((*diag,ires[ix],extpart[children[ix].first]));
	- isize+=1;
	- }
	- else {
	- int iloc = -children[ix].first-1;
	- isize+=1;
	- ires[iloc] = isize;
	- diag = new_ptr((*diag,ires[ix],intermediate[iloc]));
	- }
	- if(children[ix].second>0) {
	- diag = new_ptr((*diag,ires[ix],extpart[children[ix].second]));
	- isize+=1;
	- }
	- else {
	- int iloc = -children[ix].second-1;
	- isize+=1;
	- ires[iloc] = isize;
	- diag = new_ptr((*diag,ires[ix],intermediate[iloc]));
	- }
	- }
	- diag = new_ptr((*diag,-ndiag));
	- add(diag);
	- }
	- }
	-}
	+MEee2Mesons::MEee2Mesons() {}

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

	-int MEee2Mesons::nDim() const {
	- return 3*nmult_-5;
	+unsigned int MEee2Mesons::orderInAlphaS() const {
	+ return 0;
	}

	-void MEee2Mesons::setKinematics() {
	- HwMEBase::setKinematics();
	+unsigned int MEee2Mesons::orderInAlphaEW() const {
	+ return 0;
	}

	-bool MEee2Mesons::generateKinematics(const double * r) {
	- using Constants::pi;
	- // Save the jacobian dPS/dr for later use.
	- jacobian(1.0);
	- // set the masses of the outgoing particles
	- for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
	- meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass());
	- }
	- double ctmin = -1.0, ctmax = 1.0;
	- double cth = getCosTheta(ctmin, ctmax, r[0]);
	- phi(rnd(2.0*Constants::pi));
	- unsigned int i1(2),i2(3),i3(4);
	- Energy q = ZERO;
	- Energy e = sqrt(sHat())/2.0;
	- Energy2 pq;
	- if(nDim()==1) {
	- try {
	- q = SimplePhaseSpace::
	- getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass());
	- }
	- catch ( ImpossibleKinematics ) {
	- return false;
	- }
	- pq = 2.0eq;
	-
	- Energy pt = q*sqrt(1.0-sqr(cth));
	- meMomenta()[2].setVect(Momentum3( ptsin(phi()), ptcos(phi()), q*cth));
	- meMomenta()[3].setVect(Momentum3(-ptsin(phi()), -ptcos(phi()), -q*cth));
	-
	- meMomenta()[2].rescaleEnergy();
	- meMomenta()[3].rescaleEnergy();
	-
	- Energy2 m22 = meMomenta()[2].mass2();
	- Energy2 m32 = meMomenta()[3].mass2();
	- Energy2 e0e2 = 2.0esqrt(sqr(q) + m22);
	- tHat(pq*cth + m22 - e0e2);
	- uHat(m22 + m32 - sHat() - tHat());
	- }
	- else if(nDim()==4) {
	- double rm=r[1];
	- if(rm<1./3.) {
	- rm *=3.;
	- }
	- else if(rm<2./3.) {
	- rm = 3.*rm-1.;
	- swap(i2,i3);
	- }
	- else {
	- rm = 3.*rm-2.;
	- swap(i1,i2);
	- }
	- tcPDPtr res = getParticleData(213);
	- Energy mass = res->mass(), width = res->width();
	- Energy2 m2max = sqr(sqrt(sHat())-meMomenta()[i3].mass());
	- Energy2 m2min = sqr(meMomenta()[i1].mass()+meMomenta()[i2].mass());
	- double rhomin = atan((m2min-sqr(mass))/mass/width);
	- double rhomax = atan((m2max-sqr(mass))/mass/width);
	- double rho = rhomin+rm*(rhomax-rhomin);
	- Energy2 m2 = mass(widthtan(rho)+mass);
	- Energy mv = sqrt(m2);
	- try {
	- q = SimplePhaseSpace::
	- getMagnitude(sHat(), mv, meMomenta()[i3].mass());
	- }
	- catch ( ImpossibleKinematics ) {
	- return false;
	- }
	- pq = 2.0eq;
	-
	- Energy pt = q*sqrt(1.0-sqr(cth));
	- Lorentz5Momentum poff;
	- poff.setMass(mv);
	- poff.setVect(Momentum3( ptsin(phi()), ptcos(phi()), q*cth));
	- meMomenta()[i3].setVect(Momentum3(-ptsin(phi()), -ptcos(phi()), -q*cth));
	-
	- poff.rescaleEnergy();
	- meMomenta()[i3].rescaleEnergy();
	- // decay of the intermediate
	- bool test=Kinematics::twoBodyDecay(poff,meMomenta()[i1].mass(),
	- meMomenta()[i2].mass(),
	- -1.+2r[2],r[3]2.*pi,
	- meMomenta()[i1],meMomenta()[i2]);
	- if(!test) return false;
	- // decay piece of the jacobian
	- Energy p2 = Kinematics::pstarTwoBodyDecay(mv,meMomenta()[i1].mass(),
	- meMomenta()[i2].mass());
	- jacobian(p2/mv/8./sqr(pi)*jacobian());
	- // mass piece
	- jacobian((rhomax-rhomin)( sqr(m2-sqr(mass))+sqr(masswidth))
	- /mass/width*jacobian()/sHat());
	- }
	- else {
	- assert(false);
	- }
	- vector<LorentzMomentum> out(meMomenta().size()-2);
	- tcPDVector tout(meMomenta().size());
	- for(unsigned int ix=2;ix<meMomenta().size();++ix) {
	- out[ix-2] = meMomenta()[ix];
	- tout[ix-2] = mePartonData()[ix];
	- }
	- if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) )
	- return false;
	- jacobian((pq/sHat())Constants::pijacobian());
	- return true;
	+IBPtr MEee2Mesons::clone() const {
	+ return new_ptr(*this);
	}

	-double MEee2Mesons::me2() const {
	- // cerr << "testing in me2\n";
	- // for(unsigned int ix=0;ix<mePartonData().size();++ix)
	- // cerr << mePartonData()[ix]->PDGName() << " " << meMomenta()[ix]/GeV << "\n";
	- // wavefunctions for the incoming leptons
	+IBPtr MEee2Mesons::fullclone() const {
	+ return new_ptr(*this);
	+}
	+
	+void MEee2Mesons::doinit() {
	+ // make sure the current got initialised
	+ current_->init();
	+ // max energy
	+ Energy Emax = generator()->maximumCMEnergy();
	+ // incoming particles
	+ tPDPtr em = getParticleData(ParticleID::eminus);
	+ tPDPtr ep = getParticleData(ParticleID::eplus);
	+ // loop over the modes
	+ int nmode=0;
	+ for(unsigned int imode=0;imode<current_->numberOfModes();++imode) {
	+ // get the external particles for this mode
	+ int iq(0),ia(0);
	+ tPDVector out = current_->particles(0,imode,iq,ia);
	+ current_->decayModeInfo(imode,iq,ia);
	+ if(iq==2&&ia==-2) continue;
	+ PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(em,out,1.,ep,Emax));
	+ PhaseSpaceChannel channel(mode);
	+ if(!current_->createMode(0,tcPDPtr(), IsoSpin::IUnknown, IsoSpin::I3Unknown,
	+ imode,mode,0,-1,channel,Emax)) continue;
	+ modeMap_[imode] = nmode;
	+ addMode(mode);
	+ ++nmode;
	+ }
	+ MEMultiChannel::doinit();
	+}
	+
	+void MEee2Mesons::persistentOutput(PersistentOStream & os) const {
	+os << current_ << modeMap_;
	+}
	+
	+void MEee2Mesons::persistentInput(PersistentIStream & is, int) {
	+is >> current_ >> modeMap_;
	+}
	+
	+
	+// The following static variable is needed for the type
	+// description system in ThePEG.
	+DescribeClass<MEee2Mesons,MEMultiChannel>
	+describeHerwigMEee2Mesons("Herwig::MEee2Mesons", "HwMELeptonLowEnergy.so");
	+
	+void MEee2Mesons::Init() {
	+
	+ static ClassDocumentation<MEee2Mesons> documentation
	+ ("The MEee2Mesons class simluation the production of low multiplicity"
	+ " events via the weak current");
	+
	+ static Reference<MEee2Mesons,WeakCurrent> interfaceWeakCurrent
	+ ("WeakCurrent",
	+ "The reference for the decay current to be used.",
	+ &MEee2Mesons::current_, false, false, true, false, false);
	+
	+}
	+
	+double MEee2Mesons::me2(const int ichan) const {
	SpinorWaveFunction em_in( meMomenta()[0],mePartonData()[0],incoming);
	SpinorBarWaveFunction ep_in( meMomenta()[1],mePartonData()[1],incoming);
	vector<SpinorWaveFunction> f1;
	vector<SpinorBarWaveFunction> a1;
	for(unsigned int ix=0;ix<2;++ix) {
	em_in.reset(ix);
	f1.push_back(em_in);
	ep_in.reset(ix);
	a1.push_back(ep_in);
	}
	// compute the leptonic current
	LorentzPolarizationVectorInvE lepton[2][2];
	InvEnergy2 pre = SM().alphaEM(sHat())4.Constants::pi/sHat();
	for(unsigned ix=0;ix<2;++ix) {
	for(unsigned iy=0;iy<2;++iy) {
	lepton[ix][iy]= pre*f1[ix].dimensionedWave().vectorCurrent(a1[iy].dimensionedWave());
	}
	}
	// work out the mapping for the hadron vector
	vector<unsigned int> constants(meMomenta().size()+1);
	vector<PDT::Spin > ispin(meMomenta().size()-2);
	vector<int> hadrons(meMomenta().size()-2);
	int itemp(1);
	unsigned int ix(meMomenta().size());
	do {
	--ix;
	ispin[ix-2] = mePartonData()[ix]->iSpin();
	hadrons[ix-2] = mePartonData()[ix]->id();
	itemp *= ispin[ix-2];
	constants[ix] = itemp;
	}
	while(ix>2);
	constants[meMomenta().size()] = 1;
	constants[0] = constants[1] = constants[2];
	- // cerr << "testing constants\n";
	- // for(unsigned int ix=0;ix<constants.size();++ix)
	- // cerr << constants[ix] << " ";
	- // cerr << "\n";
	- // cerr << "testing ispn\n";
	- // for(unsigned int ix=0;ix<ispin.size();++ix)
	- // cerr << ispin[ix] << " ";
	- // cerr << "\n";
	// calculate the matrix element
	me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,ispin));
	// calculate the hadron current
	unsigned int imode = current_->decayMode(hadrons);
	Energy q = sqrt(sHat());
	- // cerr << "testing mode " << imode << " " << q/GeV << "\n";
	- ParticleVector hadpart;
	- for(unsigned int ix=2;ix<mePartonData().size();++ix) {
	- hadpart.push_back(mePartonData()[ix]->produceParticle(meMomenta()[ix]));
	- }
	+ vector<Lorentz5Momentum> momenta(meMomenta() .begin()+2, meMomenta().end());
	+ tPDVector out = mode(modeMap_.at(imode))->outgoing();
	+ if(ichan<0) iMode(modeMap_.at(imode));
	vector<LorentzPolarizationVectorE>
	- hadron(current_->current(imode,-1,q,hadpart,DecayIntegrator::Calculate));
	- // for(unsigned int ix=0;ix<hadron.size();++ix)
	- // cerr << hadron[ix].x()/GeV << " " << hadron[ix].y()/GeV << " " << hadron[ix].z()/GeV << " " << hadron[ix].t()/GeV << "\n";
	+ hadron(current_->current(tcPDPtr(), IsoSpin::IUnknown, IsoSpin::I3Unknown, imode,ichan,
	+ q,out,momenta,DecayIntegrator2::Calculate));
	// compute the matrix element
	vector<unsigned int> ihel(meMomenta().size());
	double output(0.);
	for(unsigned int hhel=0;hhel<hadron.size();++hhel) {
	// map the index for the hadrons to a helicity state
	for(unsigned int ix=meMomenta().size()-1;ix>1;--ix) {
	ihel[ix]=(hhel%constants[ix-1])/constants[ix];
	}
	// loop over the helicities of the incoming leptons
	for(ihel[1]=0;ihel[1]<2;++ihel[1]){
	for(ihel[0]=0;ihel[0]<2;++ihel[0]) {
	Complex amp = lepton[ihel[0]][ihel[1]].dot(hadron[hhel]);
	me_(ihel)= amp;
	output += std::norm(amp);
	- // cerr << "testing ME ";
	- // for(unsigned int ix=0;ix<ihel.size();++ix) cerr << ihel[ix] << " ";
	- // cerr << me_(ihel) << "\n";
	}
	}
	}
	- output = 0.25sqr(pow(sqrt(sHat())/q,int(hadpart.size()-2)));
	+ // prefactors
	+ output = 0.25sqr(pow(sqrt(sHat())/q,int(momenta.size()-2)));
	// polarization stuff
	tcPolarizedBeamPDPtr beam[2] =
	{dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[0]),
	dynamic_ptr_cast<tcPolarizedBeamPDPtr>(mePartonData()[1])};
	if( beam[0] \|\| beam[1] ) {
	RhoDMatrix rho[2] = {beam[0] ? beam[0]->rhoMatrix() : RhoDMatrix(mePartonData()[0]->iSpin()),
	beam[1] ? beam[1]->rhoMatrix() : RhoDMatrix(mePartonData()[1]->iSpin())};
	output = me_.average(rho[0],rho[1]);
	}
	return output;
	}

	-CrossSection MEee2Mesons::dSigHatDR() const {
	- return me2()jacobian()/(16.0sqr(Constants::pi)sHat())sqr(hbarc);
	+void MEee2Mesons::constructVertex(tSubProPtr) {
	}
	-
	-unsigned int MEee2Mesons::orderInAlphaS() const {
	- return 0;
	-}
	-
	-unsigned int MEee2Mesons::orderInAlphaEW() const {
	- return 0;
	-}
	-
	-Selector<MEBase::DiagramIndex>
	-MEee2Mesons::diagrams(const DiagramVector & diags) const {
	- // This example corresponds to the diagrams specified in the example
	- // in the getDiagrams() function.
	-
	- Selector<DiagramIndex> sel;
	- for ( DiagramIndex i = 0; i < diags.size(); ++i )
	- if ( diags[i]->id() == -1 ) sel.insert(1.0, i);
	- else if ( diags[i]->id() == -2 ) sel.insert(1.0, i);
	- else if ( diags[i]->id() == -3 ) sel.insert(1.0, i);
	- // You probably do not want equal weights here...
	- return sel;
	-
	- // If there is only one possible diagram you can override the
	- // MEBase::diagram function instead.
	-
	-}
	-
	-Selector<const ColourLines *>
	-MEee2Mesons::colourGeometries(tcDiagPtr) const {
	- static ColourLines none("");
	- Selector<const ColourLines *> sel;
	- sel.insert(1.0, &none);
	- return sel;
	-}
	-
	-IBPtr MEee2Mesons::clone() const {
	- return new_ptr(*this);
	-}
	-
	-IBPtr MEee2Mesons::fullclone() const {
	- return new_ptr(*this);
	-}
	-
	-void MEee2Mesons::doinit() {
	- HwMEBase::doinit();
	-}
	-
	-void MEee2Mesons::persistentOutput(PersistentOStream & os) const {
	- os << current_ << nmult_;
	-}
	-
	-void MEee2Mesons::persistentInput(PersistentIStream & is, int) {
	- is >> current_ >> nmult_;
	-}
	-
	-// The following static variable is needed for the type
	-// description system in ThePEG.
	-DescribeClass<MEee2Mesons,HwMEBase>
	- describeHerwigMEee2Mesons("Herwig::MEee2Mesons", "HwMELeptonLowEnergy.so");
	-
	-void MEee2Mesons::Init() {
	-
	- static ClassDocumentation<MEee2Mesons> documentation
	- ("The MEee2Mesons class simluation the production of low multiplicity"
	- " events via the weak current");
	-
	- static Reference<MEee2Mesons,WeakDecayCurrent> interfaceWeakCurrent
	- ("WeakCurrent",
	- "The reference for the decay current to be used.",
	- &MEee2Mesons::current_, false, false, true, false, false);
	-
	-}
	-
	-
	-void MEee2Mesons::constructVertex(tSubProPtr sub) {
	-}
	diff --git a/MatrixElement/Lepton/MEee2Mesons.h b/MatrixElement/Lepton/MEee2Mesons.h
	--- a/MatrixElement/Lepton/MEee2Mesons.h
	+++ b/MatrixElement/Lepton/MEee2Mesons.h
	@@ -1,207 +1,148 @@
	// -- C++ --
	#ifndef Herwig_MEee2Mesons_H
	#define Herwig_MEee2Mesons_H
	//
	// This is the declaration of the MEee2Mesons class.
	//

	-#include "Herwig/MatrixElement/HwMEBase.h"
	-#include "Herwig/Decay/WeakCurrents/WeakDecayCurrent.h"
	+#include "Herwig/MatrixElement/MEMultiChannel.h"
	+#include "Herwig/Decay/WeakCurrents/WeakCurrent.h"
	#include "Herwig/MatrixElement/ProductionMatrixElement.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	- * The MEee2Mesons class is designed to simulate the production of a small
	- * number of mesons via the weak current
	+ * Here is the documentation of the MEee2Mesons class.
	*
	* @see \ref MEee2MesonsInterfaces "The interfaces"
	* defined for MEee2Mesons.
	*/
	-class MEee2Mesons: public HwMEBase {
	+class MEee2Mesons: public MEMultiChannel {

	public:

	/**
	* The default constructor.
	*/
	- MEee2Mesons()
	- {}
	+ MEee2Mesons();

	public:

	/** @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;
	-
	- /**
	- * Set the typed and momenta of the incoming and outgoing partons to
	- * be used in subsequent calls to me() and colourGeometries()
	- * according to the associated XComb object. If the function is
	- * overridden in a sub class the new function must call the base
	- * class one first.
	- */
	- virtual void setKinematics();
	-
	- /**
	- * The number of internal degrees of freedom used in the matrix
	- * element.
	- */
	- virtual int nDim() const;
	-
	- /**
	- * Generate internal degrees of freedom given nDim() uniform
	- * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
	- * generator, the dSigHatDR should be a smooth function of these
	- * numbers, although this is not strictly necessary.
	- * @param r a pointer to the first of nDim() consecutive random numbers.
	- * @return true if the generation succeeded, otherwise false.
	- */
	- virtual bool generateKinematics(const double * r);
	-
	- /**
	- * Return the matrix element squared differential in the variables
	- * given by the last call to generateKinematics().
	- */
	- virtual CrossSection dSigHatDR() 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;
	+ //@}

	/**
	* Construct the vertex of spin correlations.
	*/
	virtual void constructVertex(tSubProPtr);
	- //@}


	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:
	+
	+ /**
	+ * Return the matrix element squared for a given mode and phase-space channel.
	+ * @param ichan The channel we are calculating the matrix element for.
	+ */
	+ virtual double me2(const int ichan) const;
	+
	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	virtual IBPtr fullclone() const;
	//@}

	protected:
	-
	+
	/**
	* 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();
	-
	-
	+
	private:

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

	-private:
	-
	+private :
	+
	/**
	* the hadronic current
	*/
	- WeakDecayCurrentPtr current_;
	-
	- /**
	- * Multiplicity
	- */
	- mutable unsigned int nmult_;
	+ WeakCurrentPtr current_;

	/**
	* The matrix element
	*/
	mutable ProductionMatrixElement me_;
	+
	+/**
	+ * Map for the modes
	+ */
	+map<int,int> modeMap_;
	};

	}

	#endif /* Herwig_MEee2Mesons_H */
	diff --git a/MatrixElement/MEMultiChannel.cc b/MatrixElement/MEMultiChannel.cc
	new file mode 100644
	--- /dev/null
	+++ b/MatrixElement/MEMultiChannel.cc
	@@ -0,0 +1,144 @@
	+// -- C++ --
	+//
	+// This is the implementation of the non-inlined, non-templated member
	+// functions of the MEMultiChannel class.
	+//
	+
	+#include "MEMultiChannel.h"
	+#include "Herwig/Decay/DecayIntegrator2.fh"
	+#include "Herwig/Decay/PhaseSpaceMode.h"
	+#include "ThePEG/Interface/ClassDocumentation.h"
	+#include "ThePEG/EventRecord/Particle.h"
	+#include "ThePEG/Repository/UseRandom.h"
	+#include "ThePEG/Repository/EventGenerator.h"
	+#include "ThePEG/Utilities/DescribeClass.h"
	+#include "ThePEG/Persistency/PersistentOStream.h"
	+#include "ThePEG/Persistency/PersistentIStream.h"
	+#include "ThePEG/PDT/EnumParticles.h"
	+#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
	+#include "ThePEG/Cuts/Cuts.h"
	+
	+using namespace Herwig;
	+
	+MEMultiChannel::~MEMultiChannel() {}
	+
	+void MEMultiChannel::getDiagrams() const {
	+ int ndiag=0;
	+ channelMap_.clear();
	+ for(PhaseSpaceModePtr mode : modes_) {
	+ channelMap_.push_back(vector<int>());
	+ for(const PhaseSpaceChannel & channel : mode->channels()) {
	+ ThePEG::Ptr<ThePEG::Tree2toNDiagram>::pointer diag = channel.createDiagram();
	+ ndiag+=1;
	+ diag = new_ptr((*diag,-ndiag));
	+ channelMap_.back().push_back(ndiag);
	+ add(diag);
	+ }
	+ }
	+}
	+
	+Energy2 MEMultiChannel::scale() const {
	+ return sHat();
	+}
	+
	+int MEMultiChannel::nDim() const {
	+ return modes_[0]->nRand();
	+}
	+
	+bool MEMultiChannel::generateKinematics(const double * r) {
	+ // first find the mode
	+ int imode = -1;
	+ for(unsigned int ix=0;ix<modes_.size();++ix) {
	+ bool matched=true;
	+ for(unsigned int iy=0;iy<modes_[ix]->outgoing().size();++iy) {
	+ if(mePartonData()[iy+2]!=modes_[ix]->outgoing()[iy]) {
	+ matched=false;
	+ break;
	+ }
	+ }
	+ if(matched) {
	+ imode=ix;
	+ break;
	+ }
	+ }
	+ assert(imode>=0);
	+ // fill the stack of random numbers
	+ modes_[imode]->fillStack(r);
	+ // generate the momenta
	+ int ichan;
	+ vector<Lorentz5Momentum> momenta(meMomenta().size()-2);
	+ Lorentz5Momentum pcm=meMomenta()[0]+meMomenta()[1];
	+ Energy wgt = modes_[imode]->weight(ichan,pcm,momenta);
	+ // set the jacobian
	+ jacobian(wgt/sqrt(sHat()));
	+ // and momenta
	+ for(unsigned int ix=0;ix<momenta.size();++ix)
	+ meMomenta()[ix+2]=momenta[ix];
	+ // check the cuts
	+ tcPDVector tout(momenta.size());
	+ vector<LorentzMomentum> out(meMomenta().size()-2);
	+ for(unsigned int ix=2;ix<meMomenta().size();++ix) {
	+ tout[ix-2] = mePartonData()[ix];
	+ out[ix-2]=meMomenta()[ix];
	+ }
	+ if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) )
	+ return false;
	+ // passed cuts
	+ return true;
	+}
	+
	+CrossSection MEMultiChannel::dSigHatDR() const {
	+ return sqr(hbarc)me2()jacobian()/sHat();
	+}
	+
	+Selector<MEBase::DiagramIndex>
	+MEMultiChannel::diagrams(const DiagramVector & diags) const {
	+ Selector<DiagramIndex> sel;
	+ for ( DiagramIndex i = 0; i < diags.size(); ++i ) {
	+double wgt = me2(channelMap_[iMode_][-diags[i]->id()-1] );
	+sel.insert(wgt, i);
	+ }
	+ return sel;
	+}
	+
	+Selector<const ColourLines *>
	+MEMultiChannel::colourGeometries(tcDiagPtr ) const {
	+ static ColourLines none("");
	+ Selector<const ColourLines *> sel;
	+ sel.insert(1.0, &none);
	+ return sel;
	+}
	+
	+void MEMultiChannel::persistentOutput(PersistentOStream & os) const {
	+ os << modes_ << channelMap_;
	+}
	+
	+void MEMultiChannel::persistentInput(PersistentIStream & is, int) {
	+ is >> modes_ >> channelMap_;
	+}
	+
	+
	+// The following static variable is needed for the type
	+// description system in ThePEG.
	+DescribeAbstractClass<MEMultiChannel,MEBase>
	+ describeHerwigMEMultiChannel("Herwig::MEMultiChannel", "Herwig.so");
	+
	+void MEMultiChannel::Init() {
	+
	+ static ClassDocumentation<MEMultiChannel> documentation
	+ ("There is no documentation for the MEMultiChannel class");
	+
	+}
	+
	+
	+void MEMultiChannel::doinit() {
	+ MEBase::doinit();
	+for(tPhaseSpaceModePtr mode : modes_)
	+ mode->init();
	+}
	+
	+void MEMultiChannel::doinitrun() {
	+ MEBase::doinitrun();
	+for(tPhaseSpaceModePtr mode : modes_)
	+ mode->initrun();
	+}
	diff --git a/MatrixElement/MEMultiChannel.h b/MatrixElement/MEMultiChannel.h
	new file mode 100644
	--- /dev/null
	+++ b/MatrixElement/MEMultiChannel.h
	@@ -0,0 +1,214 @@
	+// -- C++ --
	+#ifndef Herwig_MEMultiChannel_H
	+#define Herwig_MEMultiChannel_H
	+//
	+// This is the declaration of the MEMultiChannel class.
	+//
	+
	+#include "ThePEG/MatrixElement/MEBase.h"
	+#include "Herwig/Decay/PhaseSpaceMode.fh"
	+
	+namespace Herwig {
	+
	+using namespace ThePEG;
	+
	+/**
	+ * Here is the documentation of the MEMultiChannel class.
	+ *
	+ * @see \ref MEMultiChannelInterfaces "The interfaces"
	+ * defined for MEMultiChannel.
	+ */
	+class MEMultiChannel: public MEBase {
	+
	+public:
	+
	+ /**
	+ * The default constructor.
	+ */
	+ MEMultiChannel() {};
	+
	+ /**
	+ * The destructor.
	+ */
	+ virtual ~MEMultiChannel();
	+ //@}
	+
	+public:
	+
	+ /** @name Virtual functions required by the MEBase class. */
	+ //@{
	+
	+ /**
	+ * 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.
	+ */
	+ double me2() const {return me2(-1);}
	+
	+ /**
	+ * Return the scale associated with the last set phase space point.
	+ */
	+ virtual Energy2 scale() const;
	+
	+ /**
	+ * The number of internal degrees of freedom used in the matrix
	+ * element.
	+ */
	+ virtual int nDim() const;
	+
	+ /**
	+ * Generate internal degrees of freedom given nDim() uniform
	+ * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space
	+ * generator, the dSigHatDR should be a smooth function of these
	+ * numbers, although this is not strictly necessary.
	+ * @param r a pointer to the first of nDim() consecutive random numbers.
	+ * @return true if the generation succeeded, otherwise false.
	+ */
	+ virtual bool generateKinematics(const double * r);
	+
	+ /**
	+ * Return the matrix element squared differential in the variables
	+ * given by the last call to generateKinematics().
	+ */
	+ virtual CrossSection dSigHatDR() 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;
	+ //@}
	+
	+protected :
	+
	+ /**
	+ * Add a phase-space mode
	+ */
	+ void addMode(PhaseSpaceModePtr mode) {
	+ modes_.push_back(mode);
	+ }
	+
	+ /**
	+ * Access to the modes
	+ */
	+ tPhaseSpaceModePtr mode(unsigned int imode) const {
	+ return modes_[imode];
	+ }
	+
	+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 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();
	+
	+ /**
	+ * Initialize this object. Called in the run phase just before
	+ * a run begins.
	+ */
	+ virtual void doinitrun();
	+ //@}
	+
	+protected:
	+
	+ /**
	+ * Return the matrix element squared for a given mode and phase-space channel.
	+ * @param ichan The channel we are calculating the matrix element for.
	+ * @param part The decaying Particle.
	+ * @param outgoing The particles produced in the decay
	+ * @param momenta The momenta of the particles produced in the decay
	+ * @param meopt Option for the calculation of the matrix element
	+ * @return The matrix element squared for the phase-space configuration.
	+ */
	+ virtual double me2(const int ichan) const = 0;
	+
	+ /**
	+ * Access the mode currently being used
	+ */
	+ unsigned int iMode() {return iMode_;}
	+
	+ /**
	+ * Set the mode currently begining used
	+ */
	+ void iMode(unsigned int in) const {iMode_=in;}
	+
	+private:
	+
	+ /**
	+ * The assignment operator is private and must never be called.
	+ * In fact, it should not even be implemented.
	+ */
	+ MEMultiChannel & operator=(const MEMultiChannel &);
	+
	+private:
	+
	+ /**
	+ * The phase-space modes
	+ */
	+ vector<PhaseSpaceModePtr> modes_;
	+
	+ /**
	+ * Map from the phase space channels to the modes
	+ */
	+ mutable vector<vector <int> > channelMap_;
	+
	+ /**
	+ * The mode currently begining used
	+ */
	+ mutable unsigned int iMode_;
	+};
	+
	+}
	+
	+#endif /* Herwig_MEMultiChannel_H */
	diff --git a/MatrixElement/Makefile.am b/MatrixElement/Makefile.am
	--- a/MatrixElement/Makefile.am
	+++ b/MatrixElement/Makefile.am
	@@ -1,12 +1,13 @@
	SUBDIRS = General Lepton Hadron DIS Powheg Gamma Matchbox Reweighters

	noinst_LTLIBRARIES = libHwME.la

	libHwME_la_SOURCES = \
	HwMEBase.h HwMEBase.fh HwMEBase.cc \
	+MEMultiChannel.h MEMultiChannel.cc \
	MEfftoVH.h MEfftoVH.cc \
	MEfftoffH.h MEfftoffH.cc \
	HardVertex.fh HardVertex.h HardVertex.cc \
	ProductionMatrixElement.h ProductionMatrixElement.cc \
	DrellYanBase.h DrellYanBase.cc \
	BlobME.h BlobME.cc

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 9:13 PM (1 d, 33 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806217
Default Alt Text: (32 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions