diff --git a/Decay/DecayIntegrator.cc b/Decay/DecayIntegrator.cc --- a/Decay/DecayIntegrator.cc +++ b/Decay/DecayIntegrator.cc @@ -1,292 +1,302 @@ // -*- C++ -*- // // DecayIntegrator.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 DecayIntegrator class. // #include "DecayIntegrator.h" #include "PhaseSpaceMode.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.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 "Herwig/PDT/WidthCalculatorBase.h" using namespace Herwig; void DecayIntegrator::persistentOutput(PersistentOStream & os) const { os << modes_ << nIter_ << nPoint_ << nTry_ - << photonGen_ << generateInter_ << ounit(eps_,GeV); + << photonGen_ << generateInter_ << ounit(eps_,GeV) << warnings_; } void DecayIntegrator::persistentInput(PersistentIStream & is, int) { is >> modes_ >> nIter_ >> nPoint_ >> nTry_ - >> photonGen_ >> generateInter_ >> iunit(eps_,GeV); + >> photonGen_ >> generateInter_ >> iunit(eps_,GeV) >> warnings_; } // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractClass describeHerwigDecayIntegrator("Herwig::DecayIntegrator", "Herwig.so"); void DecayIntegrator::Init() { static ClassDocumentation documentation ("The DecayIntegrator class is a base decayer class " "including a multi-channel integrator."); static Parameter interfaceIteration ("Iteration", "Number of iterations for the initialization of the phase space", &DecayIntegrator::nIter_, 10, 0, 100, false, false, true); static Parameter interfacePoints ("Points", "number of phase space points to generate in the initialisation.", &DecayIntegrator::nPoint_, 10000, 1, 1000000000, false, false, true); static Parameter interfaceNtry ("Ntry", "Number of attempts to generate the decay", &DecayIntegrator::nTry_, 500, 0, 100000, false, false, true); static Reference interfacePhotonGenerator ("PhotonGenerator", "Object responsible for generating photons in the decay.", &DecayIntegrator::photonGen_, false, false, true, true, false); static Switch interfaceGenerateIntermediates ("GenerateIntermediates", "Whether or not to include intermediate particles in the output", &DecayIntegrator::generateInter_, false, false, false); static SwitchOption interfaceGenerateIntermediatesNoIntermediates (interfaceGenerateIntermediates, "No", "Don't include the intermediates", false); static SwitchOption interfaceGenerateIntermediatesIncludeIntermediates (interfaceGenerateIntermediates, "Yes", "include the intermediates", true); + + static Switch InterfacePhaseSpaceWarning + ("PhaseSpaceWarning", + "Switch on/off text warnings in PhaseSpaceMode class", + &DecayIntegrator::warnings_, true, false, false); + static SwitchOption on + (InterfacePhaseSpaceWarning,"on","turn on the warnings", true); + static SwitchOption off + (InterfacePhaseSpaceWarning,"off","turn off the warnings", false); + } double DecayIntegrator::oneLoopVirtualME(unsigned int , const Particle &, const ParticleVector &) { throw Exception() << "DecayIntegrator::oneLoopVirtualME() called. This should" << " have been overidden in an inheriting class if it is used" << Exception::runerror; } InvEnergy2 DecayIntegrator::realEmissionME(unsigned int, const Particle &, ParticleVector &, unsigned int, double, double, const LorentzRotation &, const LorentzRotation &) { throw Exception() << "DecayIntegrator::realEmmisionME() called. This should" << " have been overidden in an inheriting class if it is used" << Exception::runerror; } ParticleVector DecayIntegrator::decay(const Particle & parent, const tPDVector & children) const { // return empty vector if products heavier than parent Energy mout(ZERO); for(tPDPtr pd : children) mout += pd->massMin(); if(mout>parent.mass()) return ParticleVector(); // generate the decay bool cc; iMode_ = modeNumber(cc,parent.dataPtr(),children); if(numberModes()==0) return ParticleVector(); return modes_[iMode_]->generateDecay(parent,this,generateInter_,cc); } void DecayIntegrator::doinitrun() { HwDecayerBase::doinitrun(); if ( initialize() && Debug::level > 1 ) CurrentGenerator::current().log() << "Start of the initialisation for " << name() << "\n"; for(unsigned int ix=0;ixinitrun(); iMode_=ix; modes_[ix]->initializePhaseSpace(initialize(),this); } } // output the information for the database void DecayIntegrator::dataBaseOutput(ofstream & output,bool header) const { // header for MySQL if(header) output << "update decayers set parameters=\""; output << "newdef " << name() << ":Iteration " << nIter_ << "\n"; output << "newdef " << name() << ":Ntry " << nTry_ << "\n"; output << "newdef " << name() << ":Points " << nPoint_ << "\n"; //if(_photongen){;} output << "newdef " << name() << ":GenerateIntermediates " << generateInter_ << " \n"; // footer for MySQL if(header) { output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";\n"; } } // set the code for the partial width void DecayIntegrator::setPartialWidth(const DecayMode & dm, int imode) { int ifound = findMode(dm); if(ifound>=0) modes_[ifound]->setPartialWidth(imode); } WidthCalculatorBasePtr DecayIntegrator::threeBodyMEIntegrator(const DecayMode &) const { return WidthCalculatorBasePtr(); } int DecayIntegrator::findMode(const DecayMode & dm) { int imode(-1); vector extid; bool found(false); int id; unsigned int ix(0),iy,N,iz,tmax,nmatched; if(modes_.size()==0) return -1; do { if(!modes_[ix]) { ++ix; continue; } tcPDPtr in = modes_[ix]->incoming().first; tcPDPtr cc = modes_[ix]->incoming().first->CC(); tmax=1;if(!cc){++tmax;} for(iz=0;izid()==in->id()&&iz==0) { for(iy=0,N=modes_[ix]->numberOfParticles();iyoutgoing()[iy]->id()); } } else if(dm.parent()->id()==in->id()&&iz==1) { for(iy=0,N=modes_[ix]->numberOfParticles();iyoutgoing()[iy]->CC(); extid.push_back( cc2 ? cc2->id() : modes_[ix]->outgoing()[iy]->id()); } } else if(cc&&dm.parent()->id()==cc->id()) { for(iy=0,N=modes_[ix]->numberOfParticles();iyoutgoing()[iy]->CC(); extid.push_back( cc2 ? cc2->id() : modes_[ix]->outgoing()[iy]->id()); } } // if the parents match if(!extid.empty()) { vector matched(extid.size(),false); bool done; nmatched=0; ParticleMSet::const_iterator pit = dm.products().begin(); do { id=(**pit).id(); done=false; iy=0; do { if(id==extid[iy]&&!matched[iy]) { matched[iy]=true; ++nmatched; done=true; } ++iy; } while(iygenerateDecay(inpart,this,inter,cc); } void DecayIntegrator::addMode(PhaseSpaceModePtr mode) const { modes_.push_back(mode); if(mode) mode->init(); } ostream & Herwig::operator<<(ostream & os, const DecayIntegrator & decay) { os << "The integrator has " << decay.modes_.size() << " modes" << endl; for(unsigned int ix=0;ixresetIntermediate(part,mass,width); } } Energy DecayIntegrator::initializePhaseSpaceMode(unsigned int imode,bool init, bool onShell) const{ tcPhaseSpaceModePtr cmodeptr=mode(imode); tPhaseSpaceModePtr modeptr = const_ptr_cast(cmodeptr); modeptr->init(); return modeptr->initializePhaseSpace(init,this,onShell); } diff --git a/Decay/DecayIntegrator.h b/Decay/DecayIntegrator.h --- a/Decay/DecayIntegrator.h +++ b/Decay/DecayIntegrator.h @@ -1,515 +1,526 @@ // -*- C++ -*- // // DecayIntegrator.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_DecayIntegrator_H #define Herwig_DecayIntegrator_H // // This is the declaration of the DecayIntegrator class. // #include "DecayIntegrator.fh" #include "HwDecayerBase.h" #include "PhaseSpaceMode.fh" #include "Herwig/PDT/WidthCalculatorBase.fh" #include "Radiation/DecayRadiationGenerator.h" #include namespace Herwig { using namespace ThePEG; /** \ingroup Decay * \class DecayIntegrator * \brief Main class for Decayers implementing multi-channel phase space integration. * \author Peter Richardson * * This class is designed to be the base class for Herwig decays including * the implementation of a multichannel decayer or n-body phase space decays. * * The DecayIntegrator class inherits from ThePEG's Decayer class * and makes use of the PhaseSpaceMode class to specify a number * of decay modes. * * Additional modes can be added using the addMode method. In practice the * phase space channels for a particular mode are usually constructed in the * doinit member of a Decayer and then the modes added to the Decayer. * * For the majority of the decays currently implemented the * phase-space integration has been optimised and the maximum weight set. * If the parameters of the decay model are changed the Initialize interface * can be used to optimise the integration and calculate the maximum weight. * * In classes inheriting from this the me2() member which gives the matrix element * squared must be implemented. This should be combined with the setting of the * phase space channels, and the setting of which channels to use and their * initial weights in the doinit() member. The different decay modes should then * be initialized in the initrun() member if needed. The generate member can then * be called from the decay() member to generate a phase-space configuration for a * decay. * * @see DecayPhaseSpaceMode * @see DecayPhaseSpaceChannel * @see \ref DecayIntegratorInterfaces "The interfaces" * defined for DecayIntegrator. */ class DecayIntegrator: public HwDecayerBase { public: /** * and DecayPhaseMode */ friend class PhaseSpaceMode; /** * Enum for the matrix element option */ enum MEOption {Initialize,Calculate,Terminate}; public: /** * The default constructor. */ DecayIntegrator() : nIter_(10), nPoint_(10000), nTry_(500), generateInter_(false), iMode_(-1), - realME_(false), virtualME_(false), eps_(ZERO) + realME_(false), virtualME_(false), eps_(ZERO), warnings_(true) {} public: /** * Check if this decayer can perfom the decay for a particular mode. * Uses the modeNumber member but can be overridden * @param parent The decaying particle * @param children The decay products */ virtual bool accept(tcPDPtr parent, const tPDVector & children) const { bool cc; return modeNumber(cc,parent,children)>=0; } /** * For a given decay mode and a given particle instance, perform the * decay and return the decay products. As this is the base class this * is not implemented. * @return The vector of particles produced in the decay. */ virtual ParticleVector decay(const Particle & parent, const tPDVector & children) const; /** * 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 = 0; /** * The mode being used for this decay */ int imode() const {return iMode_;} /** * Add a phase-space mode to the list * @param mode The mode being added. */ void addMode(PhaseSpaceModePtr mode) 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. */ virtual double me2(const int ichan, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const = 0; /** * Construct the SpinInfos for the particles produced in the decay */ virtual void constructSpinInfo(const Particle & part, ParticleVector outgoing) const = 0; /** * 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; /** * Set the code for the partial width. Finds the partial width in the * GenericWidthGenerator class which corresponds to the decay mode. * @param dm The DecayMode * @param imode The mode. */ void setPartialWidth(const DecayMode & dm, int imode); /** * Specify the \f$1\to2\f$ matrix element to be used in the running width calculation. * @param mecode The code for the matrix element as described * in the GenericWidthGenerator class. * @param coupling The coupling for the matrix element. * @return True or False if this mode can be handled. */ virtual bool twoBodyMEcode(const DecayMode &, int & mecode, double & coupling) const { coupling = 1.; mecode = -1; return false; } /** * Method to return an object to calculate the 3 (or higher body) partial width * @param dm The DecayMode * @return A pointer to a WidthCalculatorBase object capable of calculating the width */ virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const; /** * The matrix element to be integrated for the three-body decays as a function * of the invariant masses of pairs of the outgoing particles. * @param imode The mode for which the matrix element is needed. * @param q2 The scale, \e i.e. the mass squared of the decaying particle. * @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$. * @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$. * @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$. * @param m1 The mass of the first outgoing particle. * @param m2 The mass of the second outgoing particle. * @param m3 The mass of the third outgoing particle. * @return The matrix element */ virtual double threeBodyMatrixElement(const int imode, const Energy2 q2, const Energy2 s3, const Energy2 s2, const Energy2 s1, const Energy m1, const Energy m2, const Energy m3) const; /** * The differential three body decay rate with one integral performed. * @param imode The mode for which the matrix element is needed. * @param q2 The scale, \e i.e. the mass squared of the decaying particle. * @param s The invariant mass which still needs to be integrate over. * @param m1 The mass of the first outgoing particle. * @param m2 The mass of the second outgoing particle. * @param m3 The mass of the third outgoing particle. * @return The differential rate \f$\frac{d\Gamma}{ds}\f$ */ virtual InvEnergy threeBodydGammads(const int imode, const Energy2 q2, const Energy2 s, const Energy m1, const Energy m2, const Energy m3) const; /** * Finds the phase-space mode corresponding to a given decay mode * @param dm The DecayMode */ int findMode(const DecayMode & dm); public: /** * Members for the generation of QED radiation in the decays */ //@{ /** * Use the DecayRadiationGenerator to generate photons in the decay. * @param p The Particle instance being decayed * @param children The decay products * @return A particle vector containing the decay products after the generation * of photons. */ ParticleVector generatePhotons(const Particle & p,ParticleVector children) { return photonGen_->generatePhotons(p,children,this); } /** * check if photons can be generated in the decay */ bool canGeneratePhotons() {return photonGen_;} /** * The one-loop virtual correction. * @param imode The mode required. * @param part The decaying particle. * @param products The decay products including the radiated photon. * @return Whether the correction is implemented */ virtual double oneLoopVirtualME(unsigned int imode, const Particle & part, const ParticleVector & products); /** * Whether or not the one loop matrix element is implemented */ bool hasOneLoopME() {return virtualME_;} /** * The real emission matrix element * @param imode The mode required * @param part The decaying particle * @param products The decay products including the radiated photon * @param iemitter The particle which emitted the photon * @param ctheta The cosine of the polar angle between the photon and the * emitter * @param stheta The sine of the polar angle between the photon and the * emitter * @param rot1 Rotation from rest frame to frame for real emission * @param rot2 Rotation to place emitting particle along z */ virtual InvEnergy2 realEmissionME(unsigned int imode, const Particle & part, ParticleVector & products, unsigned int iemitter, double ctheta, double stheta, const LorentzRotation & rot1, const LorentzRotation & rot2); /** * Whether or not the real emission matrix element is implemented */ bool hasRealEmissionME() {return realME_;} //@} public: /** * The output operator is a friend, this is mainly for debugging */ friend ostream & operator<<(ostream & os, const DecayIntegrator & decay); 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. Called in the run phase just before * a run begins. */ virtual void doinitrun(); //@} protected: /** * Generate the momenta for the decay * @param inter Generate the intermediates produced in the decay as well as the * final particles. * @param cc Is this the mode defined or its charge conjugate. * @param imode The mode being generated. * @param inpart The decaying particle. * @return The particles produced inthe decay. */ ParticleVector generate(bool inter,bool cc, const unsigned int & imode, const Particle & inpart) const; /** * Set the mode being use for this decay. */ void imode(int in) { iMode_ = in;} /** * Set the helicity matrix element for the decay. */ void ME(DecayMEPtr in) const { matrixElement_ = in;} /** * The helicity amplitude matrix element for spin correlations. */ DecayMEPtr ME() const {return matrixElement_;} /** * Reset the properities of all intermediates. * @param part The intermediate particle being reset. * @param mass The mass of the particle. * @param width The width of the particle. */ void resetIntermediate(tcPDPtr part, Energy mass, Energy width); /** * Initialize the phase-space mode * @param imode The mode * @param init Whether or not to perform the initialization */ Energy initializePhaseSpaceMode(unsigned int imode,bool init, bool onShell=false) const; protected: /** * Methods to set variables in inheriting classes */ //@{ /** * Set whether or not the intermediates are included */ void generateIntermediates(bool in) {generateInter_=in;} /** * Set whether or not the intermediates are included */ bool generateIntermediates() const {return generateInter_;} /** * Whether or not the one loop matrix element is implemented */ void hasOneLoopME(bool in) {virtualME_=in;} /** * Whether or not the real emission matrix element is implemented */ void hasRealEmissionME(bool in) {realME_=in;} /** * Set the epsilon parameter */ void epsilonPS(Energy in) {eps_=in;} /** * Clear the models */ void clearModes() {modes_.clear();} protected: /** * Number of decay modes */ unsigned int numberModes() const {return modes_.size();} /** * Pointer to a mode */ tPhaseSpaceModePtr mode(unsigned int ix) { return modes_[ix]; } /** * Pointer to a mode */ tcPhaseSpaceModePtr mode(unsigned int ix) const { return modes_[ix]; } +public: + +bool warnings() const { + return warnings_; +} + private: /** * Private and non-existent assignment operator. */ DecayIntegrator & operator=(const DecayIntegrator &) = delete; /** * Parameters for the integration */ //@{ /** * Number of iterations for th initialization. */ unsigned int nIter_; /** * Number of points for initialisation */ unsigned int nPoint_; /** * number of attempts to generate the decay */ unsigned int nTry_; /** * List of the decay modes */ mutable vector modes_; //@} /** * Whether to include the intermediates whne outputing the results. */ bool generateInter_; /** * Pointer to the object generating the QED radiation in the decay */ DecayRadiationGeneratorPtr photonGen_; /** * mode currently being generated */ mutable int iMode_; /** * The helicity matrix element for the current decay */ mutable DecayMEPtr matrixElement_; /** * Whether or not the real photon emission matrix element exists */ bool realME_; /** * Whether or not the one-loop matrix element exists */ bool virtualME_; /** * Epsilon parameter for phase-space integration */ Energy eps_; + /** + * option for turinh on/off log warnings in Phase class + */ + bool warnings_; + protected: /** * Exception for this class and those inheriting from it */ class DecayIntegratorError: public Exception {}; }; /** * Output information on the DecayIntegrator for debugging purposes */ ostream & operator<<(ostream &, const DecayIntegrator &); } #endif /* Herwig_DecayIntegrator_H */ diff --git a/Decay/PhaseSpaceMode.cc b/Decay/PhaseSpaceMode.cc --- a/Decay/PhaseSpaceMode.cc +++ b/Decay/PhaseSpaceMode.cc @@ -1,522 +1,527 @@ // -*- C++ -*- // // PhaseSpaceMode.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 PhaseSpaceMode class. // #include "PhaseSpaceMode.h" #include "ThePEG/Interface/ClassDocumentation.h" +#include "ThePEG/Interface/Switch.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/Utilities/Kinematics.h" using namespace Herwig; using namespace ThePEG::Helicity; void PhaseSpaceMode::persistentOutput(PersistentOStream & os) const { os << incoming_ << channels_ << maxWeight_ << outgoing_ << outgoingCC_ << partial_ << widthGen_ << massGen_ << BRsum_ << testOnShell_ << ounit(eMax_,GeV) << nRand_; } void PhaseSpaceMode::persistentInput(PersistentIStream & is, int) { is >> incoming_ >> channels_ >> maxWeight_ >> outgoing_ >> outgoingCC_ >> partial_ >> widthGen_ >> massGen_ >> BRsum_ >> testOnShell_ >> iunit(eMax_,GeV) >> nRand_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigPhaseSpaceMode("Herwig::PhaseSpaceMode", "libHerwig.so"); void PhaseSpaceMode::Init() { static ClassDocumentation documentation ("The PhaseSpaceMode classs contains a number of phase space" " channels for the integration of a particular decay mode"); } ParticleVector PhaseSpaceMode::generateDecay(const Particle & inpart, tcDecayIntegratorPtr decayer, bool intermediates,bool cc) { decayer->ME(DecayMEPtr()); eps_=decayer->eps_; // compute the prefactor Energy prewid = (widthGen_&&partial_>=0) ? widthGen_->partialWidth(partial_,inpart.mass()) : incoming_.first->width(); InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV; // boosts to/from rest Boost bv =-inpart.momentum().boostVector(); double gammarest = inpart.momentum().e()/inpart.momentum().mass(); LorentzRotation boostToRest( bv,gammarest); LorentzRotation boostFromRest(-bv,gammarest); // construct a new particle which is at rest Particle inrest(inpart); inrest.transform(boostToRest); double wgt(0.); vector momenta(outgoing_.size()); int ichan; unsigned int ncount(0); try { do { // phase-space pieces of the weight fillStack(); wgt = pre*weight(ichan,inrest.momentum(),momenta); // matrix element piece wgt *= decayer->me2(-1,inrest,!cc ? outgoing_ : outgoingCC_, momenta, ncount ==0 ? DecayIntegrator::Initialize : DecayIntegrator::Calculate); ++ncount; if(wgt>maxWeight_) { + if(decayer->warnings()) { CurrentGenerator::log() << "Resetting max weight for decay " << inrest.PDGName() << " -> "; for(tcPDPtr part : outgoing_) CurrentGenerator::log() << " " << part->PDGName(); CurrentGenerator::log() << " " << maxWeight_ << " " << wgt << " " << inrest.mass()/MeV << "\n"; + } maxWeight_=wgt; } } while(maxWeight_*UseRandom::rnd()>wgt&&ncountnTry_); } catch (Veto) { // restore the incoming particle to its original state inrest.transform(boostFromRest); while(!rStack_.empty()) rStack_.pop(); throw Veto(); } if(ncount>=decayer->nTry_) { + if(decayer->warnings()) { CurrentGenerator::log() << "The decay " << inrest.PDGName() << " -> "; for(tcPDPtr part : outgoing_) CurrentGenerator::log() << " " << part->PDGName(); CurrentGenerator::log() << " " << maxWeight_ << " " << decayer->nTry_ << " is too inefficient for the particle " << inpart << "vetoing the decay \n"; + } momenta.clear(); throw Veto(); } // construct the particles ParticleVector output(outgoing_.size()); for(unsigned int ix=0;ixproduceParticle(momenta[ix]); // set up the spin correlations decayer->constructSpinInfo(inrest,output); const_ptr_cast(&inpart)->spinInfo(inrest.spinInfo()); constructVertex(inpart,output,decayer); // return if intermediate particles not required if(channels_.empty()||!intermediates) { for(tPPtr part : output) part->transform(boostFromRest); } // find the intermediate particles else { // select the channel vector mewgts(channels_.size(),0.0); double total=0.; for(unsigned int ix=0,N=channels_.size();ixme2(ix,inrest,!cc ? outgoing_ : outgoingCC_, momenta,DecayIntegrator::Calculate); total+=mewgts[ix]; } // randomly pick a channel total *= UseRandom::rnd(); int iChannel = -1; do { ++iChannel; total-=mewgts[iChannel]; } while(iChannel0.); iChannel_ = iChannel; // apply boost for(tPPtr part : output) part->transform(boostFromRest); channels_[iChannel].generateIntermediates(cc,inpart,output); } decayer->ME(DecayMEPtr()); // return particles; return output; } // flat phase space generation and weight Energy PhaseSpaceMode::flatPhaseSpace(const Lorentz5Momentum & in, vector & momenta, bool onShell) const { double wgt(1.); // masses of the particles vector mass = externalMasses(in.mass(),wgt,onShell); // two body decay double ctheta = 2.*rStack_.top() - 1.; rStack_.pop(); double phi = Constants::twopi*rStack_.top(); rStack_.pop(); if(! Kinematics::twoBodyDecay(in, mass[0], mass[1], ctheta, phi, momenta[0],momenta[1])) throw Exception() << "Incoming mass - Outgoing mass negative in " << "PhaseSpaceMode::flatPhaseSpace()" << Exception::eventerror; wgt *= Kinematics::pstarTwoBodyDecay(in.mass(),mass[0],mass[1]) /8./Constants::pi/in.mass(); return wgt*in.mass(); } // generate the masses of the external particles vector PhaseSpaceMode::externalMasses(Energy inmass,double & wgt, bool onShell) const { vector mass(outgoing_.size()); vector notdone; Energy mlow(ZERO); // set masses of stable particles and limits for(unsigned int ix=0;ixmass(); if(massGen_[ix]) rStack_.pop(); } else if(!massGen_[ix] || outgoing_[ix]->stable()) { if(massGen_[ix]) rStack_.pop(); mass[ix] = outgoing_[ix]->generateMass(); mlow += mass[ix]; } else { mass[ix] = ZERO; notdone.push_back(ix); mlow+=max(outgoing_[ix]->mass()-outgoing_[ix]->widthLoCut(),ZERO); } } if(mlow>inmass) throw Veto(); // now we need to generate the masses for the particles we haven't for( ;!notdone.empty();) { unsigned int iloc=long(UseRandom::rnd()*(notdone.size()-1)); Energy low = max(outgoing_[notdone[iloc]]->mass()-outgoing_[notdone[iloc]]->widthLoCut(),ZERO); mlow-=low; double wgttemp; mass[notdone[iloc]]= massGen_[notdone[iloc]]->mass(wgttemp,*outgoing_[notdone[iloc]],low,inmass-mlow,rStack_.top()); wgttemp /= BRsum_[notdone[iloc]]; rStack_.pop(); assert(mass[notdone[iloc]]>=low&&mass[notdone[iloc]]<=inmass-mlow); wgt *= wgttemp; mlow += mass[notdone[iloc]]; notdone.erase(notdone.begin()+iloc); } return mass; } // construct the vertex for spin corrections void PhaseSpaceMode::constructVertex(const Particle & inpart, const ParticleVector & decay, tcDecayIntegratorPtr decayer) const { // construct the decay vertex VertexPtr vertex(new_ptr(DecayVertex())); DVertexPtr Dvertex(dynamic_ptr_cast(vertex)); // set the incoming particle for the decay vertex (inpart.spinInfo())->decayVertex(vertex); for(unsigned int ix=0;ixspinInfo()->productionVertex(vertex); } // set the matrix element Dvertex->ME(decayer->ME()); decayer->ME(DecayMEPtr()); } // initialise the phase space Energy PhaseSpaceMode::initializePhaseSpace(bool init, tcDecayIntegratorPtr decayer, bool onShell) { decayer->ME(DecayMEPtr()); eps_=decayer->eps_; Energy output(ZERO); // ensure that the weights add up to one if(!channels_.empty()) { double temp=0.; for(const PhaseSpaceChannel & channel : channels_) temp+= channel.weight(); for(PhaseSpaceChannel & channel : channels_) channel.weight(channel.weight()/temp); } if(!init) return ZERO; ThePEG::PPtr inpart=incoming_.first->produceParticle(); // now if using flat phase space maxWeight_=0.; if(channels_.empty()) { vector momenta(outgoing_.size()); double wsum=0.,wsqsum=0.; Energy m0,mmin(ZERO); for(tcPDPtr part : outgoing_) mmin += part->massMin(); for(unsigned int ix=0;ixnPoint_;++ix) { // set the mass of the decaying particle m0 = !onShell ? inpart->dataPtr()->generateMass() : inpart->dataPtr()->mass(); double wgt=0.; if(m0<=mmin) continue; inpart->set5Momentum(Lorentz5Momentum(m0)); // compute the prefactor Energy prewid = (widthGen_&&partial_>=0) ? widthGen_->partialWidth(partial_,inpart->mass()) : incoming_.first->width(); InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV; // generate the weight for this point try { int dummy; // phase-space piece fillStack(); wgt = pre*weight(dummy,inpart->momentum(),momenta,onShell); // matrix element piece wgt *= decayer->me2(-1,*inpart,outgoing_,momenta,DecayIntegrator::Initialize); } catch (Veto) { while(!rStack_.empty()) rStack_.pop(); wgt=0.; } if(wgt>maxWeight_) maxWeight_ = wgt; wsum += wgt; wsqsum += sqr(wgt); } wsum /= decayer->nPoint_; wsqsum=wsqsum/decayer->nPoint_-sqr(wsum); if(wsqsum<0.) wsqsum=0.; wsqsum=sqrt(wsqsum/decayer->nPoint_); Energy fact = (widthGen_&&partial_>=0) ? widthGen_->partialWidth(partial_,inpart->nominalMass()) : inpart->dataPtr()->width(); if(fact==ZERO) fact=MeV; // factor for the weight with spin correlations maxWeight_ *= inpart->dataPtr()->iSpin()==1 ? 1.1 : 1.6; if ( Debug::level > 1 ) { // ouptut the information on the initialisation CurrentGenerator::log() << "Initialized the phase space for the decay " << incoming_.first->PDGName() << " -> "; for(tPDPtr part : outgoing_) CurrentGenerator::log() << part->PDGName() << " "; CurrentGenerator::log() << "\n"; if(fact!=MeV) CurrentGenerator::log() << "The branching ratio is " << wsum << " +/- " << wsqsum << "\n"; CurrentGenerator::log() << "The partial width is " << wsum*fact/MeV << " +/- " << wsqsum*fact/MeV << " MeV\n"; CurrentGenerator::log() << "The partial width is " << wsum*fact/6.58212E-22/MeV << " +/- " << wsqsum*fact/6.58212E-22/MeV<< " s-1\n"; CurrentGenerator::log() << "The maximum weight is " << maxWeight_ << endl; } output=wsum*fact; } else { vector momenta(outgoing_.size()); double totsum(0.),totsq(0.); for(unsigned int iy=0;iynIter_;++iy) { // zero the maximum weight maxWeight_=0.; vector wsum(channels_.size(),0.),wsqsum(channels_.size(),0.); vector nchan(channels_.size(),0); totsum = 0.; totsq = 0.; Energy mmin(ZERO); for(tcPDPtr part : outgoing_) mmin += part->massMin(); for(unsigned int ix=0;ixnPoint_;++ix) { Energy m0 = !onShell ? incoming_.first->generateMass() : incoming_.first->mass(); double wgt=0.; int ichan(-1); if(m0>mmin) { inpart->set5Momentum(Lorentz5Momentum(m0)); // compute the prefactor Energy prewid= (widthGen_&&partial_>=0) ? widthGen_->partialWidth(partial_,inpart->mass()) : inpart->dataPtr()->width(); InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV; // generate the weight for this point try { fillStack(); wgt = pre*weight(ichan,inpart->momentum(),momenta,onShell); // matrix element piece wgt *= decayer->me2(-1,*inpart,outgoing_,momenta,DecayIntegrator::Initialize); } catch (Veto) { wgt=0.; } } if(wgt>maxWeight_) maxWeight_=wgt; if(ichan>=0) { wsum[ichan] += wgt; wsqsum[ichan] += sqr(wgt); ++nchan[ichan]; } totsum+=wgt; totsq+=wgt*wgt; } totsum=totsum/decayer->nPoint_; totsq=totsq/decayer->nPoint_-sqr(totsum); if(totsq<0.) totsq=0.; totsq=sqrt(totsq/decayer->nPoint_); if ( Debug::level > 1 ) CurrentGenerator::log() << "The branching ratio is " << iy << " " << totsum << " +/- " << totsq << maxWeight_ << "\n"; // compute the individual terms double total(0.); for(unsigned int ix=0;ix0.) { for(unsigned int ix=0;ixdataPtr()->iSpin()==1 ? 1.1 : 1.6; // output the information on the initialisation Energy fact = (widthGen_&&partial_>=0) ? widthGen_->partialWidth(partial_,inpart->nominalMass()) : inpart->dataPtr()->width(); if(fact==ZERO) fact=MeV; output=totsum*fact; if ( Debug::level > 1 ) { CurrentGenerator::log() << "Initialized the phase space for the decay " << incoming_.first->PDGName() << " -> "; for(tcPDPtr part : outgoing_) CurrentGenerator::log() << part->PDGName() << " "; CurrentGenerator::log() << "\n"; if(fact!=MeV) CurrentGenerator::log() << "The branching ratio is " << totsum << " +/- " << totsq << "\n"; CurrentGenerator::log() << "The partial width is " << totsum*fact/MeV << " +/- " << totsq*fact/MeV << " MeV\n"; CurrentGenerator::log() << "The partial width is " << totsum*fact/6.58212E-22/MeV << " +/- " << totsq*fact/6.58212E-22/MeV << " s-1\n"; CurrentGenerator::log() << "The maximum weight is " << maxWeight_ << "\n"; CurrentGenerator::log() << "The weights for the different phase" << " space channels are \n"; for(unsigned int ix=0,N=channels_.size();ixME(DecayMEPtr()); return output; } // generate a phase-space point using multichannel phase space Energy PhaseSpaceMode::channelPhaseSpace(int & ichan, const Lorentz5Momentum & in, vector & momenta, bool onShell) const { double wgt(rStack_.top()); rStack_.pop(); // select a channel ichan=-1; do { ++ichan; wgt-=channels_[ichan].weight(); } while(ichan0.); // generate the momenta if(ichan==int(channels_.size())) { throw Exception() << "PhaseSpaceMode::channelPhaseSpace()" << " failed to select a channel" << Exception::abortnow; } // generate the masses of the external particles double masswgt(1.); vector mass(externalMasses(in.mass(),masswgt,onShell)); momenta=channels_[ichan].generateMomenta(in,mass); // compute the denominator of the weight wgt=0.; for(const PhaseSpaceChannel & channel : channels_) wgt += channel.generateWeight(momenta); // return the weight return wgt!=0. ? in.mass()*masswgt/wgt : ZERO; } void PhaseSpaceMode::init() { // get mass and width generators outgoingCC_.clear(); for(tPDPtr part : outgoing_) { outgoingCC_.push_back(part->CC() ? part->CC() : part); } massGen_.resize(outgoing_.size()); BRsum_.resize(outgoing_.size()); widthGen_ = dynamic_ptr_cast(incoming_.first->widthGenerator()); for(unsigned int ix=0;ix(outgoing_[ix]->massGenerator()); if(outgoing_[ix]->stable()) { BRsum_[ix] = 1.; } else { double total(0.); for(tDMPtr mode : outgoing_[ix]->decayModes()) { if(mode->on()) total += mode->brat(); } BRsum_[ix] = total; } } // get max energy for decays if(!incoming_.second) eMax_ = testOnShell_ ? incoming_.first->mass() : incoming_.first->massMax(); // work out how many random numbers we need nRand_ = 3*outgoing_.size()-4; for(unsigned int ix=0;ixwidthGenerator()!=widthGen_) widthGen_ = dynamic_ptr_cast(incoming_.first->widthGenerator()); for(unsigned int ix=0;ixmassGenerator()) massGen_[ix] = dynamic_ptr_cast(outgoing_[ix]->massGenerator()); } for(PhaseSpaceChannel & channel : channels_) channel.initrun(this); if(widthGen_) const_ptr_cast(widthGen_)->initrun(); tcGenericWidthGeneratorPtr wtemp; for(unsigned int ix=0;ix(outgoing_[ix]->widthGenerator()); if(wtemp) const_ptr_cast(wtemp)->initrun(); } // ensure weights sum to one double sum(0.); for(const PhaseSpaceChannel & channel : channels_) sum+=channel.weight(); for(PhaseSpaceChannel & channel : channels_) channel.weight(channel.weight()/sum); nRand_ = 3*outgoing_.size()-4; for(unsigned int ix=0;ix