diff --git a/Decay/DecayIntegrator.cc b/Decay/DecayIntegrator.cc --- a/Decay/DecayIntegrator.cc +++ b/Decay/DecayIntegrator.cc @@ -1,386 +1,381 @@ // -*- C++ -*- // // DecayIntegrator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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. // // Author: Peter Richardson // #include "DecayIntegrator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Utilities/Debug.h" #include "Herwig/Utilities/Kinematics.h" #include "Herwig/PDT/GenericMassGenerator.h" #include "DecayPhaseSpaceMode.h" #include "Herwig/PDT/WidthCalculatorBase.h" #include "ThePEG/Interface/Reference.h" using namespace Herwig; ParticleVector DecayIntegrator::decay(const Particle & parent, const tPDVector & children) const { // return empty vector if products heavier than parent Energy mout(ZERO); for(tPDVector::const_iterator it=children.begin(); it!=children.end();++it) mout+=(**it).massMin(); if(mout>parent.mass()) return ParticleVector(); // generate the decay bool cc; _imode = modeNumber(cc,parent.dataPtr(),children); return _modes[_imode]->generate(_generateinter,cc,parent); } void DecayIntegrator::persistentOutput(PersistentOStream & os) const { os << _modes << _niter << _npoint << _ntry << _photongen << _generateinter << ounit(_eps,GeV); } void DecayIntegrator::persistentInput(PersistentIStream & is, int) { is >> _modes >> _niter >> _npoint >> _ntry >> _photongen >> _generateinter >> iunit(_eps,GeV); } // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractClass describeHerwigDecayIntegrator("Herwig::DecayIntegrator", "Herwig.so"); void DecayIntegrator::Init() { - - static RefVector interfaceModes - ("Modes", - "The phase space integration modes.", - &DecayIntegrator::_modes, 0, false, false, true, true); 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); } // output info on the integrator ostream & Herwig::operator<<(ostream & os, const DecayIntegrator & decay) { os << "The integrator has " << decay._modes.size() << " modes" << endl; for(unsigned int ix=0;ixgenerate(inter,cc,inpart); } // initialization for a run void DecayIntegrator::doinitrun() { HwDecayerBase::doinitrun(); if ( initialize() && Debug::level > 1 ) CurrentGenerator::current().log() << "Start of the initialisation for " << this->name() << "\n"; for(unsigned int ix=0;ix<_modes.size();++ix) { if(!_modes[ix]) continue; _modes[ix]->initrun(); _imode=ix; _modes[ix]->initializePhaseSpace(initialize()); } // CurrentGenerator::current().log() << *this << "\n"; } // add a new mode void DecayIntegrator::addMode(DecayPhaseSpaceModePtr in,double maxwgt, const vector inwgt) const { _modes.push_back(in); if(in) { in->setMaxWeight(maxwgt); in->setWeights(inwgt); in->setIntegrate(_niter,_npoint,_ntry); in->init(); } } // reset the properities of all intermediates void DecayIntegrator::resetIntermediate(tcPDPtr part, Energy mass, Energy width) { if(!part) return; for(unsigned int ix=0,N=_modes.size();ixresetIntermediate(part,mass,width); } } WidthCalculatorBasePtr DecayIntegrator::threeBodyMEIntegrator(const DecayMode &) const { return WidthCalculatorBasePtr(); } // set the code for the partial width void DecayIntegrator::setPartialWidth(const DecayMode & dm, int imode) { vector extid; tcPDPtr cc,cc2; int nfound(0),ifound,nmax(1),id; unsigned int ix(0),iy,N,iz,tmax,nmatched; if(dm.parent()->CC()) nmax=2; if(_modes.size()==0) return; do { if(!_modes[ix]) { ++ix; continue; } cc = _modes[ix]->externalParticles(0)->CC(); tmax = cc ? 1 : 2; for(iz=0;izid()==_modes[ix]->externalParticles(0)->id()&&iz==0) { for(iy=0,N=_modes[ix]->numberofParticles();iyexternalParticles(iy)->id()); } } else if(dm.parent()->id()==_modes[ix]->externalParticles(0)->id()&&iz==1) { for(iy=0,N=_modes[ix]->numberofParticles();iyexternalParticles(iy)->CC(); extid.push_back( cc2 ? cc2->id() : _modes[ix]->externalParticles(iy)->id()); } } else if(cc&&dm.parent()->id()==cc->id()) { for(iy=0,N=_modes[ix]->numberofParticles();iyexternalParticles(iy)->CC(); extid.push_back( cc2 ? cc2->id() : _modes[ix]->externalParticles(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=1; do { if(id==extid[iy]&&!matched[iy]) { matched[iy]=true; ++nmatched; done=true; } ++iy; } while(iy=0) _modes[ifound]->setPartialWidth(imode); } ++ix; } while(nfound extid; tcPDPtr cc,cc2; 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; } cc = _modes[ix]->externalParticles(0)->CC(); tmax=1;if(!cc){++tmax;} for(iz=0;izid()==_modes[ix]->externalParticles(0)->id()&&iz==0) { for(iy=0,N=_modes[ix]->numberofParticles();iyexternalParticles(iy)->id()); } } else if(dm.parent()->id()==_modes[ix]->externalParticles(0)->id()&&iz==1) { for(iy=0,N=_modes[ix]->numberofParticles();iyexternalParticles(iy)->CC(); extid.push_back( cc2 ? cc2->id() : _modes[ix]->externalParticles(iy)->id()); } } else if(cc&&dm.parent()->id()==cc->id()) { for(iy=0,N=_modes[ix]->numberofParticles();iyexternalParticles(iy)->CC(); extid.push_back( cc ? cc->id() : _modes[ix]->externalParticles(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=1; do { if(id==extid[iy]&&!matched[iy]) { matched[iy]=true; ++nmatched; done=true; } ++iy; } while(iy(cmodeptr); modeptr->init(); return modeptr->initializePhaseSpace(init,onShell); } // the matrix element to be integrated for the me double DecayIntegrator::threeBodyMatrixElement(const int,const Energy2, const Energy2, const Energy2,const Energy2, const Energy, const Energy, const Energy) const { throw DecayIntegratorError() << "Calling the virtual DecayIntegrator::threeBodyMatrixElement" << "method. This must be overwritten in the classes " << "inheriting from DecayIntegrator where it is needed" << Exception::runerror; } // the differential three body decay rate with one integral performed InvEnergy DecayIntegrator::threeBodydGammads(const int, const Energy2, const Energy2, const Energy, const Energy, const Energy) const { throw DecayIntegratorError() << "Calling the virtual DecayIntegrator::threeBodydGammads()" <<"method. This must be overwritten in the classes " << "inheriting from DecayIntegrator where it is needed" << Exception::runerror; } double DecayIntegrator::oneLoopVirtualME(unsigned int , const Particle &, const ParticleVector &) { throw DecayIntegratorError() << "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 DecayIntegratorError() << "DecayIntegrator::realEmmisionME() called. This should" << " have been overidden in an inheriting class if it is used" << Exception::runerror; } diff --git a/Decay/DecayIntegrator.h b/Decay/DecayIntegrator.h --- a/Decay/DecayIntegrator.h +++ b/Decay/DecayIntegrator.h @@ -1,485 +1,490 @@ // -*- C++ -*- // // DecayIntegrator.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 #include "DecayPhaseSpaceChannel.h" #include #include #include "DecayPhaseSpaceMode.fh" #include "Herwig/PDT/WidthCalculatorBase.fh" #include "Radiation/DecayRadiationGenerator.h" #include "HwDecayerBase.h" #include "DecayIntegrator.fh" 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 DecayPhaseSpaceMode 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: /** * The output operator is a friend, this is mainly for debugging */ friend ostream & operator<<(ostream &, const DecayIntegrator &); /** * and DecayPhaseMode */ friend class DecayPhaseSpaceMode; /** * Enum for the matrix element option */ enum MEOption {Initialize,Calculate,Terminate}; public: /** * Default constructor. */ DecayIntegrator() : _niter(10), _npoint(10000), _ntry(500), _generateinter(false), _imode(-1), _realME(false), _virtualME(false), _eps(ZERO) {} /** * 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. * @param maxwgt The maximum weight for the phase space integration. * @param wgts The weights of the different channels in the multichannel approach. */ void addMode(DecayPhaseSpaceModePtr mode,double maxwgt, const vector wgts) const; /** * Return the matrix element squared for a given mode and phase-space channel. * This function is purely virtual and must be implemented in classes inheriting * from DecayIntegrator. * @param ichan The channel we are calculating the matrix element for. * @param part The decaying Particle. * @param decay The particles produced in the decay. * @param opt 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 ParticleVector & decay,MEOption opt) const=0; /** * The helicity amplitude matrix element for spin correlations. */ DecayMEPtr ME() const {return _matrixelement;} /** * 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; /** * 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); /** * Finds the phase-space mode corresponding to a given decay mode * @param dm The DecayMode */ int findMode(const DecayMode & dm); /** * 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; /** * Access to the epsilon parameter */ Energy epsilonPS() const {return _eps;} 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: /** @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: /** * 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;} /** * 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); /** * Number of decay modes */ unsigned int numberModes() const {return _modes.size();} /** * Pointer to a mode */ tDecayPhaseSpaceModePtr mode(unsigned int); /** * Pointer to a mode */ tcDecayPhaseSpaceModePtr mode(unsigned int) const; /** * Get whether or not the intermediates are included */ bool generateIntermediates() const {return _generateinter;} /** * Set whether or not the intermediates are included */ void generateIntermediates(bool in) {_generateinter=in;} /** * 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; /** * 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: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object to the begining of the run phase. */ virtual void doinitrun(); //@} private: /** * Private and non-existent assignment operator. */ DecayIntegrator & operator=(const DecayIntegrator &); private: /** * Number of iterations for th initialization. */ int _niter; /** * Number of points for initialisation */ int _npoint; /** * number of attempts to generate the decay */ 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; }; /** * Output information on the DecayIntegrator for debugging purposes */ ostream & operator<<(ostream &, const DecayIntegrator &); /** * Exception for this class and those inheriting from it */ class DecayIntegratorError: public Exception {}; } #endif /* HERWIG_DecayIntegrator_H */ diff --git a/Decay/DecayPhaseSpaceChannel.cc b/Decay/DecayPhaseSpaceChannel.cc --- a/Decay/DecayPhaseSpaceChannel.cc +++ b/Decay/DecayPhaseSpaceChannel.cc @@ -1,612 +1,595 @@ // -*- C++ -*- // // DecayPhaseSpaceChannel.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 DecayPhaseSpaceChannel class. // // Author: Peter Richardson // #include "DecayPhaseSpaceChannel.h" #include "ThePEG/Utilities/DescribeClass.h" #include "DecayPhaseSpaceMode.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include using namespace Herwig; DecayPhaseSpaceChannel::DecayPhaseSpaceChannel(tcDecayPhaseSpaceModePtr in) : _mode(in) {} void DecayPhaseSpaceChannel::persistentOutput(PersistentOStream & os) const { os << _intpart << _jactype << ounit(_intmass,GeV) << ounit(_intwidth,GeV) << ounit(_intmass2,GeV2) << ounit(_intmwidth,GeV2) << _intpower << _intdau1 << _intdau2 << _intext << _mode; } void DecayPhaseSpaceChannel::persistentInput(PersistentIStream & is, int) { is >> _intpart >> _jactype >> iunit(_intmass,GeV) >> iunit(_intwidth,GeV) >> iunit(_intmass2,GeV2) >> iunit(_intmwidth,GeV2) >> _intpower >> _intdau1 >> _intdau2 >> _intext >> _mode; } - + // The following static variable is needed for the type // description system in ThePEG. -DescribeClass +DescribeClass describeHerwigDecayPhaseSpaceChannel("Herwig::DecayPhaseSpaceChannel", "Herwig.so"); void DecayPhaseSpaceChannel::Init() { static RefVector interfaceIntermediateParticles ("IntermediateParticles", "The intermediate particles in the decay chain.", &DecayPhaseSpaceChannel::_intpart, 0, false, false, true, false); - - static ParVector interfacejactype - ("Jacobian", - "The type of Jacobian to use for the intermediate particle", - &DecayPhaseSpaceChannel::_jactype, - 0, 0, 0, 0, 1, false, false, true); - - static ParVector interfaceIntermediatePower - ("IntermediatePower", - "The power to use in the Jacobian", - &DecayPhaseSpaceChannel::_intpower, - 0, 0, 0, -10, 10, false, false, true); - - static ParVector interfaceIntermediateDau1 - ("IntermediateDaughter1", - "First Daughter of the intermediate", - &DecayPhaseSpaceChannel::_intdau1, - 0, 0, 0, -10, 10, false, false, true); - - static ParVector interfaceIntermediateDau2 - ("IntermediateDaughter2", - "Second Daughter of the intermediate", - &DecayPhaseSpaceChannel::_intdau2, - 0, 0, 0, -10, 10, false, false, true); - - static ClassDocumentation documentation - ("The DecayPhaseSpaceChannel class defines a channel" - " for the multichannel integration of the phase space for a decay."); - } // generate the momenta of the external particles vector DecayPhaseSpaceChannel::generateMomenta(const Lorentz5Momentum & pin, const vector & massext) { // integers for loops unsigned int ix,iy,idau[2],iz; // storage of the momenta of the external particles vector pexternal; // and the internal particles vector pinter; // copy the momentum of the incoming particle pexternal.push_back(pin); pinter.push_back(pin); pexternal.resize(_mode->numberofParticles()); pinter.resize(_intpart.size()); // masses of the intermediate particles vector massint(_intpart.size()); massint[0]=pin.mass(); // generate all the decays in the chain Energy lower,upper,lowerb[2]; for(ix=0;ix<_intpart.size();++ix) { idau[0] = abs(_intdau1[ix]); idau[1] = abs(_intdau2[ix]); // if both decay products off-shell if(_intdau1[ix]<0&&_intdau2[ix]<0) { // lower limits on the masses of the two resonances for(iy=0;iy<2;++iy) { lowerb[iy]=ZERO; bool massless=true; for(iz=0;iz<_intext[idau[iy]].size();++iz) { if(massext[_intext[idau[iy]][iz]]!=ZERO) massless = false; lowerb[iy]+=massext[_intext[idau[iy]][iz]]; } if(massless) lowerb[iy] = _mode->epsilonPS(); } // randomize the order if(UseRandom::rnd()<0.5) { // mass of the first resonance upper = massint[ix]-lowerb[1]; lower = lowerb[0]; massint[idau[0]]=generateMass(idau[0],lower,upper); // mass of the second resonance upper = massint[ix]-massint[idau[0]]; lower = lowerb[1]; massint[idau[1]]=generateMass(idau[1],lower,upper); } else { // mass of the second resonance upper = massint[ix]-lowerb[0]; lower = lowerb[1]; massint[idau[1]]=generateMass(idau[1],lower,upper); // mass of the first resonance upper = massint[ix]-massint[idau[1]]; lower = lowerb[0]; massint[idau[0]]=generateMass(idau[0],lower,upper); } // generate the momenta of the decay products twoBodyDecay(pinter[ix],massint[idau[0]],massint[idau[1]], pinter[idau[0]],pinter[idau[1]]); } // only first off-shell else if(_intdau1[ix]<0) { // compute the limits of integration upper = massint[ix]-massext[idau[1]]; lower = ZERO; bool massless=true; for(iy=0;iy<_intext[idau[0]].size();++iy) { if(massext[_intext[idau[0]][iy]]!=ZERO) massless = false; lower+=massext[_intext[idau[0]][iy]]; } if(massless) lower = _mode->epsilonPS(); massint[idau[0]]=generateMass(idau[0],lower,upper); // generate the momenta of the decay products twoBodyDecay(pinter[ix],massint[idau[0]],massext[idau[1]], pinter[idau[0]],pexternal[idau[1]]); } // only second off-shell else if(_intdau2[ix]<0) { // compute the limits of integration upper = massint[ix]-massext[idau[0]]; lower = ZERO; bool massless=true; for(iy=0;iy<_intext[idau[1]].size();++iy) { if(massext[_intext[idau[0]][iy]]!=ZERO) massless = false; lower+=massext[_intext[idau[1]][iy]]; } if(massless) lower = _mode->epsilonPS(); massint[idau[1]]=generateMass(idau[1],lower,upper); // generate the momenta of the decay products twoBodyDecay(pinter[ix],massext[idau[0]],massint[idau[1]], pexternal[idau[0]],pinter[idau[1]]); } // both on-shell else { // generate the momenta of the decay products twoBodyDecay(pinter[ix],massext[idau[0]],massext[idau[1]], pexternal[idau[0]],pexternal[idau[1]]); } } // return the external momenta return pexternal; } // generate the weight for this channel given a phase space configuration double DecayPhaseSpaceChannel::generateWeight(const vector & output) { using Constants::pi; // integers for loops unsigned int ix,iy,idau[2],iz; // include the prefactor due to the weight of the channel double wgt=1.; // work out the masses of the intermediate particles static vector intmass; intmass.clear(); Lorentz5Momentum pinter; for(ix=0;ix<_intpart.size();++ix) { pinter=output[_intext[ix][0]]; for(iz=1;iz<_intext[ix].size();++iz) pinter+=output[_intext[ix][iz]]; pinter.rescaleMass(); intmass.push_back( pinter.mass() ); } Energy2 scale(sqr(intmass[0])); // calculate the terms for each of the decays Energy lower,upper,lowerb[2]; for(ix=0;ix<_intpart.size();++ix) { idau[0] = abs(_intdau1[ix]); idau[1] = abs(_intdau2[ix]); // if both decay products off-shell Energy pcm; if(_intdau1[ix]<0&&_intdau2[ix]<0) { // lower limits on the masses of the two resonances for(iy=0;iy<2;++iy) { lowerb[iy]=ZERO; for(iz=0;iz<_intext[idau[iy]].size();++iz) lowerb[iy]+=output[_intext[idau[iy]][iz]].mass(); } // undo effect of randomising // weight for the first order // contribution of first resonance upper = intmass[ix]-lowerb[1]; lower = lowerb[0]; InvEnergy2 wgta=massWeight(idau[0],intmass[idau[0]],lower,upper); // contribution of second resonance upper = intmass[ix]-intmass[idau[0]]; lower = lowerb[1]; InvEnergy4 wgta2 = wgta*massWeight(idau[1],intmass[idau[1]],lower,upper); // weight for the second order upper = intmass[ix]-lowerb[0]; lower = lowerb[1]; InvEnergy2 wgtb=massWeight(idau[1],intmass[idau[1]],lower,upper); upper = intmass[ix]-intmass[idau[1]]; lower = lowerb[0]; InvEnergy4 wgtb2=wgtb*massWeight(idau[0],intmass[idau[0]],lower,upper); // weight factor wgt *=0.5*sqr(scale)*(wgta2+wgtb2); // factor for the kinematics pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],intmass[idau[0]], intmass[idau[1]]); if(pcm>ZERO) wgt *= intmass[ix]*8.*pi*pi/pcm; else wgt = 0.; } // only first off-shell else if(_intdau1[ix]<0) { // compute the limits of integration upper = intmass[ix]-output[idau[1]].mass(); lower = ZERO; for(iy=0;iy<_intext[idau[0]].size();++iy) lower+=output[_intext[idau[0]][iy]].mass(); wgt *=scale*massWeight(idau[0],intmass[idau[0]],lower,upper); pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],intmass[idau[0]], output[idau[1]].mass()); if(pcm>ZERO) wgt *= intmass[ix]*8.*pi*pi/pcm; else wgt = 0.; } // only second off-shell else if(_intdau2[ix]<0) { // compute the limits of integration upper = intmass[ix]-output[idau[0]].mass(); lower = ZERO; for(iy=0;iy<_intext[idau[1]].size();++iy) lower+=output[_intext[idau[1]][iy]].mass(); wgt *=scale*massWeight(idau[1],intmass[idau[1]],lower,upper); pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],intmass[idau[1]], output[idau[0]].mass()); if(pcm>ZERO) wgt *=intmass[ix]*8.*pi*pi/pcm; else wgt=0.; } // both on-shell else { pcm = Kinematics::pstarTwoBodyDecay(intmass[ix],output[idau[1]].mass(), output[idau[0]].mass()); if(pcm>ZERO) wgt *=intmass[ix]*8.*pi*pi/pcm; else wgt = 0.; } } // finally the overall factor wgt /= pi; // return the answer return wgt; } // output the information to a stream ostream & Herwig::operator<<(ostream & os, const DecayPhaseSpaceChannel & channel) { // output of the external particles os << "Channel for the decay of " << channel._mode->externalParticles(0)->PDGName() << " -> "; for(unsigned int ix=1;ixnumberofParticles();++ix) os << channel._mode->externalParticles(ix)->PDGName() << " "; os << endl; os << "Decay proceeds in following steps "; for(unsigned int ix=0;ixPDGName() << " -> "; if(channel._intdau1[ix]>0) { os << channel._mode->externalParticles(channel._intdau1[ix])->PDGName() << "(" << channel._intdau1[ix]<< ") "; } else { os << channel._intpart[-channel._intdau1[ix]]->PDGName() << "(" << channel._intdau1[ix]<< ") "; } if(channel._intdau2[ix]>0) { os << channel._mode->externalParticles(channel._intdau2[ix])->PDGName() << "(" <PDGName() << "(" <mass()); _intwidth.push_back(_intpart[ix]->width()); _intmass2.push_back(_intpart[ix]->mass()*_intpart[ix]->mass()); _intmwidth.push_back(_intpart[ix]->mass()*_intpart[ix]->width()); } // external particles for each intermediate particle vector temp; _intext.resize(_intpart.size()); // loop over the intermediate particles for(int ix=_intpart.size()-1;ix>=0;--ix) { temp.clear(); // add the first daughter if(_intdau1[ix]>=0) { temp.push_back(_intdau1[ix]); } else { int iy = -_intdau1[ix]; vector::iterator istart=_intext[iy].begin(); vector::iterator iend=_intext[iy].end(); for(;istart!=iend;++istart) temp.push_back(*istart); } // add the second daughter if(_intdau2[ix]>=0) { temp.push_back(_intdau2[ix]); } else { int iy = -_intdau2[ix]; vector::iterator istart=_intext[iy].begin(); vector::iterator iend=_intext[iy].end(); for(;istart!=iend;++istart) temp.push_back(*istart); } _intext[ix]=temp; } // ensure intermediates either have the width set, or // can't possibly be on-shell Energy massmax; if(_mode->testOnShell()) { massmax = _mode->externalParticles(0)->mass(); for(unsigned int ix=1;ix<_mode->numberofParticles();++ix) massmax -= _mode->externalParticles(ix)->mass(); } else { massmax = _mode->externalParticles(0)->massMax(); for(unsigned int ix=1;ix<_mode->numberofParticles();++ix) massmax -= _mode->externalParticles(ix)->massMin(); } for(unsigned int ix=0;ix<_intpart.size();++ix) { if(_intwidth[ix]==ZERO && ix>0 && _jactype[ix]==0 ) { Energy massmin(ZERO); for(unsigned int iy=0;iy<_intext[ix].size();++iy) massmin += _mode->testOnShell() ? _mode->externalParticles(_intext[ix][iy])->mass() : _mode->externalParticles(_intext[ix][iy])->massMin(); // check if can be on-shell if(_intmass[ix]>=massmin&&_intmass[ix]<=massmax+massmin) { string modeout; for(unsigned int iy=0;iy<_mode->numberofParticles();++iy) { modeout += _mode->externalParticles(iy)->PDGName() + " "; } throw InitException() << "Width zero for " << _intpart[ix]->PDGName() << " in DecayPhaseSpaceChannel::doinit() " << modeout << Exception::runerror; } } } } -void DecayPhaseSpaceChannel::doinitrun() { - Interfaced::doinitrun(); +void DecayPhaseSpaceChannel::initrun() { if(!_mode->testOnShell()) return; _intmass.clear(); _intwidth.clear(); _intmass2.clear(); _intmwidth.clear(); // masses and widths of the intermediate particles for(unsigned int ix=0;ix<_intpart.size();++ix) { _intmass.push_back(_intpart[ix]->mass()); _intwidth.push_back(_intpart[ix]->width()); _intmass2.push_back(_intpart[ix]->mass()*_intpart[ix]->mass()); _intmwidth.push_back(_intpart[ix]->mass()*_intpart[ix]->width()); } // ensure intermediates either have the width set, or // can't possibly be on-shell Energy massmax = _mode->externalParticles(0)->massMax(); for(unsigned int ix=1;ix<_mode->numberofParticles();++ix) massmax -= _mode->externalParticles(ix)->massMin(); for(unsigned int ix=0;ix<_intpart.size();++ix) { if(_intwidth[ix]==0.*MeV && ix>0 && _jactype[ix]==0 ) { Energy massmin(0.*GeV); - for(unsigned int iy=0;iy<_intext[ix].size();++iy) + for(unsigned int iy=0;iy<_intext[ix].size();++iy) { massmin += _mode->externalParticles(_intext[ix][iy])->massMin(); + } // check if can be on-shell if(_intmass[ix]>=massmin&&_intmass[ix]<=massmax+massmin) { string modeout; for(unsigned int iy=0;iy<_mode->numberofParticles();++iy) { modeout += _mode->externalParticles(iy)->PDGName() + " "; } throw Exception() << "Width zero for " << _intpart[ix]->PDGName() << " in DecayPhaseSpaceChannel::doinitrun() " << modeout << Exception::runerror; } } } } // generate the final-state particles including the intermediate resonances void DecayPhaseSpaceChannel::generateIntermediates(bool cc, const Particle & in, ParticleVector & out) { // integers for the loops unsigned int ix,iz; // create the particles // incoming particle ParticleVector external; external.push_back(const_ptr_cast(&in)); // outgoing for(ix=0;ixmomentum(); for(iz=1;iz<_intext[ix].size();++iz) pinter+=external[_intext[ix][iz]]->momentum(); pinter.rescaleMass(); respart = (cc&&_intpart[ix]->CC()) ? _intpart[ix]->CC()->produceParticle(pinter) : _intpart[ix] ->produceParticle(pinter); resonance.push_back(respart); } // set up the mother daughter relations for(ix=1;ix<_intpart.size();++ix) { resonance[ix]->addChild( _intdau1[ix]<0 ? resonance[-_intdau1[ix]] : external[_intdau1[ix]]); resonance[ix]->addChild( _intdau2[ix]<0 ? resonance[-_intdau2[ix]] : external[_intdau2[ix]]); if(resonance[ix]->dataPtr()->stable()) resonance[ix]->setLifeLength(Lorentz5Distance()); } // construct the output with the particles in the first step out.push_back( _intdau1[0]>0 ? external[_intdau1[0]] : resonance[-_intdau1[0]]); out.push_back( _intdau2[0]>0 ? external[_intdau2[0]] : resonance[-_intdau2[0]]); } double DecayPhaseSpaceChannel::atanhelper_(int ires, Energy limit) { return atan2( limit*limit-_intmass2[ires], _intmwidth[ires] ); } // return the weight for a given resonance InvEnergy2 DecayPhaseSpaceChannel::massWeight(int ires, Energy moff, Energy lower,Energy upper) { InvEnergy2 wgt = ZERO; if(lower>upper) { throw DecayPhaseSpaceError() << "DecayPhaseSpaceChannel::massWeight not allowed " << ires << " " << _intpart[ires]->id() << " " << moff/GeV << " " << lower/GeV << " " << upper/GeV << Exception::eventerror; } // use a Breit-Wigner if ( _jactype[ires] == 0 ) { double rhomin = atanhelper_(ires,lower); double rhomax = atanhelper_(ires,upper) - rhomin; if ( rhomax != 0.0 ) { Energy2 moff2=moff*moff-_intmass2[ires]; wgt = _intmwidth[ires]/rhomax/(moff2*moff2+_intmwidth[ires]*_intmwidth[ires]); } else { wgt = 1./((sqr(upper)-sqr(lower))*sqr(sqr(moff)-_intmass2[ires])/ (sqr(lower)-_intmass2[ires])/(sqr(upper)-_intmass2[ires])); } } // power law else if(_jactype[ires]==1) { double rhomin = pow(sqr(lower/MeV),_intpower[ires]+1.); double rhomax = pow(sqr(upper/MeV),_intpower[ires]+1.)-rhomin; wgt = (_intpower[ires]+1.)/rhomax*pow(sqr(moff/MeV),_intpower[ires]) /MeV/MeV; } else if(_jactype[ires]==2) { wgt = 1./Constants::pi/_intmwidth[ires]; } else { throw DecayPhaseSpaceError() << "Unknown type of Jacobian in " << "DecayPhaseSpaceChannel::massWeight" << Exception::eventerror; } return wgt; } Energy DecayPhaseSpaceChannel::generateMass(int ires,Energy lower,Energy upper) { static const Energy eps=1e-9*MeV; if(lowerupper) throw DecayPhaseSpaceError() << "DecayPhaseSpaceChannel::generateMass" << " not allowed" << Exception::eventerror; if(abs(lower-upper)/(lower+upper)>2e-10) { lower +=1e-10*(lower+upper); upper -=1e-10*(lower+upper); } else return 0.5*(lower+upper); // use a Breit-Wigner if(_jactype[ires]==0) { if(_intmwidth[ires]!=ZERO) { Energy2 lower2 = sqr(lower); Energy2 upper2 = sqr(upper); double rhomin = atan2((lower2 - _intmass2[ires]),_intmwidth[ires]); double rhomax = atan2((upper2 - _intmass2[ires]),_intmwidth[ires])-rhomin; double rho = rhomin+rhomax*UseRandom::rnd(); Energy2 mass2 = max(lower2,min(upper2,_intmass2[ires]+_intmwidth[ires]*tan(rho))); if(mass2upper-1e-10*(lower+upper)) mass=upper-1e-10*(lower+upper); return mass; } void DecayPhaseSpaceChannel::twoBodyDecay(const Lorentz5Momentum & p, const Energy m1, const Energy m2, Lorentz5Momentum & p1, Lorentz5Momentum & p2 ) { static const double eps=1e-6; double ctheta,phi; Kinematics::generateAngles(ctheta,phi); Axis unitDir1=Kinematics::unitDirection(ctheta,phi); Momentum3 pstarVector; Energy min=p.mass(); if ( min >= m1 + m2 && m1 >= ZERO && m2 >= ZERO ) { pstarVector = unitDir1 * Kinematics::pstarTwoBodyDecay(min,m1,m2); } else if( m1 >= ZERO && m2 >= ZERO && (m1+m2-min)/(min+m1+m2) " << m1/GeV << ' ' << m2/GeV << Exception::eventerror; } p1 = Lorentz5Momentum(m1, pstarVector); p2 = Lorentz5Momentum(m2,-pstarVector); // boost from CM to LAB Boost bv = p.boostVector(); double gammarest = p.e()/p.mass(); p1.boost( bv , gammarest ); p2.boost( bv , gammarest ); } + +bool DecayPhaseSpaceChannel::checkKinematics() { + Energy massmax = _mode->externalParticles(0)->massMax(); + for(unsigned int ix=1;ix<_mode->numberofParticles();++ix) + massmax -= _mode->externalParticles(ix)->massMin(); + for(unsigned int ix=1;ix<_intpart.size();++ix) { + Energy massmin(0.*GeV); + for(unsigned int iy=0;iy<_intext[ix].size();++iy) { + massmin += _mode->externalParticles(_intext[ix][iy])->massMin(); + } + if(_intmass[ix]>=massmin&&_intmass[ix]<=massmax+massmin) + return false; + } + return true; +} diff --git a/Decay/DecayPhaseSpaceChannel.h b/Decay/DecayPhaseSpaceChannel.h --- a/Decay/DecayPhaseSpaceChannel.h +++ b/Decay/DecayPhaseSpaceChannel.h @@ -1,361 +1,349 @@ // -*- C++ -*- // // DecayPhaseSpaceChannel.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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_DecayPhaseSpaceChannel_H #define HERWIG_DecayPhaseSpaceChannel_H // // This is the declaration of the DecayPhaseSpaceChannel class. // #include #include #include #include "DecayPhaseSpaceChannel.fh" #include "DecayIntegrator.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Repository/UseRandom.h" #include "DecayPhaseSpaceMode.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * * This class is designed to store the information needed for a given * phase-space channel for use by the multi-channel phase space decayer * and perform the generation of the phase space for that channel. * * The decay channel is specified as a series of \f$1\to2\f$ decays to either * external particles or other intermediates. For each intermediate * the Jacobian to be used can be either a Breit-Wigner(0) or a power-law * (1). * * The class is then capable of generating a phase-space point using this * channel and computing the weight of a given point for use in a multi-channel * phase space integration using the DecayPhaseSpaceMode class. * * The class is designed so that the phase-space channels can either by specified * using the addIntermediate method directly or via the repository. * (In practice at the moment all the channels are constructed by the relevant decayers * using the former method at the moment.) * * @see DecayPhaseSpaceMode * @see DecayIntegrator * * @author Peter Richardson */ -class DecayPhaseSpaceChannel: public Interfaced { +class DecayPhaseSpaceChannel: public Base { /** * A friend output operator to allow the channel to be outputted for * debugging purposes */ friend ostream & operator<<(ostream &, const DecayPhaseSpaceChannel &); /** * DecayPhaseSpaceMode is a friend to avoid making many of the phase space * generation members public. */ friend class DecayPhaseSpaceMode; public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ DecayPhaseSpaceChannel() {} /** * Constructor with a pointer to a DecayPhaseSpaceMode. This * is the constructor which should normally be used in decayers. */ DecayPhaseSpaceChannel(tcDecayPhaseSpaceModePtr); //@} public: /** @name Set-up Members */ //@{ /** * Add a new intermediate particle * @param part A pointer to the particle data object for the intermediate. * @param jac The jacobian to be used for the generation of the particle's mass * 0 is a Breit-Wigner and 1 is a power-law * @param power The power to beb used for the mass generation if a power law * mass distribution is chosen. * @param dau1 The first daughter. If this is postive it is the \f$dau1\f$ th * outgoing particle (0 is the incoming particle), if it is negative it is the * \f$|dau1|\f$ intermediate. The intermediates are specified in the order they * are added with 0 being the incoming particle. * @param dau2 The first daughter. If this is postive it is the \f$dau2\f$ th * outgoing particle (0 is the incoming particle), if it is negative it is the * \f$|dau2|\f$ intermediate. The intermediates are specified in the order they * are added with 0 being the incoming particle. */ void addIntermediate(PDPtr part,int jac,double power,int dau1,int dau2) { _intpart.push_back(part); _jactype.push_back(jac); _intpower.push_back(power); _intdau1.push_back(dau1); _intdau2.push_back(dau2); } /** * Reset the properties of an intermediate particle. This member is used * when a Decayer is used a different value for either the mass or width * of a resonace to that in the ParticleData object. This improves the * efficiency of the integration. * @param part A pointer to the particle data object for the intermediate. * @param mass The new mass of the intermediate * @param width The new width of the intermediate. */ void resetIntermediate(tcPDPtr part,Energy mass,Energy width) { if(!part) return; int idin=part->id(); for(unsigned int ix=0;ix<_intpart.size();++ix) { if(_intpart[ix] && _intpart[ix]->id()==idin) { _intmass[ix] =mass;_intwidth[ix]=width; _intmass2[ix]=mass*mass;_intmwidth[ix]=mass*width; } } } /* * Reset the one of the daughters * @param oldp The id of the particle being reset * @param newp The id of the particle replacing it */ void resetDaughter(int oldp, int newp) { for(unsigned int ix=0;ix<_intdau1.size();++ix) { if(_intdau1[ix]==oldp) _intdau1[ix]=newp; } for(unsigned int ix=0;ix<_intdau2.size();++ix) { if(_intdau2[ix]==oldp) _intdau2[ix]=newp; } } + + /** + * Check the kinematics + */ + bool checkKinematics(); //@} protected: /** @name Phase-Space Generation Members */ //@{ /** * Generate the momenta of the external particles. This member uses the masses * of the external particles generated by the DecayPhaseMode class and the * intermediates for the channel to generate the momenta of the decay products. * @param pin The momenta of the decay products. This is outputed by the member. * @param massext The masses of the particles. This is to allow inclusion of * off-shell effects for the external particles. */ vector generateMomenta(const Lorentz5Momentum & pin, const vector & massext); /** * Generate the weight for this channel given a phase space configuration. * This member generates the weight for a given phase space configuration * and is used by the DecayPhaseSpaceMode class to compute the denominator * of the weight in the multi-channel integration. * @param output The momenta of the outgoing particles. */ double generateWeight(const vector & output); /** * Generate the final-state particles including the intermediate resonances. * This method takes the outgoing particles and adds the intermediate particles * specified by this phase-space channel. This is to allow a given set of * intermediates to be specified even when there is interference between different * intermediate states. * @param cc Whether the particles are the mode specified or its charge conjugate. * @param in The incoming particles. * @param out The outgoing particles. * */ void generateIntermediates(bool cc,const Particle & in, ParticleVector & out); /** * Calculate the momenta for a two body decay * The return value indicates success or failure. * @param p The momentum of the decaying particle * @param m1 The mass of the first decay product * @param m2 The mass of the second decay product * @param p1 The momentum of the first decay product * @param p2 The momentum of the second decay product */ void twoBodyDecay(const Lorentz5Momentum & p, const Energy m1, const Energy m2, Lorentz5Momentum & p1, Lorentz5Momentum & p2); //@} 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: +public: /** @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(); + virtual void init(); /** * Initialize this object. Called in the run phase just before * a run begins. */ - virtual void doinitrun(); + virtual void initrun(); //@} private: /** * Private and non-existent assignment operator. */ DecayPhaseSpaceChannel & operator=(const DecayPhaseSpaceChannel &); private: /** @name Mass Generation Members */ //@{ /** * Generate the mass of a resonance. * @param ires The resonance to be generated. * @param lower The lower limit on the particle's mass. * @param upper The upper limit on the particle's mass. */ Energy generateMass(int ires,Energy lower,Energy upper); /** * Return the weight for a given resonance. * @param ires The resonance to be generated. * @param moff The mass of the resonance. * @param lower The lower limit on the particle's mass. * @param upper The upper limit on the particle's mass. */ InvEnergy2 massWeight(int ires, Energy moff,Energy lower,Energy upper); //@} private: /** * pointer to the mode */ tcDecayPhaseSpaceModePtr _mode; /** * Pointers to the particle data objects of the intermediate particles */ vector _intpart; /** * The type of jacobian to be used for the intermediates. */ vector _jactype; /** * The mass of the intermediates. */ vector _intmass; /** * The width of the intermediates. */ vector _intwidth; /** * The mass squared of the intermediate particles. */ vector _intmass2; /** * The mass times the width for the intermediate particles. */ vector_intmwidth; /** * The power for the intermediate resonance. */ vector _intpower; /** * The first daughter of the intermediate resonance. */ vector _intdau1; /** * The second daughter of the intermediate resonance. */ vector _intdau2; /** * The external particles that an intermediate particle final decays to. */ vector< vector > _intext; /** * Helper function for the weight calculation. * @param ires The resonance to be generated. * @param limit The limit on the particle's mass. */ double atanhelper_(int ires, Energy limit); }; /** * write the phase space channel to a stream */ ostream & operator<<(ostream &, const DecayPhaseSpaceChannel &); /** * exception for this class and those inheriting from it */ class DecayPhaseSpaceError: public Exception {}; } #endif /* HERWIG_DecayPhaseSpaceChannel_H */ diff --git a/Decay/DecayPhaseSpaceMode.cc b/Decay/DecayPhaseSpaceMode.cc --- a/Decay/DecayPhaseSpaceMode.cc +++ b/Decay/DecayPhaseSpaceMode.cc @@ -1,534 +1,527 @@ // -*- C++ -*- // // DecayPhaseSpaceMode.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 DecayPhaseSpaceMode class. // #include "DecayPhaseSpaceMode.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/PDT/GenericWidthGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/RefVector.h" #include "DecayIntegrator.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/EventRecord/HelicityVertex.h" #include "Herwig/Decay/DecayVertex.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/Utilities/Debug.h" using namespace Herwig; using namespace ThePEG::Helicity; void DecayPhaseSpaceMode::persistentOutput(PersistentOStream & os) const { os << _integrator << _channels << _channelwgts << _maxweight << _niter << _npoint << _ntry << _extpart << _partial << _widthgen << _massgen << _testOnShell << ounit(_eps,GeV); } void DecayPhaseSpaceMode::persistentInput(PersistentIStream & is, int) { is >> _integrator >> _channels >> _channelwgts >> _maxweight >> _niter >> _npoint >> _ntry >> _extpart >> _partial >> _widthgen >> _massgen >> _testOnShell >> iunit(_eps,GeV); } // The following static variable is needed for the type // description system in ThePEG. -DescribeClass +DescribeClass describeHerwigDecayPhaseSpaceMode("Herwig::DecayPhaseSpaceMode", "Herwig.so"); void DecayPhaseSpaceMode::Init() { static ClassDocumentation documentation ("The DecayPhaseSpaceMode class contains a number of phase space" " channels for the integration of a particular decay mode"); - - static RefVector interfaceChannels - ("Channels", - "The phase space integration channels.", - &DecayPhaseSpaceMode::_channels, 0, false, false, true, true); } // flat phase space generation and weight Energy DecayPhaseSpaceMode::flatPhaseSpace(bool cc, const Particle & inpart, ParticleVector & outpart, bool onShell) const { double wgt(1.); if(outpart.empty()) { outpart.reserve(_extpart.size()-1); for(unsigned int ix=1;ix<_extpart.size();++ix) { if(cc&&_extpart[ix]->CC()) { outpart.push_back((_extpart[ix]->CC())->produceParticle()); } else { outpart.push_back(_extpart[ix]->produceParticle()); } } } // masses of the particles Energy inmass(inpart.mass()); vector mass = externalMasses(inmass,wgt,onShell); // momenta of the particles vector part(outpart.size()); // two body decay assert(outpart.size()==2); double ctheta,phi; Kinematics::generateAngles(ctheta,phi); if(! Kinematics::twoBodyDecay(inpart.momentum(), mass[1], mass[2], ctheta, phi,part[0],part[1])) throw Exception() << "Incoming mass - Outgoing mass negative in " << "DecayPhaseSpaceMode::flatPhaseSpace()" << Exception::eventerror; wgt *= Kinematics::pstarTwoBodyDecay(inmass,mass[1],mass[2])/8./Constants::pi/inmass; outpart[0]->set5Momentum(part[0]); outpart[1]->set5Momentum(part[1]); return wgt*inmass; } // initialise the phase space Energy DecayPhaseSpaceMode::initializePhaseSpace(bool init, bool onShell) { _integrator->ME(DecayMEPtr()); Energy output(ZERO); // ensure that the weights add up to one if(!_channels.empty()) { double temp=0.; for(unsigned int ix=0;ix<_channels.size();++ix) temp+=_channelwgts[ix]; for(unsigned int ix=0;ix<_channels.size();++ix) _channelwgts[ix]/=temp; } if(!init) return ZERO; // create a particle vector from the particle data one ThePEG::PPtr inpart=_extpart[0]->produceParticle(); ParticleVector particles; // now if using flat phase space _maxweight=0.; double totsum(0.),totsq(0.); InvEnergy pre=1./MeV; Energy prewid; if(_channels.empty()) { double wsum=0.,wsqsum=0.; Energy m0,mmin(ZERO); for(unsigned int ix=1;ix<_extpart.size();++ix) { mmin+=_extpart[ix]->massMin(); } for(int ix=0;ix<_npoint;++ix) { // set the mass of the decaying particle m0 = !onShell ? inpart->dataPtr()->generateMass() : inpart->dataPtr()->mass(); double wgt=0.; if(m0>mmin) { inpart->set5Momentum(Lorentz5Momentum(m0)); // compute the prefactor prewid = (_widthgen&&_partial>=0) ? _widthgen->partialWidth(_partial,inpart->mass()) : inpart->dataPtr()->width(); pre = prewid>ZERO ? 1./prewid : 1./MeV; // generate the weight for this point try { int dummy; wgt = pre*weight(false,dummy,*inpart,particles,true,onShell); } catch (Veto) { wgt=0.; } } if(wgt>_maxweight) _maxweight=wgt; wsum += wgt; wsqsum += sqr(wgt); } wsum=wsum/_npoint; wsqsum=wsqsum/_npoint-sqr(wsum); if(wsqsum<0.) wsqsum=0.; wsqsum=sqrt(wsqsum/_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 " << _extpart[0]->PDGName() << " -> "; for(unsigned int ix=1,N=_extpart.size();ixPDGName() << " "; 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 { for(int iy=0;iy<_niter;++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 m0,mmin(ZERO); for(unsigned int ix=1;ix<_extpart.size();++ix) { mmin+=_extpart[ix]->massMin(); } for(int ix=0;ix<_npoint;++ix) { m0 = !onShell ? inpart->dataPtr()->generateMass() : inpart->dataPtr()->mass(); double wgt=0.; int ichan(-1); if(m0>mmin) { inpart->set5Momentum(Lorentz5Momentum(m0)); // compute the prefactor prewid= (_widthgen&&_partial>=0) ? _widthgen->partialWidth(_partial,inpart->mass()) : inpart->dataPtr()->width(); pre = prewid>ZERO ? 1./prewid : 1./MeV; // generate the weight for this point try { wgt = pre*weight(false,ichan,*inpart,particles,true,onShell); } 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/_npoint; totsq=totsq/_npoint-sqr(totsum); if(totsq<0.) totsq=0.; totsq=sqrt(totsq/_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;ix<_channels.size();++ix) { if(nchan[ix]!=0) { wsum[ix]=wsum[ix]/nchan[ix]; wsqsum[ix]=wsqsum[ix]/nchan[ix]; if(wsqsum[ix]<0.) wsqsum[ix]=0.; wsqsum[ix]=sqrt(wsqsum[ix]/nchan[ix]); } else { wsum[ix]=0; wsqsum[ix]=0; } total+=sqrt(wsqsum[ix])*_channelwgts[ix]; } if(total>0.) { double temp; for(unsigned int ix=0;ix<_channels.size();++ix) { temp=sqrt(wsqsum[ix])*_channelwgts[ix]/total; _channelwgts[ix]=temp; } } } // factor for the weight with spin correlations _maxweight*= inpart->dataPtr()->iSpin()==1 ? 1.1 : 1.6; // ouptut the information on the initialisation Energy fact = (_widthgen&&_partial>=0) ? _widthgen->partialWidth(_partial,inpart->nominalMass()) : inpart->dataPtr()->width(); if(fact==ZERO) fact=MeV; if ( Debug::level > 1 ) { CurrentGenerator::log() << "Initialized the phase space for the decay " << _extpart[0]->PDGName() << " -> "; for(unsigned int ix=1,N=_extpart.size();ixPDGName() << " "; 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();ix momenta; double wgt(UseRandom::rnd()); // select a channel ichan=-1; do{++ichan;wgt-=_channelwgts[ichan];} while(ichan0.); // generate the momenta if(ichan==int(_channels.size())) { throw DecayIntegratorError() << "DecayPhaseSpaceMode::channelPhaseSpace()" << " failed to select a channel" << Exception::abortnow; } // generate the masses of the external particles double masswgt(1.); vector mass(externalMasses(inpart.mass(),masswgt,onShell)); momenta=_channels[ichan]->generateMomenta(inpart.momentum(),mass); // compute the denominator of the weight wgt=0.; unsigned int ix; for(ix=0;ix<_channels.size();++ix) { wgt+=_channelwgts[ix]*_channels[ix]->generateWeight(momenta); } // now we need to set the momenta of the particles // create the particles if they don't exist if(outpart.empty()) { for(ix=1;ix<_extpart.size();++ix) { if(cc&&_extpart[ix]->CC()) { outpart.push_back((_extpart[ix]->CC())->produceParticle()); } else { outpart.push_back(_extpart[ix]->produceParticle()); } } } // set up the momenta for(ix=0;ixset5Momentum(momenta[ix+1]); // return the weight return wgt!=0. ? inpart.mass()*masswgt/wgt : ZERO; } // generate the decay ParticleVector DecayPhaseSpaceMode::generate(bool intermediates,bool cc, const Particle & inpart) const { _integrator->ME(DecayMEPtr()); // compute the prefactor InvEnergy pre(1./MeV); Energy prewid; if(_widthgen&&_partial>=0) prewid=_widthgen->partialWidth(_partial,inpart.mass()); else prewid=(inpart.dataPtr()->width()); pre = prewid>ZERO ? 1./prewid : 1./MeV; // Particle vector for the output ParticleVector particles; // 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); int ncount(0),ichan; double wgt(0.); unsigned int ix; try { do { wgt=pre*weight(cc,ichan,inrest,particles,ncount==0); ++ncount; if(wgt>_maxweight) { CurrentGenerator::log() << "Resetting max weight for decay " << inrest.PDGName() << " -> "; for(ix=0;ixPDGName(); CurrentGenerator::log() << " " << _maxweight << " " << wgt << " " << inrest.mass()/MeV << "\n"; _maxweight=wgt; } } while(_maxweight*UseRandom::rnd()>wgt&&ncount<_ntry); if(ncount>=_ntry) { CurrentGenerator::log() << "The decay " << inrest.PDGName() << " -> "; for(ix=0;ixPDGName(); CurrentGenerator::log() << " " << _maxweight << " " << _ntry << " is too inefficient for the particle " << inpart << "vetoing the decay \n"; particles.clear(); throw Veto(); } } catch (Veto) { // restore the incoming particle to its original state inrest.transform(boostFromRest); throw Veto(); } // set up the vertex for spin correlations me2(-1,inrest,particles,DecayIntegrator::Terminate); const_ptr_cast(&inpart)->spinInfo(inrest.spinInfo()); constructVertex(inpart,particles); // return if intermediate particles not required if(_channelwgts.empty()||!intermediates) { for(ix=0;ixtransform(boostFromRest); } // find the intermediate particles else { // select the channel _ichannel = selectChannel(inpart,particles); for(ix=0;ixtransform(boostFromRest); // generate the particle vector _channels[_ichannel]->generateIntermediates(cc,inpart,particles); } _integrator->ME(DecayMEPtr()); return particles; } // construct the vertex for spin corrections void DecayPhaseSpaceMode::constructVertex(const Particle & inpart, const ParticleVector & decay) 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(_integrator->ME()); _integrator->ME(DecayMEPtr()); } // output info on the mode ostream & Herwig::operator<<(ostream & os, const DecayPhaseSpaceMode & decay) { os << "The mode has " << decay._channels.size() << " channels\n"; os << "This is a mode for the decay of " << decay._extpart[0]->PDGName() << " to "; for(unsigned int iz=1,N=decay._extpart.size();izPDGName() << " "; } os << "\n"; for(unsigned int ix=0;ix(_extpart[0]->widthGenerator()); for(unsigned int ix=0;ix<_extpart.size();++ix) { assert(_extpart[ix]); _massgen[ix]= dynamic_ptr_cast(_extpart[ix]->massGenerator()); } for(unsigned int ix=0;ix<_channels.size();++ix) { _channels[ix]->init(); } // set the phase-space cut off _eps = _integrator->epsilonPS(); } -void DecayPhaseSpaceMode::doinitrun() { +void DecayPhaseSpaceMode::initrun() { // update the mass and width generators if(_extpart[0]->widthGenerator()!=_widthgen) { _widthgen= dynamic_ptr_cast(_extpart[0]->widthGenerator()); } for(unsigned int ix=0;ix<_extpart.size();++ix) { if(_massgen[ix]!=_extpart[ix]->massGenerator()) _massgen[ix] = dynamic_ptr_cast(_extpart[ix]->massGenerator()); } // check the size of the weight vector is the same as the number of channels if(_channelwgts.size()!=numberChannels()) { throw Exception() << "Inconsistent number of channel weights and channels" << " in DecayPhaseSpaceMode " << Exception::abortnow; } for(unsigned int ix=0;ix<_channels.size();++ix) { _channels[ix]->initrun(); } if(_widthgen) const_ptr_cast(_widthgen)->initrun(); tcGenericWidthGeneratorPtr wtemp; for(unsigned int ix=1;ix<_extpart.size();++ix) { wtemp= dynamic_ptr_cast(_extpart[ix]->widthGenerator()); if(wtemp) const_ptr_cast(wtemp)->initrun(); } - Interfaced::doinitrun(); } // generate the masses of the external particles vector DecayPhaseSpaceMode::externalMasses(Energy inmass,double & wgt, bool onShell) const { vector mass(1,inmass); mass.reserve(_extpart.size()); vector notdone; Energy mlow(ZERO); // set masses of stable particles and limits for(unsigned int ix=1;ix<_extpart.size();++ix) { // get the mass of the particle if can't use weight if(onShell) { mass.push_back(_extpart[ix]->mass()); } else if(!_massgen[ix] || _extpart[ix]->stable()) { mass.push_back(_extpart[ix]->generateMass()); mlow+=mass[ix]; } else { mass.push_back(ZERO); notdone.push_back(ix); mlow+=max(_extpart[ix]->mass()-_extpart[ix]->widthLoCut(),ZERO); } } if(mlow>inmass) { - if(_integrator->state()==runready) { - CurrentGenerator::log() << "Decay mode " << _extpart[0]->PDGName() << " -> "; - for(unsigned int ix=1;ix<_extpart.size();++ix) { - CurrentGenerator::log() << _extpart[ix]->PDGName() << " "; - } - CurrentGenerator::log() << "is below threshold in DecayPhaseMode::externalMasses()" - << "the threshold is " << mlow/GeV - << "GeV and the parent mass is " << inmass/GeV - << " GeV\n"; - } + // if(_integrator->state()==runready) { + // CurrentGenerator::log() << "Decay mode " << _extpart[0]->PDGName() << " -> "; + // for(unsigned int ix=1;ix<_extpart.size();++ix) { + // CurrentGenerator::log() << _extpart[ix]->PDGName() << " "; + // } + // CurrentGenerator::log() << "is below threshold in DecayPhaseMode::externalMasses()" + // << "the threshold is " << mlow/GeV + // << "GeV and the parent mass is " << inmass/GeV + // << " GeV\n"; + // } throw Veto(); } // now we need to generate the masses for the particles we haven't unsigned int iloc; double wgttemp; Energy low=ZERO; for( ;!notdone.empty();) { iloc=long(UseRandom::rnd()*(notdone.size()-1)); low=max(_extpart[notdone[iloc]]->mass()-_extpart[notdone[iloc]]->widthLoCut(),ZERO); mlow-=low; mass[notdone[iloc]]= _massgen[notdone[iloc]]->mass(wgttemp,*_extpart[notdone[iloc]],low,inmass-mlow); assert(mass[notdone[iloc]]>=low&&mass[notdone[iloc]]<=inmass-mlow); wgt *= wgttemp; mlow += mass[notdone[iloc]]; notdone.erase(notdone.begin()+iloc); } return mass; } diff --git a/Decay/DecayPhaseSpaceMode.h b/Decay/DecayPhaseSpaceMode.h --- a/Decay/DecayPhaseSpaceMode.h +++ b/Decay/DecayPhaseSpaceMode.h @@ -1,471 +1,454 @@ // -*- C++ -*- // // DecayPhaseSpaceMode.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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_DecayPhaseSpaceMode_H #define HERWIG_DecayPhaseSpaceMode_H // // This is the declaration of the DecayPhaseSpaceMode class. // #include "ThePEG/Interface/Interfaced.h" #include "DecayPhaseSpaceMode.fh" #include "DecayPhaseSpaceChannel.h" #include "Herwig/PDT/GenericWidthGenerator.h" #include "Herwig/PDT/GenericMassGenerator.h" #include #include "DecayIntegrator.h" namespace Herwig { using namespace ThePEG; /** \ingroup Decay * * The DecayPhaseSpaceMode class is designed to store a group * of phase-space channels for use by the DecayIntegrator class to * generate the phase-space for a given decay mode. * * Additional phase-space channels can be added using the addChannel member. * * In practice the modes are usually constructed together with the a number of * DecayPhaseSpaceChannel objects. In classes inheriting from the * DecayIntegrator class. * * @see DecayIntegrator * @see DecayPhaseSpaceChannel * * @author Peter Richardson * */ -class DecayPhaseSpaceMode: public Interfaced { +class DecayPhaseSpaceMode: public Base { /** * A friend operator to allow the mode to be outputted for debugging purposes. */ friend ostream & operator<<(ostream &, const DecayPhaseSpaceMode &); /** * DecayIntegrator is a friend to avoid making many of the phase space * generation members public. */ friend class DecayIntegrator; /** * DecayPhaseSpaceChannel is a friend to avoid making many of the phase space * generation members public */ friend class DecayPhaseSpaceChannel; public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ DecayPhaseSpaceMode() : _maxweight(0.),_niter(10), _npoint(10000), _ntry(500), _partial(-1), _testOnShell(false), _ichannel(999) {} /** * Constructor with a pointer to a DecayPhaseIntegrator and a vector * of particle data objects for the external particles. This * is the constructor which should normally be used in decayers. * @param in The particle data objects for the external particles * @param intin A pointer to the DecayIntegrator class using this mode. * @param onShell Whether or not to perform tests for zero width on-shell particles */ DecayPhaseSpaceMode(tPDVector in, tcDecayIntegratorPtr intin,bool onShell=false) : _integrator(intin), _maxweight(0.), _niter(10), _npoint(10000), _ntry(500), _extpart(in), _partial(-1), _testOnShell(onShell), _ichannel(999) {} //@} /** * Access to the external particles. * @param ix The external particle required. * @return A pointer to the ParticleData object. */ tcPDPtr externalParticles(int ix) const {return _extpart[ix];} /** * Number of external particles. * @return The number of external particles. */ unsigned int numberofParticles() const {return _extpart.size();} /** * Number of channels * @return The number of channels. */ unsigned int numberChannels() const {return _channels.size();} /** * Add a new channel. * @param channel A pointer to the new DecayPhaseChannel */ void addChannel(DecayPhaseSpaceChannelPtr channel) { channel->init(); _channels.push_back(channel); } /** * Reset the properties of one of the intermediate particles. Only a specific channel * is reset. * @param ichan The channel to reset. * @param part The ParticleData object of the particle to reset * @param mass The mass of the intermediate. * @param width The width of gthe intermediate. */ void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width) { if(!part) return; _channels[ichan]->resetIntermediate(part,mass,width); } /** * Reset the properties of one of the intermediate particles. All the channels * are reset. * @param part The ParticleData object of the particle to reset * @param mass The mass of the intermediate. * @param width The width of gthe intermediate. */ void resetIntermediate(tcPDPtr part, Energy mass, Energy width) { for(unsigned int ix=0,N=_channels.size();ixresetIntermediate(part,mass,width); } /** * Get the maximum weight for the decay. * @return The maximum weight. */ double maxWeight() const {return _maxweight;} /** * Set the maximum weight for the decay. * @param wgt The maximum weight. */ void setMaxWeight(double wgt) const {_maxweight=wgt;} /** * Get the weight for a channel. This is the weight for the multi-channel approach. * @param ichan The channel. * @return The weight for the channel. */ double channelWeight(unsigned int ichan) const {return _channelwgts[ichan];} /** * Set the weights for the different channels. * @param in The weights for the different channels in the multi-channel approach. */ void setWeights(const vector & in) const {_channelwgts=in;} /** * Access to the selected channel */ unsigned int selectedChannel() const {return _ichannel;} /** * test on/off-shell kinematics */ bool testOnShell() const { return _testOnShell; } /** * Access to the epsilon parameter */ Energy epsilonPS() const {return _eps;} protected: /** @name Set-up, Initialization and Access Members */ //@{ /** * Initialise the phase space. * @param init Perform the initialization. */ Energy initializePhaseSpace(bool init, bool onShell=false); /** * Set the integration parameters * @param iter The number of iterations to use for initialization. * @param points The number of points to use for each iteration during initialization. * @param ntry The number of tries to generate a decay. */ void setIntegrate(int iter,int points,int ntry) { _niter=iter; _npoint=points; _ntry=ntry; } /** * Set the partial width to use for normalization. This is the partial width * in the WidthGenerator object. * @param in The partial width to use. */ void setPartialWidth(int in) {_partial=in;} //@} /** @name Phase-Space Generation Members */ //@{ /** * Access to the matrix element from the decayer. * @param ichan The channel, this is to allow the matrix element to be used to * select the intermediates * @param inpart The incoming particle. * @param opt The option for what to calculate * @param outpart The outgoing particles. */ double me2(const int ichan ,const Particle & inpart, const ParticleVector &outpart,DecayIntegrator::MEOption opt) const { return _integrator->me2(ichan,inpart,outpart,opt); } /** * Generate the decay. * @param intermediates Whether or not to generate the intermediate particle * in the decay channel. * @param cc Whether we are generating the mode specified or the charge * conjugate mode. * @param inpart The incoming particle. * @return The outgoing particles. */ ParticleVector generate(bool intermediates,bool cc,const Particle & inpart) const; /** * Select which channel to use to output the particles. * @param inpart The incoming particles. * @param outpart The outgoing particles. */ int selectChannel(const Particle & inpart, ParticleVector & outpart) const { // if using flat phase-space don't need to do this if(_channelwgts.empty()) return 0; vector mewgts(_channels.size(),0.0); double total=0.; for(unsigned int ix=0,N=_channels.size();ix0.); return ichan; } /** * Return the weight for a given phase-space point. * @param cc Whether we are generating the mode specified or the charge * conjugate mode. * @param ichan The channel to generate the weight for. * @param in The incoming particle. * @param particles The outgoing particles. * @param first Whether or not this is the first call for initialisation * @return The weight. */ Energy weight(bool cc,int & ichan, const Particle & in, ParticleVector & particles,bool first, bool onShell=false) const { ichan=0; Energy phwgt = (_channels.size()==0) ? flatPhaseSpace(cc,in,particles,onShell) : channelPhaseSpace(cc,ichan,in,particles,onShell); // generate the matrix element return me2(-1,in,particles, first ? DecayIntegrator::Initialize : DecayIntegrator::Calculate)*phwgt; } /** * Return the weight and momenta for a flat phase-space decay. * @param cc Whether we are generating the mode specified or the charge * conjugate mode. * @param inpart The incoming particle. * @param outpart The outgoing particles. * @return The weight. */ Energy flatPhaseSpace(bool cc,const Particle & inpart, ParticleVector & outpart, bool onShell=false) const; /** * Generate a phase-space point using multichannel phase space. * @param cc Whether we are generating the mode specified or the charge * conjugate mode. * @param ichan The channel to generate the weight for. * @param in The incoming particle. * @param particles The outgoing particles. * @return The weight. */ Energy channelPhaseSpace(bool cc,int & ichan, const Particle & in, ParticleVector & particles, bool onShell=false) const; /** * Construct the vertex for spin corrections * @param in The incoming particle. * @param out The outgoing particles. */ void constructVertex(const Particle & in, const ParticleVector & out) const; /** * Generate the masses of the external particles. * @param inmass The mass of the decaying particle. * @param wgt The weight for the masses. * @return The masses. */ vector externalMasses(Energy inmass,double & wgt, bool onShell) 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(); + virtual void init(); /** * Initialize this object to the begining of the run phase. */ - virtual void doinitrun(); + virtual void initrun(); //@} private: /** * Private and non-existent assignment operator. */ DecayPhaseSpaceMode & operator=(const DecayPhaseSpaceMode &); private: /** * pointer to the decayer */ tcDecayIntegratorPtr _integrator; /** * pointers to the phase-space channels */ vector _channels; /** * the weights for the different channels */ mutable vector _channelwgts; /** * the maximum weight for the decay */ mutable double _maxweight; /** * Number of iterations for the initialization. */ int _niter; /** * Number of weights for each iteration of the initialization. */ int _npoint; /** * Number of attempts to generate the decay */ int _ntry; /** * External particles */ tPDVector _extpart; /** * Which of the partial widths of the incoming particle to use */ int _partial; /** * The width generator for the incoming particle. */ cGenericWidthGeneratorPtr _widthgen; /** * The mass generators for the outgoing particles. */ vector _massgen; /** * Whether to check on-shell or off-shell kinematics * in doinit, if on-shell off-shell is tested in initrun */ bool _testOnShell; /** * The selected channel */ mutable unsigned int _ichannel; /** * Epsilon parameter for phase-space integration */ Energy _eps; }; /** * The output operator which is used for debugging. */ ostream & operator<<(ostream &, const DecayPhaseSpaceMode &); } #endif /* HERWIG_DecayPhaseSpaceMode_H */ diff --git a/Decay/Perturbative/SMTopDecayer.cc b/Decay/Perturbative/SMTopDecayer.cc --- a/Decay/Perturbative/SMTopDecayer.cc +++ b/Decay/Perturbative/SMTopDecayer.cc @@ -1,766 +1,764 @@ // -*- C++ -*- // // SMTopDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 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 SMTopDecayer class. // #include "SMTopDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "Herwig/Decay/DecayVertex.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/PDT/ThreeBodyAllOn1IntegralCalculator.h" #include "Herwig/Shower/RealEmissionProcess.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; SMTopDecayer::SMTopDecayer() : _wquarkwgt(6,0.),_wleptonwgt(3,0.), _xg_sampling(1.5), _initialenhance(1.), _finalenhance(2.3) { _wleptonwgt[0] = 0.302583; _wleptonwgt[1] = 0.301024; _wleptonwgt[2] = 0.299548; _wquarkwgt[0] = 0.851719; _wquarkwgt[1] = 0.0450162; _wquarkwgt[2] = 0.0456962; _wquarkwgt[3] = 0.859839; _wquarkwgt[4] = 3.9704e-06; _wquarkwgt[5] = 0.000489657; generateIntermediates(true); } bool SMTopDecayer::accept(tcPDPtr parent, const tPDVector & children) const { if(abs(parent->id()) != ParticleID::t) return false; int id0(0),id1(0),id2(0); for(tPDVector::const_iterator it = children.begin(); it != children.end();++it) { int id=(**it).id(),absid(abs(id)); if(absid==ParticleID::b&&double(id)/double(parent->id())>0) { id0=id; } else { switch (absid) { case ParticleID::nu_e: case ParticleID::nu_mu: case ParticleID::nu_tau: id1 = id; break; case ParticleID::eminus: case ParticleID::muminus: case ParticleID::tauminus: id2 = id; break; case ParticleID::b: case ParticleID::d: case ParticleID::s: id1 = id; break; case ParticleID::u: case ParticleID::c: id2=id; break; default : break; } } } if(id0==0||id1==0||id2==0) return false; if(double(id1)/double(id2)>0) return false; return true; } ParticleVector SMTopDecayer::decay(const Particle & parent, const tPDVector & children) const { int id1(0),id2(0); for(tPDVector::const_iterator it = children.begin(); it != children.end();++it) { int id=(**it).id(),absid=abs(id); if(absid == ParticleID::b && double(id)/double(parent.id())>0) continue; //leptons if(absid > 10 && absid%2==0) id1=absid; if(absid > 10 && absid%2==1) id2=absid; //quarks if(absid < 10 && absid%2==0) id2=absid; if(absid < 10 && absid%2==1) id1=absid; } unsigned int imode(0); if(id2 >=11 && id2<=16) imode = (id1-12)/2; else imode = id1+1+id2/2; bool cc = parent.id() == ParticleID::tbar; ParticleVector out(generate(true,cc,imode,parent)); //arrange colour flow PPtr pparent=const_ptr_cast(&parent); out[1]->incomingColour(pparent,out[1]->id()<0); ParticleVector products = out[0]->children(); if(products[0]->hasColour()) products[0]->colourNeighbour(products[1],true); else if(products[0]->hasAntiColour()) products[0]->colourNeighbour(products[1],false); return out; } void SMTopDecayer::persistentOutput(PersistentOStream & os) const { os << FFWVertex_ << FFGVertex_ << FFPVertex_ << WWWVertex_ << _wquarkwgt << _wleptonwgt << _wplus << _initialenhance << _finalenhance << _xg_sampling; } void SMTopDecayer::persistentInput(PersistentIStream & is, int) { is >> FFWVertex_ >> FFGVertex_ >> FFPVertex_ >> WWWVertex_ >> _wquarkwgt >> _wleptonwgt >> _wplus >> _initialenhance >> _finalenhance >> _xg_sampling; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigSMTopDecayer("Herwig::SMTopDecayer", "HwPerturbativeDecay.so"); void SMTopDecayer::Init() { static ClassDocumentation documentation ("This is the implementation of the SMTopDecayer which " "decays top quarks into bottom quarks and either leptons " "or quark-antiquark pairs including the matrix element for top decay", "The matrix element correction for top decay \\cite{Hamilton:2006ms}.", "%\\cite{Hamilton:2006ms}\n" "\\bibitem{Hamilton:2006ms}\n" " K.~Hamilton and P.~Richardson,\n" " ``A simulation of QCD radiation in top quark decays,''\n" " JHEP {\\bf 0702}, 069 (2007)\n" " [arXiv:hep-ph/0612236].\n" " %%CITATION = JHEPA,0702,069;%%\n"); static ParVector interfaceQuarkWeights ("QuarkWeights", "Maximum weights for the hadronic decays", &SMTopDecayer::_wquarkwgt, 6, 1.0, 0.0, 10.0, false, false, Interface::limited); static ParVector interfaceLeptonWeights ("LeptonWeights", "Maximum weights for the semi-leptonic decays", &SMTopDecayer::_wleptonwgt, 3, 1.0, 0.0, 10.0, false, false, Interface::limited); static Parameter interfaceEnhancementFactor ("InitialEnhancementFactor", "The enhancement factor for initial-state radiation in the shower to ensure" " the weight for the matrix element correction is less than one.", &SMTopDecayer::_initialenhance, 1.0, 1.0, 10000.0, false, false, Interface::limited); static Parameter interfaceFinalEnhancementFactor ("FinalEnhancementFactor", "The enhancement factor for final-state radiation in the shower to ensure" " the weight for the matrix element correction is less than one", &SMTopDecayer::_finalenhance, 1.6, 1.0, 1000.0, false, false, Interface::limited); static Parameter interfaceSamplingTopHardMEC ("SamplingTopHardMEC", "The importance sampling power for choosing an initial xg, " "to sample xg according to xg^-_xg_sampling", &SMTopDecayer::_xg_sampling, 1.5, 1.2, 2.0, false, false, Interface::limited); } double SMTopDecayer::me2(const int, const Particle & inpart, const ParticleVector & decay, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half, PDT::Spin1Half,PDT::Spin1Half))); // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(inpart.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&inpart), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&inpart), incoming); // fix rho if no correlations fixRho(_rho); } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&inpart),incoming,true); SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[1],outgoing,true); SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[2],outgoing,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&inpart),incoming,true); SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[1],outgoing,true); SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[2],outgoing,true); } } if ( ( decay[1]->momentum() + decay[2]->momentum() ).m() < decay[1]->data().constituentMass() + decay[2]->data().constituentMass() ) return 0.0; // spinors for the decay product if(inpart.id()>0) { SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar ,decay[0],outgoing); SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[1],outgoing); SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[2],outgoing); } else { SpinorWaveFunction ::calculateWaveFunctions(_inHalf ,decay[0],outgoing); SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[1],outgoing); SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[2],outgoing); } Energy2 scale(sqr(inpart.mass())); if(inpart.id() == ParticleID::t) { //Define intermediate vector wave-function for Wplus tcPDPtr Wplus(getParticleData(ParticleID::Wplus)); VectorWaveFunction inter; unsigned int thel,bhel,fhel,afhel; for(thel = 0;thel<2;++thel){ for(bhel = 0;bhel<2;++bhel){ inter = FFWVertex_->evaluate(scale,1,Wplus,_inHalf[thel], _inHalfBar[bhel]); for(afhel=0;afhel<2;++afhel){ for(fhel=0;fhel<2;++fhel){ (*ME())(thel,bhel,afhel,fhel) = FFWVertex_->evaluate(scale,_outHalf[afhel], _outHalfBar[fhel],inter); } } } } } else if(inpart.id() == ParticleID::tbar) { VectorWaveFunction inter; tcPDPtr Wminus(getParticleData(ParticleID::Wminus)); unsigned int tbhel,bbhel,afhel,fhel; for(tbhel = 0;tbhel<2;++tbhel){ for(bbhel = 0;bbhel<2;++bbhel){ inter = FFWVertex_-> evaluate(scale,1,Wminus,_inHalf[bbhel],_inHalfBar[tbhel]); for(afhel=0;afhel<2;++afhel){ for(fhel=0;fhel<2;++fhel){ (*ME())(tbhel,bbhel,fhel,afhel) = FFWVertex_->evaluate(scale,_outHalf[afhel], _outHalfBar[fhel],inter); } } } } } double output = (ME()->contract(_rho)).real(); if(abs(decay[1]->id())<=6) output *=3.; return output; } void SMTopDecayer::doinit() { PerturbativeDecayer::doinit(); //get vertices from SM object tcHwSMPtr hwsm = dynamic_ptr_cast(standardModel()); if(!hwsm) throw InitException() << "Must have Herwig::StandardModel in " << "SMTopDecayer::doinit()"; FFWVertex_ = hwsm->vertexFFW(); FFGVertex_ = hwsm->vertexFFG(); FFPVertex_ = hwsm->vertexFFP(); WWWVertex_ = hwsm->vertexWWW(); //initialise FFWVertex_->init(); FFGVertex_->init(); FFPVertex_->init(); WWWVertex_->init(); //set up decay modes _wplus = getParticleData(ParticleID::Wplus); DecayPhaseSpaceModePtr mode; DecayPhaseSpaceChannelPtr Wchannel; tPDVector extpart(4); vector wgt(1,1.0); extpart[0] = getParticleData(ParticleID::t); extpart[1] = getParticleData(ParticleID::b); //lepton modes for(int i=11; i<17;i+=2) { extpart[2] = getParticleData(-i); extpart[3] = getParticleData(i+1); mode = new_ptr(DecayPhaseSpaceMode(extpart,this)); Wchannel = new_ptr(DecayPhaseSpaceChannel(mode)); Wchannel->addIntermediate(extpart[0],0,0.0,-1,1); Wchannel->addIntermediate(_wplus,0,0.0,2,3); - Wchannel->init(); mode->addChannel(Wchannel); addMode(mode,_wleptonwgt[(i-11)/2],wgt); } //quark modes unsigned int iz=0; for(int ix=1;ix<6;ix+=2) { for(int iy=2;iy<6;iy+=2) { // check that the combination of particles is allowed if(FFWVertex_->allowed(-ix,iy,ParticleID::Wminus)) { extpart[2] = getParticleData(-ix); extpart[3] = getParticleData( iy); mode = new_ptr(DecayPhaseSpaceMode(extpart,this)); Wchannel = new_ptr(DecayPhaseSpaceChannel(mode)); Wchannel->addIntermediate(extpart[0],0,0.0,-1,1); Wchannel->addIntermediate(_wplus,0,0.0,2,3); - Wchannel->init(); mode->addChannel(Wchannel); addMode(mode,_wquarkwgt[iz],wgt); ++iz; } else { throw InitException() << "SMTopDecayer::doinit() the W vertex" << "cannot handle all the quark modes" << Exception::abortnow; } } } } void SMTopDecayer::dataBaseOutput(ofstream & os,bool header) const { if(header) os << "update decayers set parameters=\""; // parameters for the PerturbativeDecayer base class for(unsigned int ix=0;ix<_wquarkwgt.size();++ix) { os << "newdef " << name() << ":QuarkWeights " << ix << " " << _wquarkwgt[ix] << "\n"; } for(unsigned int ix=0;ix<_wleptonwgt.size();++ix) { os << "newdef " << name() << ":LeptonWeights " << ix << " " << _wleptonwgt[ix] << "\n"; } PerturbativeDecayer::dataBaseOutput(os,false); if(header) os << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } void SMTopDecayer::doinitrun() { PerturbativeDecayer::doinitrun(); if(initialize()) { for(unsigned int ix=0;ixmaxWeight(); else _wquarkwgt [ix-3] = mode(ix)->maxWeight(); } } } WidthCalculatorBasePtr SMTopDecayer::threeBodyMEIntegrator(const DecayMode & dm) const { // identify W decay products int sign = dm.parent()->id() > 0 ? 1 : -1; int iferm(0),ianti(0); for(ParticleMSet::const_iterator pit=dm.products().begin(); pit!=dm.products().end();++pit) { int id = (**pit).id(); if(id*sign != ParticleID::b) { if (id*sign > 0 ) iferm = id*sign; else ianti = id*sign; } } assert(iferm!=0&&ianti!=0); // work out which mode we are doing int imode(-1); for(unsigned int ix=0;ixexternalParticles(2)->id() == ianti && mode(ix)->externalParticles(3)->id() == iferm ) { imode = ix; break; } } assert(imode>=0); // get the masses we need Energy m[3] = {mode(imode)->externalParticles(1)->mass(), mode(imode)->externalParticles(3)->mass(), mode(imode)->externalParticles(2)->mass()}; return new_ptr(ThreeBodyAllOn1IntegralCalculator (3,_wplus->mass(),_wplus->width(),0.0,*this,imode,m[0],m[1],m[2])); } InvEnergy SMTopDecayer::threeBodydGammads(const int imode, const Energy2 mt2, const Energy2 mffb2, const Energy mb, const Energy mf, const Energy mfb) const { Energy mffb(sqrt(mffb2)); Energy mw(_wplus->mass()); Energy2 mw2(sqr(mw)),gw2(sqr(_wplus->width())); Energy mt(sqrt(mt2)); Energy Eb = 0.5*(mt2-mffb2-sqr(mb))/mffb; Energy Ef = 0.5*(mffb2-sqr(mfb)+sqr(mf))/mffb; Energy Ebm = sqrt(sqr(Eb)-sqr(mb)); Energy Efm = sqrt(sqr(Ef)-sqr(mf)); Energy2 upp = sqr(Eb+Ef)-sqr(Ebm-Efm); Energy2 low = sqr(Eb+Ef)-sqr(Ebm+Efm); InvEnergy width=(dGammaIntegrand(mffb2,upp,mt,mb,mf,mfb,mw)- dGammaIntegrand(mffb2,low,mt,mb,mf,mfb,mw)) /32./mt2/mt/8/pow(Constants::pi,3)/(sqr(mffb2-mw2)+mw2*gw2); // couplings width *= 0.25*sqr(4.*Constants::pi*generator()->standardModel()->alphaEM(mt2)/ generator()->standardModel()->sin2ThetaW()); width *= generator()->standardModel()->CKM(*mode(imode)->externalParticles(0), *mode(imode)->externalParticles(1)); if(abs(mode(imode)->externalParticles(2)->id())<=6) { width *=3.; if(abs(mode(imode)->externalParticles(2)->id())%2==0) width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(2), *mode(imode)->externalParticles(3)); else width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(3), *mode(imode)->externalParticles(2)); } // final spin average assert(!std::isnan(double(width*MeV))); return 0.5*width; } Energy6 SMTopDecayer::dGammaIntegrand(Energy2 mffb2, Energy2 mbf2, Energy mt, Energy mb, Energy mf, Energy mfb, Energy mw) const { Energy2 mt2(sqr(mt)) ,mb2(sqr(mb)) ,mf2(sqr(mf )),mfb2(sqr(mfb )),mw2(sqr(mw )); Energy4 mt4(sqr(mt2)),mb4(sqr(mb2)),mf4(sqr(mf2)),mfb4(sqr(mfb2)),mw4(sqr(mw2)); return -mbf2 * ( + 6 * mb2 * mf2 * mfb2 * mffb2 + 6 * mb2 * mt2 * mfb2 * mffb2 + 6 * mb2 * mt2 * mf2 * mffb2 + 12 * mb2 * mt2 * mf2 * mfb2 - 3 * mb2 * mfb4 * mffb2 + 3 * mb2 * mf2 * mffb2 * mffb2 - 3 * mb2 * mf4 * mffb2 - 6 * mb2 * mt2 * mfb4 - 6 * mb2 * mt2 * mf4 - 3 * mb4 * mfb2 * mffb2 - 3 * mb4 * mf2 * mffb2 - 6 * mb4 * mf2 * mfb2 + 3 * mt4 * mf4 + 3 * mb4 * mfb4 + 3 * mb4 * mf4 + 3 * mt4 * mfb4 + 3 * mb2 * mfb2 * mffb2 * mffb2 + 3 * mt2 * mfb2 * mffb2 * mffb2 - 3 * mt2 * mfb4 * mffb2 + 3 * mt2 * mf2 * mffb2 * mffb2 - 3 * mt2 * mf4 * mffb2 - 3 * mt4 * mfb2 * mffb2 - 3 * mt4 * mf2 * mffb2 - 6 * mt4 * mf2 * mfb2 + 6 * mt2 * mf2 * mfb2 * mffb2 + 12 * mt2 * mf2 * mw4 + 12 * mb2 * mfb2 * mw4 + 12 * mb2 * mt2 * mw4 + 6 * mw2 * mt2 * mfb2 * mbf2 - 12 * mw2 * mt2 * mf2 * mffb2 - 6 * mw2 * mt2 * mf2 * mbf2 - 12 * mw2 * mt2 * mf2 * mfb2 - 12 * mw2 * mb2 * mfb2 * mffb2 - 6 * mw2 * mb2 * mfb2 * mbf2 + 6 * mw2 * mb2 * mf2 * mbf2 - 12 * mw2 * mb2 * mf2 * mfb2 - 12 * mw2 * mb2 * mt2 * mfb2 - 12 * mw2 * mb2 * mt2 * mf2 + 12 * mf2 * mfb2 * mw4 + 4 * mbf2 * mbf2 * mw4 - 6 * mfb2 * mbf2 * mw4 - 6 * mf2 * mbf2 * mw4 - 6 * mt2 * mbf2 * mw4 - 6 * mb2 * mbf2 * mw4 + 12 * mw2 * mt2 * mf4 + 12 * mw2 * mt4 * mf2 + 12 * mw2 * mb2 * mfb4 + 12 * mw2 * mb4 * mfb2) /mw4 / 3.; } void SMTopDecayer::initializeMECorrection(RealEmissionProcessPtr born, double & initial, double & final) { // check the outgoing particles PPtr part[2]; for(unsigned int ix=0;ixbornOutgoing().size();++ix) { part[ix]= born->bornOutgoing()[ix]; } // check the final-state particles and get the masses if(abs(part[0]->id())==ParticleID::Wplus&&abs(part[1]->id())==ParticleID::b) { _ma=part[0]->mass(); _mc=part[1]->mass(); } else if(abs(part[1]->id())==ParticleID::Wplus&&abs(part[0]->id())==ParticleID::b) { _ma=part[1]->mass(); _mc=part[0]->mass(); } else { return; } // set the top mass _mt=born->bornIncoming()[0]->mass(); // set the gluon mass _mg=getParticleData(ParticleID::g)->constituentMass(); // set the radiation enhancement factors initial = _initialenhance; final = _finalenhance; // reduced mass parameters _a=sqr(_ma/_mt); _g=sqr(_mg/_mt); _c=sqr(_mc/_mt); double lambda = sqrt(1.+sqr(_a)+sqr(_c)-2.*_a-2.*_c-2.*_a*_c); _ktb = 0.5*(3.-_a+_c+lambda); _ktc = 0.5*(1.-_a+3.*_c+lambda); useMe(); } bool SMTopDecayer::softMatrixElementVeto(PPtr parent, PPtr progenitor, const bool & , const Energy & highestpT, const vector &, const double & z, const Energy & scale, const Energy & pt) { // check if we need to apply the full correction // the initial-state correction if(abs(progenitor->id())==ParticleID::t&&abs(parent->id())==ParticleID::t) { // check if hardest so far // if not just need to remove effect of enhancement bool veto(false); // if not hardest so far if(ptxgbcut(_ktb)) wgt = 0.; if(wgt>1.) { generator()->log() << "Violation of maximum for initial-state " << " soft veto in " << "SMTopDecayer::softMatrixElementVeto" << "xg = " << xg << " xa = " << xa << "weight = " << wgt << "\n"; wgt=1.; } // compute veto from weight veto = !UseRandom::rndbool(wgt); } } // return the veto return veto; } // final-state correction else if(abs(progenitor->id())==ParticleID::b&&abs(parent->id())==ParticleID::b) { // check if hardest so far // if not just need to remove effect of enhancement // if not hardest so far if(ptlog() << "Imaginary root for final-state veto in " << "SMTopDecayer::softMatrixElementVeto" << "\nz = " << z << "\nkappa = " << kappa << "\nxa = " << xa << "\nroot^2= " << root; return true; } root=sqrt(root); double xg((2.-xa)*(1.-r)-(z-r)*root); // xfact (below) is supposed to equal xg/(1-z). double xfact(z*kappa/2./(z*(1.-z)*kappa+_c)*(2.-xa-root)+root); // calculate the full result double f(me(xa,xg)); // jacobian double J(z*root); double wgt(f*J*2.*kappa/(1.+sqr(z)-2.*_c/kappa/z)/sqr(xfact)/_finalenhance); if(wgt>1.) { generator()->log() << "Violation of maximum for final-state soft veto in " << "SMTopDecayer::softMatrixElementVeto" << "xg = " << xg << " xa = " << xa << "weight = " << wgt << "\n"; wgt=1.; } // compute veto from weight and return return !UseRandom::rndbool(wgt); } // otherwise don't veto else return !UseRandom::rndbool(1./_finalenhance); } double SMTopDecayer::me(double xw,double xg) { double prop(1.+_a-_c-xw),xg2(sqr(xg)); double lambda=sqrt(1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c); double denom=(1.-2*_a*_a+_a+_c*_a+_c*_c-2.*_c); double wgt=-_c*xg2/prop+(1.-_a+_c)*xg-(prop*(1 - xg)+xg2) +(0.5*(1.+2.*_a+_c)*sqr(prop-xg)*xg+2.*_a*prop*xg2)/denom; return wgt/(lambda*prop); } // xgbcut is the point along the xg axis where the upper bound on the // top quark (i.e. b) emission phase space goes back on itself in the // xa vs xg plane i.e. roughly mid-way along the xg axis in // the xa vs xg Dalitz plot. double SMTopDecayer::xgbcut(double kt) { double lambda2 = 1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c; double num1 = kt*kt*(1.-_a-_c); double num2 = 2.*kt*sqrt(_a*(kt*kt*_c+lambda2*(kt-1.))); return (num1-num2)/(kt*kt-4.*_a*(kt-1.)); } double SMTopDecayer::loME(const Particle & inpart, const ParticleVector & decay) { // spinors vector swave; vector awave; vector vwave; tPPtr Wboson = abs(decay[0]->id())==ParticleID::Wplus ? decay[0] : decay[1]; tPPtr bquark = abs(decay[0]->id())==ParticleID::Wplus ? decay[1] : decay[0]; // spinors if(inpart.id()>0) { SpinorWaveFunction ::calculateWaveFunctions(swave,const_ptr_cast(&inpart), incoming); SpinorBarWaveFunction::calculateWaveFunctions(awave,bquark,outgoing); } else { SpinorBarWaveFunction::calculateWaveFunctions(awave,const_ptr_cast(&inpart), incoming); SpinorWaveFunction ::calculateWaveFunctions(swave,bquark,outgoing); } // polarization vectors VectorWaveFunction::calculateWaveFunctions(vwave,Wboson,outgoing,false); Energy2 scale(sqr(inpart.mass())); double me=0.; if(inpart.id() == ParticleID::t) { for(unsigned int thel = 0; thel < 2; ++thel) { for(unsigned int bhel = 0; bhel < 2; ++bhel) { for(unsigned int whel = 0; whel < 3; ++whel) { Complex diag = FFWVertex_->evaluate(scale,swave[thel],awave[bhel],vwave[whel]); me += norm(diag); } } } } else if(inpart.id() == ParticleID::tbar) { for(unsigned int thel = 0; thel < 2; ++thel) { for(unsigned int bhel = 0; bhel < 2; ++bhel){ for(unsigned int whel = 0; whel < 3; ++whel) { Complex diag = FFWVertex_->evaluate(scale,swave[bhel],awave[thel],vwave[whel]); me += norm(diag); } } } } return me; } double SMTopDecayer::realME(const Particle & inpart, const ParticleVector & decay, ShowerInteraction inter) { // vertex for emission from fermions AbstractFFVVertexPtr vertex = inter==ShowerInteraction::QCD ? FFGVertex_ : FFPVertex_; // spinors vector swave; vector awave; vector vwave,gwave; tPPtr Wboson = abs(decay[0]->id())==ParticleID::Wplus ? decay[0] : decay[1]; tPPtr bquark = abs(decay[0]->id())==ParticleID::Wplus ? decay[1] : decay[0]; // spinors if(inpart.id()>0) { SpinorWaveFunction ::calculateWaveFunctions(swave,const_ptr_cast(&inpart), incoming); SpinorBarWaveFunction::calculateWaveFunctions(awave,bquark,outgoing); } else { SpinorBarWaveFunction::calculateWaveFunctions(awave,const_ptr_cast(&inpart), incoming); SpinorWaveFunction ::calculateWaveFunctions(swave,bquark,outgoing); } // polarization vectors VectorWaveFunction::calculateWaveFunctions(vwave,Wboson,outgoing,false); VectorWaveFunction::calculateWaveFunctions(gwave,decay[2],outgoing,true ); Energy2 scale(sqr(inpart.mass())); double me=0.; vector diag(3,0.); if(inpart.id() == ParticleID::t) { for(unsigned int thel = 0; thel < 2; ++thel) { for(unsigned int bhel = 0; bhel < 2; ++bhel) { for(unsigned int whel = 0; whel < 3; ++whel) { for(unsigned int ghel =0; ghel <3; ghel+=2) { // emission from top SpinorWaveFunction interF = vertex->evaluate(scale,3,inpart.dataPtr(),swave[thel],gwave[ghel]); diag[0] = FFWVertex_->evaluate(scale,interF,awave[bhel],vwave[whel]); // emission from bottom SpinorBarWaveFunction interB = vertex->evaluate(scale,3,bquark->dataPtr()->CC(),awave[bhel],gwave[ghel]); diag[1] = FFWVertex_->evaluate(scale,swave[thel],interB,vwave[whel]); // emission from W if(inter==ShowerInteraction::QED) { VectorWaveFunction interV = WWWVertex_->evaluate(scale,3,Wboson->dataPtr()->CC(),vwave[whel],gwave[ghel]); diag[1] = FFWVertex_->evaluate(scale,swave[thel],awave[bhel],interV); } Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.)); me += norm(sum); } } } } } else if(inpart.id() == ParticleID::tbar) { for(unsigned int thel = 0; thel < 2; ++thel) { for(unsigned int bhel = 0; bhel < 2; ++bhel){ for(unsigned int whel = 0; whel < 3; ++whel) { for(unsigned int ghel =0; ghel <3; ghel+=2) { // emission from top SpinorBarWaveFunction interB = vertex->evaluate(scale,3,inpart.dataPtr(),awave[thel],gwave[ghel]); diag[1] = FFWVertex_->evaluate(scale,swave[bhel],interB,vwave[whel]); // emission from bottom SpinorWaveFunction interF = vertex->evaluate(scale,3,bquark->dataPtr()->CC(),swave[bhel],gwave[ghel]); diag[0] = FFWVertex_->evaluate(scale,interF,awave[thel],vwave[whel]); // emission from W if(inter==ShowerInteraction::QED) { VectorWaveFunction interV = WWWVertex_->evaluate(scale,3,Wboson->dataPtr()->CC(),vwave[whel],gwave[ghel]); diag[1] = FFWVertex_->evaluate(scale,swave[bhel],awave[thel],interV); } Complex sum = std::accumulate(diag.begin(),diag.end(),Complex(0.)); me += norm(sum); } } } } } // divide out the coupling me /= norm(vertex->norm()); // return the total return me; } double SMTopDecayer::matrixElementRatio(const Particle & inpart, const ParticleVector & decay2, const ParticleVector & decay3, MEOption , ShowerInteraction inter) { double Nc = standardModel()->Nc(); double Cf = (sqr(Nc) - 1.) / (2.*Nc); // if(inter==ShowerInteraction::QED) return 0.; // double f = (1. + sqr(e2()) - 2.*sqr(s2()) + s2() + s2()*e2() - 2.*e2()); // // // double B = f/s2(); // Energy2 PbPg = decay3[0]->momentum()*decay3[2]->momentum(); // Energy2 PtPg = inpart.momentum()*decay3[2]->momentum(); // Energy2 PtPb = inpart.momentum()*decay3[0]->momentum(); // double R = Cf *((-4.*sqr(mb())*f/s2()) * ((sqr(mb())*e2()/sqr(PbPg)) + // (sqr(mb())/sqr(PtPg)) - 2.*(PtPb/(PtPg*PbPg))) + // (16. + 8./s2() + 8.*e2()/s2()) * ((PtPg/PbPg) + (PbPg/PtPg)) - // (16./s2()) * (1. + e2())); // return R/B*Constants::pi; double Bnew = loME(inpart,decay2); double Rnew = realME(inpart,decay3,inter); double output = Rnew/Bnew*4.*Constants::pi*sqr(inpart.mass())*UnitRemoval::InvE2; if(inter==ShowerInteraction::QCD) output *= Cf; return output; }