diff --git a/MatrixElement/Lepton/MEee2Mesons.cc b/MatrixElement/Lepton/MEee2Mesons.cc
--- a/MatrixElement/Lepton/MEee2Mesons.cc
+++ b/MatrixElement/Lepton/MEee2Mesons.cc
@@ -1,190 +1,190 @@
// -*- 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/DecayIntegrator.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/StandardModel/StandardModelBase.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/PDF/PolarizedBeamParticleData.h"
using namespace Herwig;
typedef LorentzVector<complex<InvEnergy> > LorentzPolarizationVectorInvE;
MEee2Mesons::MEee2Mesons() {}
Energy2 MEee2Mesons::scale() const {
return sHat();
}
unsigned int MEee2Mesons::orderInAlphaS() const {
return 0;
}
unsigned int MEee2Mesons::orderInAlphaEW() const {
- return 0;
+ return 2;
}
IBPtr MEee2Mesons::clone() const {
return new_ptr(*this);
}
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_;
+ os << current_ << modeMap_;
}
void MEee2Mesons::persistentInput(PersistentIStream & is, int) {
-is >> current_ >> modeMap_;
+ 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
int nOut = int(meMomenta().size())-2;
vector<unsigned int> constants(nOut+1);
vector<PDT::Spin > iSpin(nOut);
vector<int> hadrons(nOut);
int itemp(1);
int ix(nOut);
do {
--ix;
iSpin[ix] = mePartonData()[ix+2]->iSpin();
itemp *= iSpin[ix];
constants[ix] = itemp;
hadrons[ix] = mePartonData()[ix+2]->id();
}
while(ix>0);
constants[nOut] = 1;
// 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());
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(tcPDPtr(), IsoSpin::IUnknown, IsoSpin::I3Unknown, imode,ichan,
q,out,momenta,DecayIntegrator::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(int ix=nOut;ix>0;--ix) {
ihel[ix+1]=(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);
}
}
}
// symmetry factors
map<long,int> ncount;
double symmetry(1.);
for(tPDPtr o : out) ncount[o->id()]+=1;
for(map<long,int>::const_iterator it=ncount.begin();it!=ncount.end();++it) {
symmetry *= it->second;
}
// prefactors
output *= 0.25*sqr(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/symmetry;
}
void MEee2Mesons::constructVertex(tSubProPtr) {
}
diff --git a/MatrixElement/Lepton/MEee2Mesons.h b/MatrixElement/Lepton/MEee2Mesons.h
--- a/MatrixElement/Lepton/MEee2Mesons.h
+++ b/MatrixElement/Lepton/MEee2Mesons.h
@@ -1,148 +1,149 @@
// -*- C++ -*-
#ifndef Herwig_MEee2Mesons_H
#define Herwig_MEee2Mesons_H
//
// This is the declaration of the MEee2Mesons class.
//
#include "Herwig/MatrixElement/MEMultiChannel.h"
#include "Herwig/Decay/WeakCurrents/WeakCurrent.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
/**
- * Here is the documentation of the MEee2Mesons class.
+ * The MEee2Mesons class implements \f$e^+e^-\f$ annhilation to mesons at low energy
+ * using hadronic currents.
*
* @see \ref MEee2MesonsInterfaces "The interfaces"
* defined for MEee2Mesons.
*/
class MEee2Mesons: public MEMultiChannel {
public:
/**
* The default constructor.
*/
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;
/**
* Return the scale associated with the last set phase space point.
*/
virtual Energy2 scale() 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 &) = delete;
private :
/**
* the hadronic current
*/
WeakCurrentPtr current_;
/**
* The matrix element
*/
mutable ProductionMatrixElement me_;
/**
* Map for the modes
*/
map<int,int> modeMap_;
};
}
#endif /* Herwig_MEee2Mesons_H */
diff --git a/MatrixElement/MEDM2Mesons.cc b/MatrixElement/MEDM2Mesons.cc
new file mode 100644
--- /dev/null
+++ b/MatrixElement/MEDM2Mesons.cc
@@ -0,0 +1,194 @@
+// -*- C++ -*-
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the MEDM2Mesons class.
+//
+
+#include "MEDM2Mesons.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/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
+
+using namespace Herwig;
+typedef LorentzVector<complex<InvEnergy> > LorentzPolarizationVectorInvE;
+
+MEDM2Mesons::MEDM2Mesons() {}
+
+Energy2 MEDM2Mesons::scale() const {
+ return sHat();
+}
+
+unsigned int MEDM2Mesons::orderInAlphaS() const {
+ return 0;
+}
+
+unsigned int MEDM2Mesons::orderInAlphaEW() const {
+ return 0;
+}
+
+IBPtr MEDM2Mesons::clone() const {
+ return new_ptr(*this);
+}
+
+IBPtr MEDM2Mesons::fullclone() const {
+ return new_ptr(*this);
+}
+
+void MEDM2Mesons::doinit() {
+ // make sure the current got initialised
+ current_->init();
+ // max energy
+ Energy Emax = generator()->maximumCMEnergy();
+ // 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(incomingA_,out,1.,
+ incomingB_,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 MEDM2Mesons::persistentOutput(PersistentOStream & os) const {
+ os << current_ << modeMap_ << incomingA_ << incomingB_;
+}
+
+void MEDM2Mesons::persistentInput(PersistentIStream & is, int) {
+ is >> current_ >> modeMap_ >> incomingA_ >> incomingB_;
+}
+
+//The following static variable is needed for the type
+// description system in ThePEG.
+DescribeClass<MEDM2Mesons,MEMultiChannel>
+ describeHerwigMEDM2Mesons("Herwig::MEDM2Mesons", "Herwig.so");
+
+void MEDM2Mesons::Init() {
+
+ static ClassDocumentation<MEDM2Mesons> documentation
+ ("The MEDM2Mesons class simulates the annhilation of"
+ " DM particles to mesons at low energy");
+
+ static Reference<MEDM2Mesons,WeakCurrent> interfaceWeakCurrent
+ ("WeakCurrent",
+ "The reference for the decay current to be used.",
+ &MEDM2Mesons::current_, false, false, true, false, false);
+
+ static Reference<MEDM2Mesons,ParticleData> interfaceIncomingA
+ ("IncomingA",
+ "First incoming particle",
+ &MEDM2Mesons::incomingA_, false, false, true, false, false);
+
+ static Reference<MEDM2Mesons,ParticleData> interfaceIncomingB
+ ("IncomingB",
+ "Second incoming particle",
+ &MEDM2Mesons::incomingB_, false, false, true, false, false);
+
+}
+
+
+double MEDM2Mesons::me2(const int ichan) const {
+ // compute the incoming current
+ LorentzPolarizationVectorInvE lepton[2][2];
+ if(incomingA_->iSpin()==PDT::Spin1Half && incomingB_->iSpin()==PDT::Spin1Half) {
+ 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);
+ }
+ // this should be coupling of DM to mediator/ mediator propagator
+ InvEnergy2 pre = 1./GeV2;
+ 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());
+ }
+ }
+ }
+ // TODO think about other spins for the DM
+ else
+ assert(false);
+ // work out the mapping for the hadron vector
+ int nOut = int(meMomenta().size())-2;
+ vector<unsigned int> constants(nOut+1);
+ vector<PDT::Spin > iSpin(nOut);
+ vector<int> hadrons(nOut);
+ int itemp(1);
+ int ix(nOut);
+ do {
+ --ix;
+ iSpin[ix] = mePartonData()[ix+2]->iSpin();
+ itemp *= iSpin[ix];
+ constants[ix] = itemp;
+ hadrons[ix] = mePartonData()[ix+2]->id();
+ }
+ while(ix>0);
+ constants[nOut] = 1;
+ // 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());
+ vector<Lorentz5Momentum> momenta(meMomenta() .begin()+2, meMomenta().end());
+ tPDVector out = mode(modeMap_.at(imode))->outgoing();
+ if(ichan<0) iMode(modeMap_.at(imode));
+ // get the hadronic currents for the I=1 and I=0 components
+ vector<LorentzPolarizationVectorE>
+ hadronI0(current_->current(tcPDPtr(), IsoSpin::IZero, IsoSpin::I3Zero, imode,ichan,
+ q,out,momenta,DecayIntegrator::Calculate));
+ vector<LorentzPolarizationVectorE>
+ hadronI1(current_->current(tcPDPtr(), IsoSpin::IOne, IsoSpin::I3Zero, imode,ichan,
+ q,out,momenta,DecayIntegrator::Calculate));
+ // compute the matrix element
+ vector<unsigned int> ihel(meMomenta().size());
+ double output(0.);
+ for(unsigned int hhel=0;hhel<hadronI0.size();++hhel) {
+ // map the index for the hadrons to a helicity state
+ for(int ix=nOut;ix>0;--ix) {
+ ihel[ix+1]=(hhel%constants[ix-1])/constants[ix];
+ }
+ // loop over the helicities of the incoming particles
+ for(ihel[1]=0;ihel[1]<2;++ihel[1]){
+ for(ihel[0]=0;ihel[0]<2;++ihel[0]) {
+ // work one coefficients for the I1 and I0 bits
+ Complex amp = lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel]+hadronI1[hhel]);
+ me_(ihel)= amp;
+ output += std::norm(amp);
+ }
+ }
+ }
+ // symmetry factors
+ map<long,int> ncount;
+ double symmetry(1.);
+ for(tPDPtr o : out) ncount[o->id()]+=1;
+ for(map<long,int>::const_iterator it=ncount.begin();it!=ncount.end();++it) {
+ symmetry *= it->second;
+ }
+ // prefactors
+ output *= 0.25*sqr(pow(sqrt(sHat())/q,int(momenta.size()-2)));
+ return output/symmetry;
+}
+
+
+void MEDM2Mesons::constructVertex(tSubProPtr) {
+}
diff --git a/MatrixElement/MEDM2Mesons.h b/MatrixElement/MEDM2Mesons.h
new file mode 100644
--- /dev/null
+++ b/MatrixElement/MEDM2Mesons.h
@@ -0,0 +1,162 @@
+// -*- C++ -*-
+#ifndef Herwig_MEDM2Mesons_H
+#define Herwig_MEDM2Mesons_H
+//
+// This is the declaration of the MEDM2Mesons class.
+//
+
+#include "Herwig/MatrixElement/MEMultiChannel.h"
+#include "Herwig/Decay/WeakCurrents/WeakCurrent.h"
+#include "Herwig/MatrixElement/ProductionMatrixElement.h"
+
+namespace Herwig {
+
+using namespace ThePEG;
+
+/**
+ * The MEDM2Mesons class implements dark matter annhilation via a vector current to
+ * mesons at low energies
+ *
+ * @see \ref MEDM2MesonsInterfaces "The interfaces"
+ * defined for MEDM2Mesons.
+ */
+class MEDM2Mesons: public MEMultiChannel {
+
+public:
+
+ /**
+ * The default constructor.
+ */
+ MEDM2Mesons();
+
+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;
+
+ /**
+ * Return the scale associated with the last set phase space point.
+ */
+ virtual Energy2 scale() 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;
+ //@}
+
+private:
+
+ /**
+ * 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.
+ */
+ MEDM2Mesons & operator=(const MEDM2Mesons &);
+
+private :
+
+ /**
+ * Hadronic current etc
+ */
+ //@{
+ /**
+ * the hadronic current
+ */
+ WeakCurrentPtr current_;
+
+ /**
+ * The matrix element
+ */
+ mutable ProductionMatrixElement me_;
+
+ /**
+ * Map for the modes
+ */
+ map<int,int> modeMap_;
+ //@}
+
+ /**
+ * DM
+ */
+ //@{
+ /**
+ * Incoming Particles
+ */
+ PDPtr incomingA_, incomingB_;
+ //@}
+};
+
+}
+
+#endif /* Herwig_MEDM2Mesons_H */
diff --git a/MatrixElement/Makefile.am b/MatrixElement/Makefile.am
--- a/MatrixElement/Makefile.am
+++ b/MatrixElement/Makefile.am
@@ -1,13 +1,14 @@
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
+BlobME.h BlobME.cc \
+MEDM2Mesons.h MEDM2Mesons.cc