diff --git a/Decay/General/DMMediatorDecayer.cc b/Decay/General/DMMediatorDecayer.cc --- a/Decay/General/DMMediatorDecayer.cc +++ b/Decay/General/DMMediatorDecayer.cc @@ -1,238 +1,273 @@ // -*- C++ -*- // // DMMediatorDecayer.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 DMMediatorDecayer class. // #include "DMMediatorDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" +#include "Herwig/Models/General/BSMModel.h" using namespace Herwig; IBPtr DMMediatorDecayer::clone() const { return new_ptr(*this); } IBPtr DMMediatorDecayer::fullclone() const { return new_ptr(*this); } void DMMediatorDecayer::persistentOutput(PersistentOStream & os) const { - os << inpart_ << currentOut_ << current_ << mode_ << wgtloc_ << wgtmax_ << weights_ << cDMmed_ << cSMmed_; + os << inpart_ << currentOut_ << current_ << mode_ << wgtloc_ << wgtmax_ << weights_ << cSMmed_; } void DMMediatorDecayer::persistentInput(PersistentIStream & is, int) { - is >> inpart_ >> currentOut_ >> current_ >> mode_ >> wgtloc_ >> wgtmax_ >> weights_ >> cDMmed_ >> cSMmed_; + is >> inpart_ >> currentOut_ >> current_ >> mode_ >> wgtloc_ >> wgtmax_ >> weights_ >> cSMmed_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigDMMediatorDecayer("Herwig::DMMediatorDecayer", "Herwig.so"); void DMMediatorDecayer::Init() { static ClassDocumentation documentation ("The DMMediatorDecayer class is designed for the decays of low mass vector bosons"); } -void DMMediatorDecayer::setDecayInfo(PDPtr in, const vector & outCurrent, - WeakCurrentPtr current, double normin, vector sm) { +void DMMediatorDecayer::setDecayInfo(PDPtr in, const vector & outCurrent, WeakCurrentPtr current) { inpart_ = in; currentOut_ = outCurrent; current_ = current; - cDMmed_=normin; - cSMmed_=sm; + // cast the model + Ptr::ptr model = dynamic_ptr_cast::ptr>(generator()->standardModel()); + bool foundU(false),foundD(false),foundS(false); + // find the vertices we need and extract the couplings + for(unsigned int ix = 0; ix < model->numberOfVertices(); ++ix ) { + VertexBasePtr vertex = model->vertex(ix); + if(vertex->getNpoint()!=3) continue; + for(unsigned int iloc = 0;iloc < 3; ++iloc) { + vector ext = vertex->search(iloc, in->id()); + if(ext.empty()) continue; + for(unsigned int ioff=0;ioffsetCoupling(sqr(in->mass()),getParticleData(1),getParticleData(-1),in); + cSMmed_[0] = vertex->norm(); + } + else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) { + foundU = true; + vertex->setCoupling(sqr(in->mass()),getParticleData(2),getParticleData(-2),in); + cSMmed_[1] = vertex->norm(); + } + else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) { + foundS = true; + vertex->setCoupling(sqr(in->mass()),getParticleData(3),getParticleData(-3),in); + cSMmed_[2] = vertex->norm(); + } + } + } + } + if(!foundD) { + throw InitException() << "Cannot find down quark coupling in DMMediatorDecayer::doinit()"; + } + if(!foundU) { + throw InitException() << "Cannot find up quark coupling in DMMediatorDecayer::doinit()"; + } + if(!foundS) { + throw InitException() << "Cannot find strange quark coupling in DMMediatorDecayer::doinit()"; + } } - - int DMMediatorDecayer::modeNumber(bool & cc, tcPDPtr parent, const tPDVector & children) const { vector id; id.push_back(parent->id()); for(unsigned int ix=0;ixid()); return modeNumber(cc,id); } void DMMediatorDecayer::doinitrun() { current_->initrun(); DecayIntegrator::doinitrun(); } void DMMediatorDecayer::doinit() { DecayIntegrator::doinit(); // make sure the current got initialised current_->init(); // find the mode for(unsigned int ix=0;ixnumberOfModes();++ix) { // get the external particles for this mode int iq(0),ia(0); tPDVector ptemp = current_->particles(inpart_->iCharge(),ix,iq,ia); // check this is the right mode if(ptemp.size()!=currentOut_.size()) continue; vector matched(ptemp.size(),false); bool match = true; for(unsigned int iy=0;iycreateMode(inpart_->iCharge(),tcPDPtr(),FlavourInfo(), ix,mode,0,-1,channel,inpart_->mass()); if(done) { // the maximum weight and the channel weights // the weights for the channel if(weights_.empty()) { weights_.resize(mode->channels().size(),1./(mode->channels().size())); } mode_ = ix; // special for the two body modes if(out.size()==2) { weights_.clear(); mode=new_ptr(PhaseSpaceMode(inpart_,out,1.)); } mode->maxWeight(wgtmax_); mode->setWeights(weights_); addMode(mode); } break; } } int DMMediatorDecayer::modeNumber(bool & cc, vector id) const { // incoming particle long idtemp; tPDPtr p0=getParticleData(id[0]); idtemp = p0->CC() ? -id[0] : id[0]; if( id[0] ==inpart_->id()) cc=false; else if(idtemp==inpart_->id()) cc=true ; else return -1; vector idout; for(vector::iterator it=++id.begin();it!=id.end();++it) { idout.push_back(*it); } unsigned int icurr=current_->decayMode(idout); if(mode_==icurr) return 0; else return -1; } void DMMediatorDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast(&part), Helicity::incoming,true,false); weakCurrent()->constructSpinInfo(ParticleVector(decay.begin(),decay.end())); } double DMMediatorDecayer::me2(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { using namespace ThePEG::Helicity; // polarization vectors for the incoming particle if(meopt==Initialize) { VectorWaveFunction::calculateWaveFunctions(vectors_,rho_, const_ptr_cast(&part), incoming,false); // fix rho if no correlations fixRho(rho_); } // work out the mapping for the hadron vector int nOut = momenta.size(); vector constants(nOut+1); vector iSpin(nOut); vector hadrons(nOut); int itemp(1); int ix(nOut); do { --ix; iSpin[ix] = outgoing[ix]->iSpin(); itemp *= iSpin[ix]; constants[ix] = itemp; hadrons[ix] = outgoing[ix]->id(); } while(ix>0); constants[nOut] = 1; Energy2 scale(sqr(part.mass())); // calculate the hadron current Energy q = part.mass(); // currents for the different flavour components vector hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); vector hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); vector hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar), mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate)); // compute the matrix element GeneralDecayMEPtr newME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,iSpin))); vector ihel(momenta.size()+1); unsigned int hI0_size = hadronI0.size(); unsigned int hI1_size = hadronI1.size(); unsigned int hss_size = hadronssbar.size(); unsigned int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size); for(unsigned int hhel=0;hhel0;--ix) { ihel[ix]=(hhel%constants[ix-1])/constants[ix]; } for(ihel[0]=0;ihel[0]<3;++ihel[0]) { Complex amp = 0.; // work on coefficients for the I1 and I0 bits if(hI0_size != 0 ) amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI0[hhel])); if(hI1_size !=0) amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI1[hhel])); if(hss_size !=0) amp += cSMmed_[2] /q*(vectors_[ihel[0]].wave().dot(hadronssbar[hhel])); - (*newME)(ihel) = cDMmed_*amp; + (*newME)(ihel) = amp; } } // store the matrix element ME(newME); // return the answer double output = (ME()->contract(rho_)).real(); return output; } Energy DMMediatorDecayer::partialWidth(tPDPtr part, vector out) { vector id; id.push_back(part->id()); for(unsigned int ix=0;ixid()); bool cc; int mode=modeNumber(cc,id); imode(mode); return initializePhaseSpaceMode(mode,true,true); } diff --git a/Decay/General/DMMediatorDecayer.h b/Decay/General/DMMediatorDecayer.h --- a/Decay/General/DMMediatorDecayer.h +++ b/Decay/General/DMMediatorDecayer.h @@ -1,238 +1,238 @@ // -*- C++ -*- // // DMMediatorDecayer.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_DMMediatorDecayer_H #define HERWIG_DMMediatorDecayer_H // // This is the declaration of the DMMediatorDecayer class. // #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/Decay/DecayIntegrator.h" #include "Herwig/Decay/WeakCurrents/WeakCurrent.h" #include "Herwig/Decay/PhaseSpaceMode.h" #include "DMMediatorDecayer.fh" namespace Herwig { using namespace ThePEG; using ThePEG::Helicity::VectorWaveFunction; /** * Here is the documentation of the DMMediatorDecayer class. * * @see \ref DMMediatorDecayerInterfaces "The interfaces" * defined for DMMediatorDecayer. */ class DMMediatorDecayer: public DecayIntegrator { public: /** * The default constructor. */ - DMMediatorDecayer() : wgtmax_(0.) {} + DMMediatorDecayer() : cSMmed_({0.,0.,0.}), wgtmax_(0.) + {} /** @name Virtual functions required by the Decayer class. */ //@{ /** * 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; //@} /** @name Virtual functions required by the Decayer class. */ //@{ /** * 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; /** * Function to return partial Width * @param inpart Pointer to incoming particle data object * @param out The outgoing particles from the current */ virtual Energy partialWidth(tPDPtr inpart, vector out); //@} /** * set up the decay */ void setDecayInfo(PDPtr in, const vector & outCurrent, - WeakCurrentPtr current, - double normin, vector sm); + WeakCurrentPtr current); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ 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: /** @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: /** * The number of the mode * @param cc Whether of not this is the charge conjugate of the defined mode * @param id The PDG codes of the particles */ int modeNumber(bool & cc, vector id) const; /** * Access to the map between the number of the mode and the modes in * the current */ unsigned int mode() const { return mode_; } /** * Access to the weak current */ WeakCurrentPtr weakCurrent() const { return current_; } private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DMMediatorDecayer & operator=(const DMMediatorDecayer &) = delete; private: /** * Incoming particle **/ PDPtr inpart_; /** * Outgoing particles from the current */ vector currentOut_; /** * DM coupling to the dark mediator */ - double cDMmed_; + Complex cDMmed_; /** * SM couplings to the dark mediator */ - vector cSMmed_; + vector cSMmed_; /** * Pointer to the current */ WeakCurrentPtr current_; /** * mapping of the modes to the currents */ unsigned int mode_; /** * location of the weights */ int wgtloc_; /** * the maximum weight */ double wgtmax_; /** * The weights for the different channels */ vector weights_; /** * Spin density matrix */ mutable RhoDMatrix rho_; /** * Polarization vectors for the decaying particle */ mutable vector vectors_; }; } #endif /* HERWIG_DMMediatorDecayer_H */ diff --git a/Models/DarkMatter/MEDM2Mesons.cc b/Models/DarkMatter/MEDM2Mesons.cc --- a/Models/DarkMatter/MEDM2Mesons.cc +++ b/Models/DarkMatter/MEDM2Mesons.cc @@ -1,268 +1,266 @@ // -*- 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/Vertex/Vector/FFVVertex.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "Herwig/Models/General/BSMModel.h" using namespace Herwig; typedef LorentzVector > LorentzPolarizationVectorInvE; MEDM2Mesons::MEDM2Mesons() : cDMmed_(0.), cSMmed_({0.0,0.0,0.0}) {} 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;imodenumberOfModes();++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(), FlavourInfo(), imode,mode,0,-1,channel,Emax)) continue; modeMap_[imode] = nmode; addMode(mode); ++nmode; } // cast the model Ptr::ptr model = dynamic_ptr_cast::ptr>(generator()->standardModel()); bool foundDM(false),foundU(false),foundD(false),foundS(false); // find the vertices we need and extract the couplings for(unsigned int ix = 0; ix < model->numberOfVertices(); ++ix ) { VertexBasePtr vertex = model->vertex(ix); if(vertex->getNpoint()!=3) continue; for(unsigned int iloc = 0;iloc < 3; ++iloc) { vector ext = vertex->search(iloc, Mediator_->id()); if(ext.empty()) continue; for(unsigned int ioff=0;ioffid() && ext[ioff+1]==incomingB_->id()) || (ext[ioff]==incomingB_->id() && ext[ioff+1]==incomingA_->id())) { foundDM = true; vertex->setCoupling(sqr(Emax),incomingA_,incomingB_,Mediator_); cDMmed_ = vertex->norm(); } else if(abs(ext[ioff])==1 && abs(ext[ioff+1])==1 && ext[ioff]==-ext[ioff+1]) { foundD = true; vertex->setCoupling(sqr(Emax),getParticleData(1),getParticleData(-1),Mediator_); cSMmed_[0] = vertex->norm(); } else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) { foundU = true; vertex->setCoupling(sqr(Emax),getParticleData(2),getParticleData(-2),Mediator_); cSMmed_[1] = vertex->norm(); } else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) { foundS = true; vertex->setCoupling(sqr(Emax),getParticleData(3),getParticleData(-3),Mediator_); cSMmed_[2] = vertex->norm(); } } } } if(!foundDM) { throw InitException() << "Cannot find DM coupling in MEDM2Mesons::doinit()"; } if(!foundD) { throw InitException() << "Cannot find down quark coupling in MEDM2Mesons::doinit()"; } if(!foundU) { throw InitException() << "Cannot find up quark coupling in MEDM2Mesons::doinit()"; } if(!foundS) { throw InitException() << "Cannot find strange quark coupling in MEDM2Mesons::doinit()"; } - cDMmed_ = 1.; - cSMmed_ = {1.0,1.0,1.0}; MEMultiChannel::doinit(); } void MEDM2Mesons::persistentOutput(PersistentOStream & os) const { os << current_ << modeMap_ << incomingA_ << incomingB_ << Mediator_ << cDMmed_ << cSMmed_; } void MEDM2Mesons::persistentInput(PersistentIStream & is, int) { is >> current_ >> modeMap_ >> incomingA_ >> incomingB_ >> Mediator_ >> cDMmed_ >> cSMmed_; } //The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigMEDM2Mesons("Herwig::MEDM2Mesons", "Herwig.so"); void MEDM2Mesons::Init() { static ClassDocumentation documentation ("The MEDM2Mesons class simulates the annhilation of" " DM particles to mesons at low energy"); static Reference interfaceWeakCurrent ("WeakCurrent", "The reference for the decay current to be used.", &MEDM2Mesons::current_, false, false, true, false, false); static Reference interfaceIncomingA ("IncomingA", "First incoming particle", &MEDM2Mesons::incomingA_, false, false, true, false, false); static Reference interfaceIncomingB ("IncomingB", "Second incoming particle", &MEDM2Mesons::incomingB_, false, false, true, false, false); static Reference interfaceMediator_ ("Mediator", "DM mediator", &MEDM2Mesons::Mediator_, 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 f1; vector 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 complex mmed = Mediator_->mass(); complex mmed2 = sqr(mmed); complex mwid = Mediator_->width(); complex prop = sHat()-mmed2+Complex(0.,1.)*mmed*mwid; complex pre = cDMmed_/prop; 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 constants(nOut+1); vector iSpin(nOut); vector 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 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 hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); vector hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); vector hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar), imode,ichan,q,out,momenta,DecayIntegrator::Calculate)); // compute the matrix element vector ihel(meMomenta().size()); double output(0.); unsigned int hI0_size = hadronI0.size(); unsigned int hI1_size = hadronI1.size(); unsigned int hss_size = hadronssbar.size(); unsigned int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size); for(unsigned int hhel=0;hhel0;--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]) { Complex amp = 0.; // work on coefficients for the I1 and I0 bits if(hI0_size != 0 ) amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel])); if(hI1_size !=0) amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI1[hhel])); if(hss_size !=0) amp += cSMmed_[2]* (lepton[ihel[0]][ihel[1]].dot(hadronssbar[hhel])); me_(ihel)= amp; output += std::norm(amp); } } } // symmetry factors map ncount; double symmetry(1.); for(tPDPtr o : out) ncount[o->id()]+=1; for(map::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/PDT/DMMediatorWidthGenerator.cc b/PDT/DMMediatorWidthGenerator.cc --- a/PDT/DMMediatorWidthGenerator.cc +++ b/PDT/DMMediatorWidthGenerator.cc @@ -1,217 +1,206 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the DMMediatorWidthGenerator class. // #include "DMMediatorWidthGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/RefVector.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 "Herwig/Decay/PhaseSpaceMode.h" #include "Herwig/Decay/General/DMMediatorDecayer.h" using namespace Herwig; namespace { struct ParticleOrdering { /** * Operator for the ordering * @param p1 The first ParticleData object * @param p2 The second ParticleData object */ bool operator() (tcPDPtr p1, tcPDPtr p2) const { return abs(p1->id()) > abs(p2->id()) || ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) || ( p1->id() == p2->id() && p1->fullName() > p2->fullName() ); } }; } -DMMediatorWidthGenerator::DMMediatorWidthGenerator() : cDMmed_(0.) { - cSMmed_ = {1.0,1.0,1.0}; -} +DMMediatorWidthGenerator::DMMediatorWidthGenerator() +{} IBPtr DMMediatorWidthGenerator::clone() const { return new_ptr(*this); } IBPtr DMMediatorWidthGenerator::fullclone() const { return new_ptr(*this); } void DMMediatorWidthGenerator::doinit() { for(tWeakCurrentPtr current : weakCurrents_) { current->init(); for(unsigned int imode=0;imodenumberOfModes();++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; // order the particles bool skip=false; for(unsigned int ix=0;ix outgoing(out.begin(),out.end()); string tag = parent_->PDGName() + "->"; bool first=false; for(tcPDPtr part : outgoing) { if(!first) first=true; else tag+=","; tag+=part->PDGName(); } tag+=";"; int charge(0); for(tcPDPtr part : outgoing) charge+=part->iCharge(); if(charge!=0) continue; // create the decayer ostringstream fullname; fullname << "/Herwig/Decays/DMMediator_" << parent_->PDGName(); for(tcPDPtr part : out) fullname << "_" << part->PDGName(); string classname = "Herwig::DMMediatorDecayer"; DMMediatorDecayerPtr decayer = dynamic_ptr_cast (generator()->preinitCreate(classname,fullname.str())); - decayer->setDecayInfo(parent_,out,current,cDMmed_,cSMmed_); + decayer->setDecayInfo(parent_,out,current); // // set decayer options from base class // setDecayerInterfaces(fullname.str()); // initialize the decayer decayer->init(); // calculate the width Energy pWidth = decayer->partialWidth(parent_,out); if(pWidth<=ZERO) { generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0"); continue; } tDMPtr ndm = generator()->preinitCreateDecayMode(tag); generator()->preinitInterface(ndm, "Decayer", "set", decayer->fullName()); parent_->stable(false); generator()->preinitInterface(ndm, "Active", "set", "Yes"); setBranchingRatio(ndm, pWidth); } } WidthGenerator::doinit(); cerr << "Width for " << parent_->PDGName() << " = " << parent_->width()/GeV << " GeV\n"; cerr << "Decay modes : \n"; for(auto mode : parent_->decayModes()) { cerr << mode->tag() << " " << mode->brat() << "\n"; } } void DMMediatorWidthGenerator::persistentOutput(PersistentOStream & os) const { - os << weakCurrents_ << cDMmed_ << cSMmed_; + os << weakCurrents_; } void DMMediatorWidthGenerator::persistentInput(PersistentIStream & is, int) { - is >> weakCurrents_ >> cDMmed_ >> cSMmed_; + is >> weakCurrents_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigDMMediatorWidthGenerator("Herwig::DMMediatorWidthGenerator", "Herwig.so"); void DMMediatorWidthGenerator::Init() { static ClassDocumentation documentation ("There is no documentation for the DMMediatorWidthGenerator class"); static RefVector interfaceWeakCurrents ("WeakCurrents", "Weak currents to use for the decays", &DMMediatorWidthGenerator::weakCurrents_, -1, false, false, true, false, false); static Reference interfaceParent ("Parent", "The decaying particle", &DMMediatorWidthGenerator::parent_, false, false, true, false, false); - - static Parameter interfacecDMmed - ("cDMmed", - "coupling of DM to dark mediator", - &DMMediatorWidthGenerator::cDMmed_, 1.0, 0., 10., false, false, Interface::limited); - - static ParVector interfacecSMmed - ("cSMmed", - "coupling of SM to dark mediator", - &DMMediatorWidthGenerator::cSMmed_, -1 , 1.0 , -10. , 10. , false, false, Interface::limited); } bool DMMediatorWidthGenerator::accept(const ParticleData & pd) const { if(pd.iSpin()!=PDT::Spin1) return false; if(parent_ && pd.id()!=parent_->id()) return false; return true; } Energy DMMediatorWidthGenerator::width(const ParticleData & pd, Energy) const { return pd.width(); } WidthGenerator::DecayMap DMMediatorWidthGenerator::rate(const ParticleData & pd) const { return pd.decaySelector(); } void DMMediatorWidthGenerator::setBranchingRatio(tDMPtr dm, Energy pwidth) { // if zero width just set BR to zero if(pwidth==ZERO) { generator()->preinitInterface(dm, "BranchingRatio","set", "0."); generator()->preinitInterface(dm, "OnOff","set", "Off"); return; } // Need width and branching ratios for all currently created decay modes DecaySet modes = parent_->decayModes(); unsigned int nmodes=0; for( auto dm : modes ) { if(dm->on()) ++nmodes; } if( nmodes==0 ) return; double dmbrat(0.); if( nmodes == 1 ) { parent_->width(pwidth); if( pwidth > ZERO ) parent_->cTau(hbarc/pwidth); dmbrat = 1.; parent_->width(pwidth); } else { Energy currentwidth(parent_->width()); Energy newWidth(currentwidth + pwidth); parent_->width(newWidth); if( newWidth > ZERO ) parent_->cTau(hbarc/newWidth); //need to reweight current branching fractions if there are any double factor = newWidth > ZERO ? double(currentwidth/newWidth) : 0.; for(DecaySet::const_iterator dit = modes.begin(); dit != modes.end(); ++dit) { if( **dit == *dm || !(**dit).on() ) continue; double newbrat = (**dit).brat()*factor; ostringstream brf; brf << setprecision(13)<< newbrat; generator()->preinitInterface(*dit, "BranchingRatio", "set", brf.str()); } //set brat for current mode dmbrat = newWidth > ZERO ? double(pwidth/newWidth) : 0.; parent_->width(newWidth); } ostringstream br; br << setprecision(13) << dmbrat; generator()->preinitInterface(dm, "BranchingRatio", "set", br.str()); } diff --git a/PDT/DMMediatorWidthGenerator.h b/PDT/DMMediatorWidthGenerator.h --- a/PDT/DMMediatorWidthGenerator.h +++ b/PDT/DMMediatorWidthGenerator.h @@ -1,158 +1,148 @@ // -*- C++ -*- #ifndef Herwig_DMMediatorWidthGenerator_H #define Herwig_DMMediatorWidthGenerator_H // // This is the declaration of the DMMediatorWidthGenerator class. // #include "ThePEG/PDT/WidthGenerator.h" #include "Herwig/Decay/WeakCurrents/WeakCurrent.h" namespace Herwig { using namespace ThePEG; /** * Here is the documentation of the DMMediatorWidthGenerator class. * * @see \ref DMMediatorWidthGeneratorInterfaces "The interfaces" * defined for DMMediatorWidthGenerator. */ class DMMediatorWidthGenerator: public WidthGenerator { public: /** * The default constructor. */ DMMediatorWidthGenerator(); /** @name Virtual functions to be overridden by sub-classes. */ //@{ /** * Return true if this object can be used for the given particle * type with the given decay map. */ virtual bool accept(const ParticleData &) const; /** * Given a particle type and a mass of an instance of that particle * type, calculate a width. */ virtual Energy width(const ParticleData &, Energy m) const; /** * Return decay map for the given particle type. */ virtual DecayMap rate(const ParticleData &) 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(); /** * Overloaded function from Interfaced */ virtual bool preInitialize() const { return true; } protected: /** * Set the branching ratio of this mode. This requires * calculating a new width for the decaying particle and reweighting * the current branching fractions. * @param dm The decaymode for which to set the branching ratio * @param pwidth The calculated width of the mode */ void setBranchingRatio(tDMPtr dm, Energy pwidth); 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: /** @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(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DMMediatorWidthGenerator & operator=(const DMMediatorWidthGenerator &); private: /** * The particle */ PDPtr parent_; /** * Weak currents to use for the decay */ vector weakCurrents_; - /** - * DM coupling to the dark mediator - */ - double cDMmed_; - - /** - * SM couplings to the dark mediator - */ - vector cSMmed_; - }; } #endif /* Herwig_DMMediatorWidthGenerator_H */ diff --git a/src/DM.in b/src/DM.in --- a/src/DM.in +++ b/src/DM.in @@ -1,82 +1,78 @@ # -*- ThePEG-repository -*- ################################################## # Read the model parameters ################################################## read DM.model ################################################## # Set the beams ################################################## read snippets/EECollider.in cd /Herwig/EventHandlers set EventHandler:BeamA /Herwig/Particles/chi set EventHandler:BeamB /Herwig/Particles/chi set /Herwig/Particles/chi:PDF /Herwig/Partons/NoPDF ################################################## # Set long-lived hadrons/muon to be unstable ################################################## set /Herwig/Particles/pi0:Stable Unstable set /Herwig/Particles/pi+:Stable Unstable set /Herwig/Particles/mu-:Stable Unstable set /Herwig/Particles/n0:Stable Unstable set /Herwig/Particles/K_L0:Stable Unstable set /Herwig/Particles/K+:Stable Unstable ################################################## # Selected the hard process ################################################## # leading-order processes ################################################## create Herwig::DMMediatorWidthGenerator /Herwig/Widths/DMWidth set /Herwig/Particles/Zp:Width_generator /Herwig/Widths/DMWidth set /Herwig/Widths/DMWidth:Parent /Herwig/Particles/Zp insert /Herwig/Widths/DMWidth:WeakCurrents 0 /Herwig/Decays/TwoKaonCzyzCurrent insert /Herwig/Widths/DMWidth:WeakCurrents 0 /Herwig/Decays/TwoPionCzyzCurrent insert /Herwig/Widths/DMWidth:WeakCurrents 0 /Herwig/Decays/PhiPiCurrent -set /Herwig/Widths/DMWidth:cDMmed 1.0 -set /Herwig/Widths/DMWidth:cSMmed[0] 1.0 -set /Herwig/Widths/DMWidth:cSMmed[1] 1.0 -set /Herwig/Widths/DMWidth:cSMmed[2] 1.0 cd /Herwig/MatrixElements # Set Processes insert SubProcess:MatrixElements 0 MEDM2Kaons ################################################## # LEP physics parameters (override defaults) ################################################## cd /Herwig/Generators # must be slightly greater than 2 x DM mass set EventGenerator:EventHandler:LuminosityFunction:Energy 2.0001 set /Herwig/Generators/EventGenerator:EventHandler:Cuts:MHatMin 0.2 cd /Herwig/Generators ################################################## ## prepare for Rivet analysis or HepMC output ## when running with parton shower ################################################## # create ThePEG::RivetAnalysis /Herwig/Analysis/Rivet RivetAnalysis.so # insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/Rivet # insert /Herwig/Analysis/Rivet:Analyses 0 ANALYSIS ################################################### # Save run for later usage with 'Herwig run' ################################################## set EventGenerator:MaxErrors 10000 set EventGenerator:EventHandler:StatLevel Full set EventGenerator:EventHandler:CascadeHandler NULL saverun DM EventGenerator