diff --git a/Decay/ScalarMeson/EtaPiPiFermionsDecayer.cc b/Decay/ScalarMeson/EtaPiPiFermionsDecayer.cc --- a/Decay/ScalarMeson/EtaPiPiFermionsDecayer.cc +++ b/Decay/ScalarMeson/EtaPiPiFermionsDecayer.cc @@ -1,359 +1,354 @@ // -*- C++ -*- // // EtaPiPiFermionsDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 EtaPiPiFermionsDecayer class. // #include "EtaPiPiFermionsDecayer.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Parameter.h" -#include "ThePEG/Interface/ParVector.h" +#include "ThePEG/Interface/Command.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/Helicity/epsilon.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" #include "ThePEG/Helicity/HelicityFunctions.h" #include "Herwig/Decay/FormFactors/OmnesFunction.h" using namespace Herwig; using namespace ThePEG::Helicity; DescribeClass describeHerwigEtaPiPiFermionsDecayer("Herwig::EtaPiPiFermionsDecayer", "HwSMDecay.so"); HERWIG_INTERPOLATOR_CLASSDESC(EtaPiPiFermionsDecayer,double,Energy) void EtaPiPiFermionsDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { - for(unsigned int ix=0;ixmaxWeight(); + } } } EtaPiPiFermionsDecayer::EtaPiPiFermionsDecayer() - : incoming_(2), lepton_(2), coupling_(2), maxWeight_(2), option_(2) { - // the pion decay constant - fPi_=130.7*MeV; - // the rho mass - mRho_=0.7711*GeV; - rhoWidth_=0.1492*GeV; + : fPi_(130.7*MeV), mRho_(0.7711*GeV), rhoWidth_(0.1492*GeV) { // the constants for the omnes function form aConst_=0.5/mRho_/mRho_; cConst_=1.0; // use local values of the parameters localParameters_=true; // the modes - // eta decay - incoming_[0] = 221; - option_[0] = 3; - coupling_[0] = 5.060e-3; - lepton_[0]=11; - maxWeight_[0] = 3.95072; - // eta' decay - incoming_[1] = 331; - option_[1] = 3; - coupling_[1] = 4.278e-3; - lepton_[1]=11; - maxWeight_[1] = 3.53141; - rhoConst_=0.; mPi_=ZERO; // intermediates generateIntermediates(false); } void EtaPiPiFermionsDecayer::doinit() { DecayIntegrator::doinit(); // check the consistence of the parameters unsigned int isize=incoming_.size(); if(isize!=coupling_.size()||isize!=option_.size()||isize!=lepton_.size()||isize!=maxWeight_.size()) throw InitException() << "Inconsistent parameters in " << "EtaPiPiFermionsDecayer::doinit()" << Exception::abortnow; // set the parameters tPDPtr rho(getParticleData(ParticleID::rho0)); if(!localParameters_) { mRho_=rho->mass(); rhoWidth_=rho->width(); } mPi_=getParticleData(ParticleID::piplus)->mass(); Energy pcm(Kinematics::pstarTwoBodyDecay(mRho_,mPi_,mPi_)); rhoConst_=sqr(mRho_)*rhoWidth_/pow<3,1>(pcm); // set up the modes tPDPtr gamma = getParticleData(ParticleID::gamma); PhaseSpaceModePtr mode; for(unsigned int ix=0;ixaddChannel(newChannel); addMode(mode); } } int EtaPiPiFermionsDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const { int imode(-1); // check number of external particles if(children.size()!=4){return imode;} // check the outgoing particles unsigned int npip(0),npim(0),nl(0); int il(0); tPDVector::const_iterator pit = children.begin(); int id; for(;pit!=children.end();++pit) { id=(**pit).id(); if(id==ParticleID::piplus) ++npip; else if(id==ParticleID::piminus) ++npim; else { il = abs(id); ++nl; } } if(!(npip==1&&npim==1&&nl==2)) return imode; unsigned int ix(0); id=parent->id(); do { if(id==incoming_[ix] && il==lepton_[ix]) imode=ix; ++ix; } while(imode<0&&ix> iunit(fPi_,MeV) >> incoming_ >> coupling_ >> maxWeight_ >> lepton_ >> option_ >> iunit(aConst_,1/MeV2) >> cConst_ >>iunit(mRho_,MeV) >> iunit(rhoWidth_,MeV) >> rhoConst_ >> iunit(mPi_,MeV) >> localParameters_ >> omnesFunction_; } void EtaPiPiFermionsDecayer::Init() { static ClassDocumentation documentation ("The EtaPiPiFermionsDecayer class is design for the decay of" " the eta and eta prime to pi+pi-gamma", "The decays of $\\eta,\\eta'\\to\\pi^+\\pi^-\\gamma$ were simulated" " using the matrix elements from \\cite{Venugopal:1998fq,Holstein:2001bt}", "\\bibitem{Venugopal:1998fq} E.~P.~Venugopal and B.~R.~Holstein,\n" "Phys.\\ Rev.\\ D {\\bf 57} (1998) 4397 [arXiv:hep-ph/9710382].\n" "%%CITATION = PHRVA,D57,4397;%%\n" "\\bibitem{Holstein:2001bt} B.~R.~Holstein,\n" " Phys.\\ Scripta {\\bf T99} (2002) 55 [arXiv:hep-ph/0112150].\n" "%%CITATION = PHSTB,T99,55;%%\n"); static Reference interfaceOmnesFunction ("OmnesFunction", "Omnes function for the matrix element", &EtaPiPiFermionsDecayer::omnesFunction_, false, false, true, false, false); static Parameter interfacefpi ("fpi", "The pion decay constant", &EtaPiPiFermionsDecayer::fPi_, MeV, 130.7*MeV, ZERO, 200.*MeV, false, false, false); - static ParVector interfaceIncoming - ("Incoming", - "The PDG code for the incoming particle", - &EtaPiPiFermionsDecayer::incoming_, - 0, 0, 0, -10000000, 10000000, false, false, true); - - static ParVector interfaceCoupling - ("Coupling", - "The coupling for the decay mode", - &EtaPiPiFermionsDecayer::coupling_, - 0, 0, 0, 0., 100., false, false, true); - - static ParVector interfaceMaxWeight - ("MaxWeight", - "The maximum weight for the decay mode", - &EtaPiPiFermionsDecayer::maxWeight_, - 0, 0, 0, 0., 100., false, false, true); + static Command interfaceSetUpDecayMode + ("SetUpDecayMode", + "Set up the particles (incoming, lepton, option for VMD, coupling and maxweight." + "There are three options for for VMD, 0 is a VMD model using M Gamma for the width term," + " 1 is a VMD model using q Gamma for the width term," + "2. analytic form of the Omnes function," + "3. experimental form of the Omnes function.", + &EtaPiPiFermionsDecayer::setUpDecayMode, false); static Parameter interfaceRhoMass ("RhoMass", "The mass of the rho", &EtaPiPiFermionsDecayer::mRho_, MeV, 771.1*MeV, 400.*MeV, 1000.*MeV, false, false, false); static Parameter interfaceRhoWidth ("RhoWidth", "The width of the rho", &EtaPiPiFermionsDecayer::rhoWidth_, MeV, 149.2*MeV, 100.*MeV, 300.*MeV, false, false, false); static Switch interfaceLocalParameters ("LocalParameters", "Use local values of the rho mass and width", &EtaPiPiFermionsDecayer::localParameters_, true, false, false); static SwitchOption interfaceLocalParametersLocal (interfaceLocalParameters, "Local", "Use local parameters", true); static SwitchOption interfaceLocalParametersParticleData (interfaceLocalParameters, "ParticleData", "Use values from the particle data objects", false); static Parameter interfaceOmnesC ("OmnesC", "The constant c for the Omnes form of the prefactor", &EtaPiPiFermionsDecayer::cConst_, 1.0, -10., 10., false, false, false); static Parameter interfaceOmnesA ("OmnesA", "The constant a for the Omnes form of the prefactor", &EtaPiPiFermionsDecayer::aConst_, 1./GeV2, 0.8409082/GeV2, ZERO, 10./GeV2, false, false, false); - static ParVector interfaceOption - ("Option", - "The form of the prefactor 0 is a VMD model using M Gamma for the width term," - "1 is a VMD model using q Gamma for the width term," - "2. analytic form of the Omnes function," - "3. experimental form of the Omnes function.", - &EtaPiPiFermionsDecayer::option_, - 0, 0, 0, 0, 4, false, false, true); - - static ParVector interfaceLepton - ("Lepton", - "The type of lepton", - &EtaPiPiFermionsDecayer::lepton_, -1, 11, 11, 15, - false, false, Interface::limited); +} -} void EtaPiPiFermionsDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // set up the spin information for the decay products ScalarWaveFunction::constructSpinInfo(const_ptr_cast(&part), incoming,true); for(unsigned int ix=0;ix<2;++ix) ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true); SpinorBarWaveFunction:: constructSpinInfo(wavebar_,decay[2],outgoing,true); SpinorWaveFunction:: constructSpinInfo(wave_ ,decay[3],outgoing,true); } double EtaPiPiFermionsDecayer::me2(const int,const Particle & part, const tPDVector &, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin0, PDT::Spin1Half,PDT::Spin1Half))); useMe(); if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(rho_,const_ptr_cast(&part),incoming); } // calculate the spinors wavebar_.resize(2); wave_ .resize(2); for(unsigned int ihel=0;ihel<2;++ihel) { wavebar_[ihel] = HelicityFunctions::dimensionedSpinorBar(-momenta[2],ihel,Helicity::outgoing); wave_ [ihel] = HelicityFunctions::dimensionedSpinor (-momenta[3],ihel,Helicity::outgoing); } Lorentz5Momentum pff(momenta[2]+momenta[3]); pff.rescaleMass(); Energy2 mff2(pff.mass()*pff.mass()); // prefactor for the matrix element complex pre(part.mass()*coupling_[imode()]*2.*sqrt(2.)/(fPi_*fPi_*fPi_)* sqrt(SM().alphaEM(mff2)*4.*Constants::pi)/mff2); Lorentz5Momentum ppipi(momenta[0]+momenta[1]); ppipi.rescaleMass(); Energy q(ppipi.mass()); Energy2 q2(q*q); Complex ii(0.,1.); // first VMD option Complex fact; if(option_[imode()]==0) { Energy pcm(Kinematics::pstarTwoBodyDecay(q,mPi_,mPi_)); Complex resfact(q2/(mRho_*mRho_-q2-ii*mRho_*pcm*pcm*pcm*rhoConst_/q2)); fact=(1.+1.5*resfact); } // second VMD option else if(option_[imode()]==1) { Energy pcm(Kinematics::pstarTwoBodyDecay(q,mPi_,mPi_)); Complex resfact(q2/(mRho_*mRho_-q2-ii*pcm*pcm*pcm*rhoConst_/q)); fact=(1.+1.5*resfact); } // omnes function else if(option_[imode()]==2 || option_[imode()]==3) { fact=(1.-cConst_+cConst_*(1.+aConst_*q2)/omnesFunction_->D(q2)); } pre = pre*fact; LorentzVector > epstemp(pre*Helicity::epsilon(momenta[0],momenta[1],pff)); // compute the matrix element vector ispin(5,0); for(ispin[3]=0;ispin[3]<2;++ispin[3]) { for(ispin[4]=0;ispin[4]<2;++ispin[4]) { LorentzPolarizationVectorE fcurrent = wave_[ispin[4]].vectorCurrent(wavebar_[ispin[3]]); (*ME())(ispin) = Complex(epstemp.dot(fcurrent)); } } // contract the whole thing return ME()->contract(rho_).real(); } void EtaPiPiFermionsDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); output << "newdef " << name() << ":fpi " << fPi_/MeV << "\n"; output << "newdef " << name() << ":RhoMass " << mRho_/MeV << "\n"; output << "newdef " << name() << ":RhoWidth " << rhoWidth_/MeV << "\n"; output << "newdef " << name() << ":LocalParameters " << localParameters_ << "\n"; output << "newdef " << name() << ":OmnesC " << cConst_ << "\n"; output << "newdef " << name() << ":OmnesA " << aConst_*GeV2 << "\n"; - for(unsigned int ix=0;ix<2;++ix) { - output << "newdef " << name() << ":Incoming " << ix << " " - << incoming_[ix] << "\n"; - output << "newdef " << name() << ":Coupling " << ix << " " - << coupling_[ix] << "\n"; - output << "newdef " << name() << ":Lepton " << ix << " " - << lepton_[ix] << "\n"; - output << "newdef " << name() << ":MaxWeight " << ix << " " - << maxWeight_[ix] << "\n"; - output << "newdef " << name() << ":Option " << ix << " " - << option_[ix] << "\n"; + for(unsigned int ix=0;ixdataBaseOutput(output,false,true); output << "newdef " << name() << ":OmnesFunction " << omnesFunction_->name() << " \n"; if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } + +string EtaPiPiFermionsDecayer::setUpDecayMode(string arg) { + // parse first bit of the string + string stype = StringUtils::car(arg); + arg = StringUtils::cdr(arg); + // extract PDG code for the incoming particle + long in = stoi(stype); + tcPDPtr pData = getParticleData(in); + if(!pData) + return "Incoming particle with id " + std::to_string(in) + "does not exist"; + if(pData->iSpin()!=PDT::Spin0) + return "Incoming particle with id " + std::to_string(in) + "does not have spin 0"; + // and outgoing particles + stype = StringUtils::car(arg); + arg = StringUtils::cdr(arg); + int out = stoi(stype); + pData = getParticleData(out); + if(!pData) + return "Outgoing fermion with id " + std::to_string(out) + "does not exist"; + if(pData->iSpin()!=PDT::Spin1Half) + return "Outgoing fermion with id " + std::to_string(out) + "does not have spin 1/2"; + // vmd option + stype = StringUtils::car(arg); + arg = StringUtils::cdr(arg); + int VMDinclude = stoi(stype); + // get the coupling + stype = StringUtils::car(arg); + arg = StringUtils::cdr(arg); + double g = stof(stype); + // and the maximum weight + stype = StringUtils::car(arg); + arg = StringUtils::cdr(arg); + double wgt = stof(stype); + // store the information + incoming_.push_back(in); + lepton_.push_back(out); + coupling_.push_back(g); + option_.push_back(VMDinclude); + maxWeight_.push_back(wgt); + // success + return ""; +} diff --git a/Decay/ScalarMeson/EtaPiPiFermionsDecayer.h b/Decay/ScalarMeson/EtaPiPiFermionsDecayer.h --- a/Decay/ScalarMeson/EtaPiPiFermionsDecayer.h +++ b/Decay/ScalarMeson/EtaPiPiFermionsDecayer.h @@ -1,255 +1,262 @@ // -*- C++ -*- // // EtaPiPiFermionsDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_EtaPiPiFermionsDecayer_H #define HERWIG_EtaPiPiFermionsDecayer_H // This is the declaration of the EtaPiPiFermionsDecayer class. #include "Herwig/Utilities/Kinematics.h" #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/PhaseSpaceMode.h" #include "Herwig/Decay/FormFactors/OmnesFunction.fh" #include "ThePEG/Helicity/LorentzSpinorBar.h" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * * The EtaPiPiFermionsDecayer class implements the decay of * the \f$\eta\f$ or \f$\eta'\f$ to \f$\pi^+\pi^-\gamma\f$ using either * a VMD type model or a model using either the theoretical or experimental * form of the Omnes function taken from hep-ph/0112150. * * The matrix element is given by * \f[\mathcal{M} = B(s_{+-},s_{+\gamma},s_{-\gamma})\epsilon^{\mu\nu\alpha\beta} * \epsilon^*_{\mu}p_{+\nu}p_{-\alpha}p_{\gamma\beta}\f] * where \f$p_{+,-}\f$ are the momenta of the positively and negatively charged pions, * \f$p_{\gamma}\f$ is the momentum of the photon and \f$s_{ij} = (p_i+p_j)^2\f$. * * The different models take * * \f[B(s_{+-},s_{+\gamma},s_{-\gamma}) = * B_0\left(1+\frac32\frac{s_{+-}}{M^2_\rho-s_{+-}-iM_\rho\Gamma_\rho(s_{+-})}\right)\f] * where \f$M_\rho\f$ and \f$\Gamma_\rho\f$ are the mass and running width * of the \f$\rho\f$ * respectively for the VMD model. * * For the Omnes function case we take * * \f[B(s_{+-},s_{+\gamma},s_{-\gamma}) = * B_0\left(1-c+c\frac{1+as_{+-}}{D_1(s_{+-})}\right)\f] * either the experimental or analytic form of the Omnes function \f$D_1(s_{+-})\f$ * taken from hep-ph/0112150 can be used. * * The coefficient \f$B_0\f$ is given in hep-ph/0112150. We use the values from this * paper and use their default choice \f$c=1\f$, \f$a=\frac1{2M_\rho}\f$. * * @see DecayIntegrator * */ class EtaPiPiFermionsDecayer: public DecayIntegrator { public: /** * Default constructor. */ EtaPiPiFermionsDecayer(); /** * Which of the possible decays is required * @param cc Is this mode the charge conjugate * @param parent The decaying particle * @param children The decay products */ virtual int modeNumber(bool & cc, tcPDPtr parent, const tPDVector & children) const; /** * 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. */ double me2(const int ichan,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const; /** * Construct the SpinInfos for the particles produced in the decay */ virtual void constructSpinInfo(const Particle & part, ParticleVector outgoing) const; /** * Output the setup information for the particle database * @param os The stream to output the information to * @param header Whether or not to output the information for MySQL */ virtual void dataBaseOutput(ofstream & os,bool header) 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); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ 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. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving and * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object to the begining of the run phase. */ virtual void doinitrun(); //@} +public: + + /** + * Set the parameters for a decay mode + */ + string setUpDecayMode(string arg); + private: /** * Private and non-existent assignment operator. */ EtaPiPiFermionsDecayer & operator=(const EtaPiPiFermionsDecayer &) = delete; private: /** * the pion decay constant, \f$F_\pi\f$. */ Energy fPi_; /** * the PDG code for the incoming particle */ vector incoming_; /** * The PDG code of the leptons */ vector lepton_; /** * Coupling for the decay, \f$B_0\f$. */ vector coupling_; /** * The maximum weight */ vector maxWeight_; /** * The option for the energy dependence of the prefactor */ vector option_; /** * The constants for the omnes function form. */ InvEnergy2 aConst_; /** * The constants for the Omnes function form. */ double cConst_; /** * The \f$\rho\f$ mass */ Energy mRho_; /** * The \f$\rho\f$ width */ Energy rhoWidth_; /** * Constant for the running \f$rho\f$ width. */ double rhoConst_; /** * The \f$m_\pi\f$. */ Energy mPi_; /** * Use local values of the parameters. */ bool localParameters_; /** * Object calculating the Omnes function */ OmnesFunctionPtr omnesFunction_; /** * Spin densit matrix */ mutable RhoDMatrix rho_; /** * Spinors for the fermions */ mutable vector > wave_; /** * Barred spinors for the fermions */ mutable vector > wavebar_; }; } #endif /* HERWIG_EtaPiPiFermionsDecayer_H */