Page MenuHomeHEPForge
Files F7879112

No OneTemporary
Actions

View Options

diff --git a/Decay/DecayIntegrator2.h b/Decay/DecayIntegrator2.h
--- a/Decay/DecayIntegrator2.h
+++ b/Decay/DecayIntegrator2.h
@@ -1,499 +1,504 @@
// -*- C++ -*-
#ifndef Herwig_DecayIntegrator2_H
#define Herwig_DecayIntegrator2_H
//
// This is the declaration of the DecayIntegrator2 class.
//
#include "DecayIntegrator2.fh"
#include "HwDecayerBase.h"
#include "PhaseSpaceMode.fh"
#include "Herwig/PDT/WidthCalculatorBase.fh"
#include "Radiation/DecayRadiationGenerator.h"
#include <Herwig/Decay/DecayVertex.h>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Decay
* \class DecayIntegrator2
* \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 <code>DecayIntegrator2</code> class inherits from ThePEG's Decayer class
* and makes use of the <code>PhaseSpaceMode</code> 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 DecayIntegrator2Interfaces "The interfaces"
* defined for DecayIntegrator2.
*/
class DecayIntegrator2: public HwDecayerBase {
public:
/**
* and DecayPhaseMode
*/
friend class PhaseSpaceMode;
/**
* Enum for the matrix element option
*/
enum MEOption {Initialize,Calculate,Terminate};
public:
/**
* The default constructor.
*/
DecayIntegrator2() : nIter_(10), nPoint_(10000), nTry_(500),
generateInter_(false), iMode_(-1),
realME_(false), virtualME_(false), eps_(ZERO)
{}
public:
/**
* Check if this decayer can perfom the decay for a particular mode.
* Uses the modeNumber member but can be overridden
* @param parent The decaying particle
* @param children The decay products
*/
virtual bool accept(tcPDPtr parent, const tPDVector & children) const {
bool cc;
return modeNumber(cc,parent,children)>=0;
}
/**
* For a given decay mode and a given particle instance, perform the
* decay and return the decay products. As this is the base class this
* is not implemented.
* @return The vector of particles produced in the decay.
*/
virtual ParticleVector decay(const Particle & parent,
const tPDVector & children) const;
/**
* Which of the possible decays is required
* @param cc Is this mode the charge conjugate
* @param parent The decaying particle
* @param children The decay products
*/
virtual int modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const = 0;
/**
* The mode being used for this decay
*/
int imode() const {return iMode_;}
/**
* Add a phase-space mode to the list
* @param mode The mode being added.
*/
void addMode(PhaseSpaceModePtr mode) const;
/**
* Return the matrix element squared for a given mode and phase-space channel.
* @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
* @param outgoing The particles produced in the decay
* @param momenta The momenta of the particles produced in the decay
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
virtual double me2(const int ichan, const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const = 0;
/**
* Construct the SpinInfos for the particles produced in the decay
*/
virtual void constructSpinInfo(const Particle & part,
ParticleVector outgoing) const = 0;
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
*/
virtual void dataBaseOutput(ofstream & os,bool header) const;
/**
* Set the code for the partial width. Finds the partial width in the
* GenericWidthGenerator class which corresponds to the decay mode.
* @param dm The DecayMode
* @param imode The mode.
*/
void setPartialWidth(const DecayMode & dm, int imode);
/**
* Specify the \f$1\to2\f$ matrix element to be used in the running width calculation.
* @param mecode The code for the matrix element as described
* in the GenericWidthGenerator class.
* @param coupling The coupling for the matrix element.
* @return True or False if this mode can be handled.
*/
virtual bool twoBodyMEcode(const DecayMode &, int & mecode,
double & coupling) const {
coupling = 1.;
mecode = -1;
return false;
}
/**
* Method to return an object to calculate the 3 (or higher body) partial width
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
/**
* The matrix element to be integrated for the three-body decays as a function
* of the invariant masses of pairs of the outgoing particles.
* @param imode The mode for which the matrix element is needed.
* @param q2 The scale, \e i.e. the mass squared of the decaying particle.
* @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
* @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
* @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
* @param m1 The mass of the first outgoing particle.
* @param m2 The mass of the second outgoing particle.
* @param m3 The mass of the third outgoing particle.
* @return The matrix element
*/
virtual double threeBodyMatrixElement(const int imode, const Energy2 q2,
const Energy2 s3, const Energy2 s2,
const Energy2 s1, const Energy m1,
const Energy m2, const Energy m3) const;
/**
* The differential three body decay rate with one integral performed.
* @param imode The mode for which the matrix element is needed.
* @param q2 The scale, \e i.e. the mass squared of the decaying particle.
* @param s The invariant mass which still needs to be integrate over.
* @param m1 The mass of the first outgoing particle.
* @param m2 The mass of the second outgoing particle.
* @param m3 The mass of the third outgoing particle.
* @return The differential rate \f$\frac{d\Gamma}{ds}\f$
*/
virtual InvEnergy threeBodydGammads(const int imode, const Energy2 q2,
const Energy2 s,
const Energy m1, const Energy m2,
const Energy m3) const;
/**
* Finds the phase-space mode corresponding to a given decay mode
* @param dm The DecayMode
*/
int findMode(const DecayMode & dm);
public:
/**
* Members for the generation of QED radiation in the decays
*/
//@{
/**
* Use the DecayRadiationGenerator to generate photons in the decay.
* @param p The Particle instance being decayed
* @param children The decay products
* @return A particle vector containing the decay products after the generation
* of photons.
*/
ParticleVector generatePhotons(const Particle & p,ParticleVector children) {
assert(false);
//return photonGen_->generatePhotons(p,children,this);
}
/**
* check if photons can be generated in the decay
*/
bool canGeneratePhotons() {return photonGen_;}
/**
* The one-loop virtual correction.
* @param imode The mode required.
* @param part The decaying particle.
* @param products The decay products including the radiated photon.
* @return Whether the correction is implemented
*/
virtual double oneLoopVirtualME(unsigned int imode,
const Particle & part,
const ParticleVector & products);
/**
* Whether or not the one loop matrix element is implemented
*/
bool hasOneLoopME() {return virtualME_;}
/**
* The real emission matrix element
* @param imode The mode required
* @param part The decaying particle
* @param products The decay products including the radiated photon
* @param iemitter The particle which emitted the photon
* @param ctheta The cosine of the polar angle between the photon and the
* emitter
* @param stheta The sine of the polar angle between the photon and the
* emitter
* @param rot1 Rotation from rest frame to frame for real emission
* @param rot2 Rotation to place emitting particle along z
*/
virtual InvEnergy2 realEmissionME(unsigned int imode,
const Particle & part,
ParticleVector & products,
unsigned int iemitter,
double ctheta, double stheta,
const LorentzRotation & rot1,
const LorentzRotation & rot2);
/**
* Whether or not the real emission matrix element is implemented
*/
bool hasRealEmissionME() {return realME_;}
//@}
public:
/**
* The output operator is a friend, this is mainly for debugging
*/
friend ostream & operator<<(ostream & os, const DecayIntegrator2 & decay);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/**
* Generate the momenta for the decay
* @param inter Generate the intermediates produced in the decay as well as the
* final particles.
* @param cc Is this the mode defined or its charge conjugate.
* @param imode The mode being generated.
* @param inpart The decaying particle.
* @return The particles produced inthe decay.
*/
ParticleVector generate(bool inter,bool cc, const unsigned int & imode,
const Particle & inpart) const;
/**
* Set the mode being use for this decay.
*/
void imode(int in) { iMode_ = in;}
/**
* Set the helicity matrix element for the decay.
*/
void ME(DecayMEPtr in) const { matrixElement_ = in;}
/**
* The helicity amplitude matrix element for spin correlations.
*/
DecayMEPtr ME() const {return matrixElement_;}
/**
* Reset the properities of all intermediates.
* @param part The intermediate particle being reset.
* @param mass The mass of the particle.
* @param width The width of the particle.
*/
void resetIntermediate(tcPDPtr part, Energy mass, Energy width);
/**
* Initialize the phase-space mode
* @param imode The mode
* @param init Whether or not to perform the initialization
*/
Energy initializePhaseSpaceMode(unsigned int imode,bool init, bool onShell=false) const;
protected:
/**
* Methods to set variables in inheriting classes
*/
//@{
/**
* Set whether or not the intermediates are included
*/
void generateIntermediates(bool in) {generateInter_=in;}
+
+ /**
+ * Set whether or not the intermediates are included
+ */
+ bool generateIntermediates() const {return generateInter_;}
/**
* Whether or not the one loop matrix element is implemented
*/
void hasOneLoopME(bool in) {virtualME_=in;}
/**
* Whether or not the real emission matrix element is implemented
*/
void hasRealEmissionME(bool in) {realME_=in;}
/**
* Set the epsilon parameter
*/
void epsilonPS(Energy in) {eps_=in;}
//@}
/**
* Number of decay modes
*/
unsigned int numberModes() const {return modes_.size();}
/**
* Pointer to a mode
*/
tPhaseSpaceModePtr mode(unsigned int ix) {
return modes_[ix];
}
/**
* Pointer to a mode
*/
tcPhaseSpaceModePtr mode(unsigned int ix) const {
return modes_[ix];
}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DecayIntegrator2 & operator=(const DecayIntegrator2 &);
/**
* Parameters for the integration
*/
//@{
/**
* Number of iterations for th initialization.
*/
unsigned int nIter_;
/**
* Number of points for initialisation
*/
unsigned int nPoint_;
/**
* number of attempts to generate the decay
*/
unsigned int nTry_;
/**
* List of the decay modes
*/
mutable vector<PhaseSpaceModePtr> 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_;
protected:
/**
* Exception for this class and those inheriting from it
*/
class DecayIntegrator2Error: public Exception {};
};
/**
* Output information on the DecayIntegrator for debugging purposes
*/
ostream & operator<<(ostream &, const DecayIntegrator2 &);
}
#endif /* Herwig_DecayIntegrator2_H */
diff --git a/Decay/General/GeneralTwoBodyDecayer.cc b/Decay/General/GeneralTwoBodyDecayer2.cc
copy from Decay/General/GeneralTwoBodyDecayer.cc
copy to Decay/General/GeneralTwoBodyDecayer2.cc
--- a/Decay/General/GeneralTwoBodyDecayer.cc
+++ b/Decay/General/GeneralTwoBodyDecayer2.cc
@@ -1,822 +1,819 @@
// -*- C++ -*-
//
-// GeneralTwoBodyDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
+// GeneralTwoBodyDecayer2.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 GeneralTwoBodyDecayer class.
+// functions of the GeneralTwoBodyDecayer2 class.
//
-#include "GeneralTwoBodyDecayer.h"
+#include "GeneralTwoBodyDecayer2.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Utilities/Exception.h"
#include "Herwig/Shower/RealEmissionProcess.h"
using namespace Herwig;
-ParticleVector GeneralTwoBodyDecayer::decay(const Particle & parent,
+ParticleVector GeneralTwoBodyDecayer2::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;
int imode=modeNumber(cc,parent.dataPtr(),children);
// generate the kinematics
ParticleVector decay=generate(generateIntermediates(),cc,imode,parent);
// make the colour connections
colourConnections(parent, decay);
// return the answer
return decay;
}
-void GeneralTwoBodyDecayer::doinit() {
- PerturbativeDecayer::doinit();
+void GeneralTwoBodyDecayer2::doinit() {
+ PerturbativeDecayer2::doinit();
assert( incoming_ && outgoing_.size()==2);
//create phase space mode
- tPDVector extpart(3);
- extpart[0] = incoming_;
- extpart[1] = outgoing_[0];
- extpart[2] = outgoing_[1];
- addMode(new_ptr(DecayPhaseSpaceMode(extpart, this)), maxWeight_, vector<double>());
+ addMode(new_ptr(PhaseSpaceMode(incoming_,{outgoing_[0],outgoing_[1]},
+ maxWeight_)));
}
-int GeneralTwoBodyDecayer::modeNumber(bool & cc, tcPDPtr parent,
+int GeneralTwoBodyDecayer2::modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const {
long parentID = parent->id();
long id1 = children[0]->id();
long id2 = children[1]->id();
cc = false;
long out1 = outgoing_[0]->id();
long out2 = outgoing_[1]->id();
if( parentID == incoming_->id() &&
((id1 == out1 && id2 == out2) ||
(id1 == out2 && id2 == out1)) ) {
return 0;
}
else if(incoming_->CC() && parentID == incoming_->CC()->id()) {
cc = true;
if( outgoing_[0]->CC()) out1 = outgoing_[0]->CC()->id();
if( outgoing_[1]->CC()) out2 = outgoing_[1]->CC()->id();
if((id1 == out1 && id2 == out2) ||
(id1 == out2 && id2 == out1)) return 0;
}
return -1;
}
-void GeneralTwoBodyDecayer::
+void GeneralTwoBodyDecayer2::
colourConnections(const Particle & parent,
const ParticleVector & out) const {
PDT::Colour incColour(parent.data().iColour());
PDT::Colour outaColour(out[0]->data().iColour());
PDT::Colour outbColour(out[1]->data().iColour());
//incoming colour singlet
if(incColour == PDT::Colour0) {
// colour triplet-colourantitriplet
if((outaColour == PDT::Colour3 && outbColour == PDT::Colour3bar) ||
(outaColour == PDT::Colour3bar && outbColour == PDT::Colour3)) {
bool ac(out[0]->id() < 0);
out[0]->colourNeighbour(out[1],!ac);
}
//colour octet
else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour8) {
out[0]->colourNeighbour(out[1]);
out[0]->antiColourNeighbour(out[1]);
}
// colour singlets
else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour0) {
}
// unknown
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour singlet in "
- << "GeneralTwoBodyDecayer::colourConnections "
+ << "GeneralTwoBodyDecayer2::colourConnections "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
//incoming colour triplet
else if(incColour == PDT::Colour3) {
// colour triplet + singlet
if(outaColour == PDT::Colour3 && outbColour == PDT::Colour0) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
}
//opposite order
else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3) {
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
}
// octet + triplet
else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[1]->antiColourNeighbour(out[0]);
}
//opposite order
else if(outaColour == PDT::Colour3 && outbColour == PDT::Colour8) {
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[0]->antiColourNeighbour(out[1]);
}
else if(outaColour == PDT::Colour3bar && outaColour == PDT::Colour3bar) {
tColinePtr col[2] = {ColourLine::create(out[0],true),
ColourLine::create(out[1],true)};
parent.colourLine()->setSinkNeighbours(col[0],col[1]);
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour triplet in "
- << "GeneralTwoBodyDecayer::colourConnections() "
+ << "GeneralTwoBodyDecayer2::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
// incoming colour anti triplet
else if(incColour == PDT::Colour3bar) {
// colour antitriplet +singlet
if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour0) {
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
//opposite order
else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3bar) {
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
//octet + antitriplet
else if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour8) {
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
out[0]->colourNeighbour(out[1]);
}
//opposite order
else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3bar) {
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
out[1]->colourNeighbour(out[0]);
}
else if(outaColour == PDT::Colour3 && outbColour == PDT::Colour3) {
tColinePtr col[2] = {ColourLine::create(out[0]),
ColourLine::create(out[1])};
parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour antitriplet "
- << "in GeneralTwoBodyDecayer::colourConnections() "
+ << "in GeneralTwoBodyDecayer2::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
//incoming colour octet
else if(incColour == PDT::Colour8) {
// triplet-antitriplet
if(outaColour == PDT::Colour3&&outbColour == PDT::Colour3bar) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
// opposite order
else if(outbColour == PDT::Colour3&&outaColour == PDT::Colour3bar) {
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
}
// neutral octet
else if(outaColour == PDT::Colour0&&outbColour == PDT::Colour8) {
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
else if(outbColour == PDT::Colour0&&outaColour == PDT::Colour8) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour octet "
- << "in GeneralTwoBodyDecayer::colourConnections() "
+ << "in GeneralTwoBodyDecayer2::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
else if(incColour == PDT::Colour6) {
if(outaColour == PDT::Colour3 && outbColour == PDT::Colour3) {
tPPtr tempParent = const_ptr_cast<tPPtr>(&parent);
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(tempParent->colourInfo());
tColinePtr line1 = const_ptr_cast<tColinePtr>(parentColour->colourLines()[0]);
line1->addColoured(dynamic_ptr_cast<tPPtr>(out[0]));
tColinePtr line2 = const_ptr_cast<tColinePtr>(parentColour->colourLines()[1]);
line2->addColoured(dynamic_ptr_cast<tPPtr>(out[1]));
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour sextet "
- << "in GeneralTwoBodyDecayer::colourConnections() "
+ << "in GeneralTwoBodyDecayer2::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
else if(incColour == PDT::Colour6bar) {
if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour3bar) {
tPPtr tempParent = const_ptr_cast<tPPtr>(&parent);
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(tempParent->colourInfo());
tColinePtr line1 = const_ptr_cast<tColinePtr>(parentColour->antiColourLines()[0]);
line1->addAntiColoured(dynamic_ptr_cast<tPPtr>(out[0]));
tColinePtr line2 = const_ptr_cast<tColinePtr>(parentColour->antiColourLines()[1]);
line2->addAntiColoured(dynamic_ptr_cast<tPPtr>(out[1]));
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour anti-sextet "
- << "in GeneralTwoBodyDecayer::colourConnections() "
+ << "in GeneralTwoBodyDecayer2::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
else
throw Exception() << "Unknown incoming colour in "
- << "GeneralTwoBodyDecayer::colourConnections() "
+ << "GeneralTwoBodyDecayer2::colourConnections() "
<< incColour
<< Exception::runerror;
}
-bool GeneralTwoBodyDecayer::twoBodyMEcode(const DecayMode & dm, int & mecode,
+bool GeneralTwoBodyDecayer2::twoBodyMEcode(const DecayMode & dm, int & mecode,
double & coupling) const {
assert(dm.parent()->id() == incoming_->id());
ParticleMSet::const_iterator pit = dm.products().begin();
long id1 = (*pit)->id();
++pit;
long id2 = (*pit)->id();
long id1t(outgoing_[0]->id()), id2t(outgoing_[1]->id());
mecode = -1;
coupling = 1.;
if( id1 == id1t && id2 == id2t ) {
return true;
}
else if( id1 == id2t && id2 == id1t ) {
return false;
}
else
assert(false);
return false;
}
-void GeneralTwoBodyDecayer::persistentOutput(PersistentOStream & os) const {
+void GeneralTwoBodyDecayer2::persistentOutput(PersistentOStream & os) const {
os << incoming_ << outgoing_ << maxWeight_;
}
-void GeneralTwoBodyDecayer::persistentInput(PersistentIStream & is, int) {
+void GeneralTwoBodyDecayer2::persistentInput(PersistentIStream & is, int) {
is >> incoming_ >> outgoing_ >> maxWeight_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeAbstractClass<GeneralTwoBodyDecayer,PerturbativeDecayer>
-describeHerwigGeneralTwoBodyDecayer("Herwig::GeneralTwoBodyDecayer", "Herwig.so");
+DescribeAbstractClass<GeneralTwoBodyDecayer2,PerturbativeDecayer2>
+describeHerwigGeneralTwoBodyDecayer2("Herwig::GeneralTwoBodyDecayer2", "Herwig.so");
-void GeneralTwoBodyDecayer::Init() {
+void GeneralTwoBodyDecayer2::Init() {
- static ClassDocumentation<GeneralTwoBodyDecayer> documentation
+ static ClassDocumentation<GeneralTwoBodyDecayer2> documentation
("This class is designed to be a base class for all 2 body decays"
"in a general model");
}
-double GeneralTwoBodyDecayer::brat(const DecayMode &, const Particle & p,
+double GeneralTwoBodyDecayer2::brat(const DecayMode &, const Particle & p,
double oldbrat) const {
ParticleVector children = p.children();
if( children.size() != 2 || !p.data().widthGenerator() )
return oldbrat;
// partial width for this mode
Energy scale = p.mass();
Energy pwidth =
partialWidth( make_pair(p.dataPtr(), scale),
make_pair(children[0]->dataPtr(), children[0]->mass()),
make_pair(children[1]->dataPtr(), children[1]->mass()) );
Energy width = p.data().widthGenerator()->width(p.data(), scale);
return pwidth/width;
}
-void GeneralTwoBodyDecayer::doinitrun() {
- PerturbativeDecayer::doinitrun();
+void GeneralTwoBodyDecayer2::doinitrun() {
+ PerturbativeDecayer2::doinitrun();
for(unsigned int ix=0;ix<numberModes();++ix) {
- double fact = pow(1.5,int(mode(ix)->externalParticles(0)->iSpin())-1);
- mode(ix)->setMaxWeight(fact*mode(ix)->maxWeight());
+ double fact = pow(1.5,int(mode(ix)->incoming().first->iSpin())-1);
+ mode(ix)->maxWeight(fact*mode(ix)->maxWeight());
}
}
-double GeneralTwoBodyDecayer::colourFactor(tcPDPtr in, tcPDPtr out1,
+double GeneralTwoBodyDecayer2::colourFactor(tcPDPtr in, tcPDPtr out1,
tcPDPtr out2) const {
// identical particle symmetry factor
double output = out1->id()==out2->id() ? 0.5 : 1.;
// colour neutral incoming particle
if(in->iColour()==PDT::Colour0) {
// both colour neutral
if(out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour0)
output *= 1.;
// colour triplet/ antitriplet
else if((out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3bar) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3 ) ) {
output *= 3.;
}
// colour octet colour octet
else if(out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour8 ) {
output *= 8.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour neutral particle in "
- << "GeneralTwoBodyDecayer::colourFactor() for "
+ << "GeneralTwoBodyDecayer2::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
// triplet
else if(in->iColour()==PDT::Colour3) {
// colour triplet + neutral
if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour3) ||
(out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour0) ) {
output *= 1.;
}
// colour triplet + octet
else if((out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour3) ||
(out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour8) ) {
output *= 4./3.;
}
// colour anti triplet anti triplet
else if(out1->iColour()==PDT::Colour3bar &&
out2->iColour()==PDT::Colour3bar) {
output *= 2.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour triplet particle in "
- << "GeneralTwoBodyDecayer::colourFactor() for "
+ << "GeneralTwoBodyDecayer2::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
// anti triplet
else if(in->iColour()==PDT::Colour3bar) {
// colour anti triplet + neutral
if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour3bar ) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour0 ) ) {
output *= 1.;
}
// colour anti triplet + octet
else if((out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour3bar ) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour8 ) ) {
output *= 4./3.;
}
// colour triplet triplet
else if(out1->iColour()==PDT::Colour3 &&
out2->iColour()==PDT::Colour3) {
output *= 2.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour anti triplet particle in "
- << "GeneralTwoBodyDecayer::colourFactor() for "
+ << "GeneralTwoBodyDecayer2::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else if(in->iColour()==PDT::Colour8) {
// colour octet + neutral
if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour8 ) ||
(out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour0 ) ) {
output *= 1.;
}
// colour triplet/antitriplet
else if((out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3bar) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3 ) ) {
output *= 0.5;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour octet particle in "
- << "GeneralTwoBodyDecayer::colourFactor() for "
+ << "GeneralTwoBodyDecayer2::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else if(in->iColour()==PDT::Colour6) {
// colour sextet -> triplet triplet
if( out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3 ) {
output *= 1.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour sextet particle in "
- << "GeneralTwoBodyDecayer::colourFactor() for "
+ << "GeneralTwoBodyDecayer2::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else if(in->iColour()==PDT::Colour6bar) {
// colour sextet -> triplet triplet
if( out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3bar ) {
output *= 1.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour anti-sextet particle in "
- << "GeneralTwoBodyDecayer::colourFactor() for "
+ << "GeneralTwoBodyDecayer2::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else
throw Exception() << "Unknown colour "
<< in->iColour() << " for the decaying particle in "
- << "GeneralTwoBodyDecayer::colourFactor() for "
+ << "GeneralTwoBodyDecayer2::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
return output;
}
-Energy GeneralTwoBodyDecayer::partialWidth(PMPair inpart, PMPair outa,
+Energy GeneralTwoBodyDecayer2::partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const {
// select the number of the mode
tPDVector children;
children.push_back(const_ptr_cast<PDPtr>(outa.first));
children.push_back(const_ptr_cast<PDPtr>(outb.first));
bool cc;
int nmode=modeNumber(cc,inpart.first,children);
- tcPDPtr newchild[2] = {mode(nmode)->externalParticles(1),
- mode(nmode)->externalParticles(2)};
+ tcPDPtr newchild[2] = {mode(nmode)->outgoing()[0],
+ mode(nmode)->outgoing()[1]};
// make the particles
Lorentz5Momentum pparent = Lorentz5Momentum(inpart.second);
PPtr parent = inpart.first->produceParticle(pparent);
Lorentz5Momentum pout[2];
double ctheta,phi;
Kinematics::generateAngles(ctheta,phi);
Kinematics::twoBodyDecay(pparent, outa.second, outb.second,
ctheta, phi,pout[0],pout[1]);
if( ( !cc && outa.first!=newchild[0]) ||
( cc && !(( outa.first->CC() && outa.first->CC() == newchild[0])||
( !outa.first->CC() && outa.first == newchild[0]) )))
swap(pout[0],pout[1]);
ParticleVector decay;
decay.push_back(newchild[0]->produceParticle(pout[0]));
decay.push_back(newchild[1]->produceParticle(pout[1]));
double me = me2(-1,*parent,decay,Initialize);
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,
outa.second, outb.second);
return me/(8.*Constants::pi)*pcm;
}
-void GeneralTwoBodyDecayer::decayInfo(PDPtr incoming, PDPair outgoing) {
+void GeneralTwoBodyDecayer2::decayInfo(PDPtr incoming, PDPair outgoing) {
incoming_=incoming;
outgoing_.clear();
outgoing_.push_back(outgoing.first );
outgoing_.push_back(outgoing.second);
}
-double GeneralTwoBodyDecayer::matrixElementRatio(const Particle & inpart,
+double GeneralTwoBodyDecayer2::matrixElementRatio(const Particle & inpart,
const ParticleVector & decay2,
const ParticleVector & decay3,
MEOption meopt,
ShowerInteraction inter) {
// calculate R/B
double B = me2 (0, inpart, decay2, meopt);
double R = threeBodyME(0, inpart, decay3, inter, meopt);
return R/B;
}
-const vector<DVector> & GeneralTwoBodyDecayer::getColourFactors(const Particle & inpart,
+const vector<DVector> & GeneralTwoBodyDecayer2::getColourFactors(const Particle & inpart,
const ParticleVector & decay,
unsigned int & nflow) {
// calculate the colour factors for the three-body decay
vector<int> sing,trip,atrip,oct,sex,asex;
for(unsigned int it=0;it<decay.size();++it) {
if (decay[it]->dataPtr()->iColour() == PDT::Colour0 ) sing. push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour3 ) trip. push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour3bar ) atrip.push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour8 ) oct. push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour6 ) sex. push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour6bar ) asex. push_back(it);
}
// identical particle symmetry factor
double symFactor=1.;
if (( sing.size()==2 && decay[ sing[0]]->id()==decay[ sing[1]]->id()) ||
( trip.size()==2 && decay[ trip[0]]->id()==decay[ trip[1]]->id()) ||
(atrip.size()==2 && decay[atrip[0]]->id()==decay[atrip[1]]->id()) ||
( oct.size()==2 && decay[ oct[0]]->id()==decay[ oct[1]]->id()) ||
( sex.size()==2 && decay[ sex[0]]->id()==decay[ sex[1]]->id()) ||
( asex.size()==2 && decay[ asex[0]]->id()==decay[ asex[1]]->id()))
symFactor /= 2.;
else if (oct.size()==3 &&
decay[oct[0]]->id()==decay[oct[1]]->id() &&
decay[oct[0]]->id()==decay[oct[2]]->id())
symFactor /= 6.;
colour_ = vector<DVector>(1,DVector(1,symFactor*1.));
// decaying colour singlet
if(inpart.dataPtr()->iColour() == PDT::Colour0) {
if(trip.size()==1 && atrip.size()==1 && oct.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*4.));
}
else if(trip.size()==1 && atrip.size()==1 && sing.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*3.));
}
else if (oct.size()==3){
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*24.));
}
else if(sing.size()==3) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor));
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour scalar particle in "
- << "GeneralTwoBodyDecayer::getColourFactors() for "
+ << "GeneralTwoBodyDecayer2::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// decaying colour triplet
else if(inpart.dataPtr()->iColour() == PDT::Colour3) {
if(trip.size()==1 && sing.size()==1 && oct.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*4./3.));
}
else if(trip.size()==1 && sing.size()==2) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor));
}
else if(trip.size()==1 && oct.size()==2) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = symFactor*16./9.; colour_[0][1] = -symFactor*2./9.;
colour_[1][0] = -symFactor*2./9.; colour_[1][1] = symFactor*16./9.;
}
else if(atrip.size()==2 && sing.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*2.));
}
else if(atrip.size()==2 && oct.size()==1) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = 8./3.*symFactor; colour_[0][1] =-4./3.*symFactor;
colour_[1][0] = -4./3.*symFactor; colour_[1][1] = 8./3.*symFactor;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour triplet particle in "
- << "GeneralTwoBodyDecayer::getColourFactors() for "
+ << "GeneralTwoBodyDecayer2::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// decaying colour anti-triplet
else if(inpart.dataPtr()->iColour() == PDT::Colour3bar) {
if(atrip.size()==1 && sing.size()==1 && oct.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*4./3.));
}
else if(atrip.size()==1 && sing.size()==2) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor));
}
else if(atrip.size()==1 && oct.size()==2){
nflow = 2;
colour_.clear();
colour_ .resize(2,DVector(2,0.));
colour_[0][0] = symFactor*16./9.; colour_[0][1] = -symFactor*2./9.;
colour_[1][0] = -symFactor*2./9.; colour_[1][1] = symFactor*16./9.;
}
else if(trip.size()==2 && sing.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*2.));
}
else if(trip.size()==2 && oct.size()==1) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = 8./3.*symFactor; colour_[0][1] =-4./3.*symFactor;
colour_[1][0] = -4./3.*symFactor; colour_[1][1] = 8./3.*symFactor;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour anti-triplet particle in "
- << "GeneralTwoBodyDecayer::getColourFactors() for "
+ << "GeneralTwoBodyDecayer2::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// decaying colour octet
else if(inpart.dataPtr()->iColour() == PDT::Colour8) {
if(oct.size()==1 && trip.size()==1 && atrip.size()==1) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = symFactor*2./3. ; colour_[0][1] = -symFactor*1./12.;
colour_[1][0] = -symFactor*1./12.; colour_[1][1] = symFactor*2./3. ;
}
else if (sing.size()==1 && trip.size()==1 && atrip.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*0.5));
}
else if (oct.size()==2 && sing.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*3.));
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for a decaying colour octet particle in "
- << "GeneralTwoBodyDecayer::getColourFactors() for "
+ << "GeneralTwoBodyDecayer2::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// Sextet
else if(inpart.dataPtr()->iColour() == PDT::Colour6) {
if(trip.size()==2 && sing.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor));
}
else if(trip.size()==2 && oct.size()==1) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = 4./3.*symFactor; colour_[0][1] = 1./3.*symFactor;
colour_[1][0] = 1./3.*symFactor; colour_[1][1] = 4./3.*symFactor;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for a decaying colour sextet particle in "
- << "GeneralTwoBodyDecayer::getColourFactors() for "
+ << "GeneralTwoBodyDecayer2::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// anti Sextet
else if(inpart.dataPtr()->iColour() == PDT::Colour6bar) {
if(atrip.size()==2 && sing.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor));
}
else if(atrip.size()==2 && oct.size()==1) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = 4./3.*symFactor; colour_[0][1] = 1./3.*symFactor;
colour_[1][0] = 1./3.*symFactor; colour_[1][1] = 4./3.*symFactor;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for a decaying colour anti-sextet particle in "
- << "GeneralTwoBodyDecayer::getColourFactors() for "
+ << "GeneralTwoBodyDecayer2::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
else
throw Exception() << "Unknown colour for the decaying particle in "
- << "GeneralTwoBodyDecayer::getColourFactors() for "
+ << "GeneralTwoBodyDecayer2::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
return colour_;
}
-const GeneralTwoBodyDecayer::CFlow &
-GeneralTwoBodyDecayer::colourFlows(const Particle & inpart,
+const GeneralTwoBodyDecayer2::CFlow &
+GeneralTwoBodyDecayer2::colourFlows(const Particle & inpart,
const ParticleVector & decay) {
// static initialization of commonly used colour structures
static const CFlow init = CFlow(3, CFlowPairVec(1, make_pair(0, 1.)));
static CFlow tripflow = init;
static CFlow atripflow = init;
static CFlow octflow = init;
static CFlow sexflow = init;
static CFlow epsflow = init;
static const CFlow fpflow = CFlow(4, CFlowPairVec(1, make_pair(0, 1.)));
static bool initialized = false;
if (! initialized) {
tripflow[2].resize(2, make_pair(0,1.));
tripflow[2][0] = make_pair(0, 1.);
tripflow[2][1] = make_pair(1,-1.);
tripflow[1][0] = make_pair(1, 1.);
atripflow[1].resize(2, make_pair(0,1.));
atripflow[1][0] = make_pair(0, 1.);
atripflow[1][1] = make_pair(1,-1.);
atripflow[2][0] = make_pair(1, 1.);
octflow[0].resize(2, make_pair(0,1.));
octflow[0][0] = make_pair(0,-1.);
octflow[0][1] = make_pair(1, 1.);
octflow[2][0] = make_pair(1, 1.);
sexflow[0].resize(2, make_pair(0,1.));
sexflow[0][1] = make_pair(1,1.);
sexflow[2][0] = make_pair(1,1.);
epsflow[0].resize(2, make_pair(0,-1.));
epsflow[0][0] = make_pair(0,-1.);
epsflow[0][1] = make_pair(1,-1.);
epsflow[2][0] = make_pair(1,1.);
initialized = true;
}
// main function body
int sing=0,trip=0,atrip=0,oct=0,sex=0,asex=0;
for (size_t it=0; it<decay.size(); ++it) {
switch ( decay[it]->dataPtr()->iColour() ) {
case PDT::Colour0: ++sing; break;
case PDT::Colour3: ++trip; break;
case PDT::Colour3bar: ++atrip; break;
case PDT::Colour8: ++oct; break;
case PDT::Colour6: ++sex; break;
case PDT::Colour6bar: ++asex; break;
/// @todo: handle these better
case PDT::ColourUndefined: break;
case PDT::Coloured: break;
}
}
const CFlow * retval = 0;
// decaying colour triplet
if(inpart.dataPtr()->iColour() == PDT::Colour3 &&
trip==1 && oct==2) {
retval = &tripflow;
}
// decaying colour anti-triplet
else if(inpart.dataPtr()->iColour() == PDT::Colour3bar &&
atrip==1 && oct==2){
retval = &atripflow;
}
// decaying colour octet
else if(inpart.dataPtr()->iColour() == PDT::Colour8 &&
oct==1 && trip==1 && atrip==1) {
retval = &octflow;
}
// decaying colour sextet
else if(inpart.dataPtr()->iColour() ==PDT::Colour6 &&
oct==1 && trip==2) {
retval = &sexflow;
}
// decaying anti colour sextet
else if(inpart.dataPtr()->iColour() ==PDT::Colour6bar &&
oct==1 && atrip==2) {
retval = &sexflow;
}
// decaying colour triplet (eps)
else if(inpart.dataPtr()->iColour() == PDT::Colour3 &&
atrip==2 && oct==1) {
retval = &epsflow;
}
// decaying colour anti-triplet (eps)
else if(inpart.dataPtr()->iColour() == PDT::Colour3bar &&
trip==2 && oct==1){
retval = &epsflow;
}
else {
retval = &fpflow;
}
return *retval;
}
-double GeneralTwoBodyDecayer::threeBodyME(const int , const Particle &,
+double GeneralTwoBodyDecayer2::threeBodyME(const int , const Particle &,
const ParticleVector &,
ShowerInteraction, MEOption) {
- throw Exception() << "Base class PerturbativeDecayer::threeBodyME() "
+ throw Exception() << "Base class GeneralTwoBodyDecayer2::threeBodyME() "
<< "called, should have an implementation in the inheriting class"
<< Exception::runerror;
return 0.;
}
diff --git a/Decay/General/GeneralTwoBodyDecayer.fh b/Decay/General/GeneralTwoBodyDecayer2.fh
copy from Decay/General/GeneralTwoBodyDecayer.fh
copy to Decay/General/GeneralTwoBodyDecayer2.fh
--- a/Decay/General/GeneralTwoBodyDecayer.fh
+++ b/Decay/General/GeneralTwoBodyDecayer2.fh
@@ -1,17 +1,17 @@
// -*- C++ -*-
//
-// This is the forward declaration of the GeneralTwoBodyDecayer class.
+// This is the forward declaration of the GeneralTwoBodyDecayer2 class.
//
-#ifndef ThePEG_GeneralTwoBodyDecayer_FH
-#define ThePEG_GeneralTwoBodyDecayer_FH
+#ifndef ThePEG_GeneralTwoBodyDecayer2_FH
+#define ThePEG_GeneralTwoBodyDecayer2_FH
#include "ThePEG/Config/Pointers.h"
namespace Herwig {
using namespace ThePEG;
-class GeneralTwoBodyDecayer;
-ThePEG_DECLARE_POINTERS(GeneralTwoBodyDecayer,GeneralTwoBodyDecayerPtr);
+class GeneralTwoBodyDecayer2;
+ThePEG_DECLARE_POINTERS(GeneralTwoBodyDecayer2,GeneralTwoBodyDecayer2Ptr);
}
#endif
diff --git a/Decay/General/GeneralTwoBodyDecayer.h b/Decay/General/GeneralTwoBodyDecayer2.h
copy from Decay/General/GeneralTwoBodyDecayer.h
copy to Decay/General/GeneralTwoBodyDecayer2.h
--- a/Decay/General/GeneralTwoBodyDecayer.h
+++ b/Decay/General/GeneralTwoBodyDecayer2.h
@@ -1,306 +1,306 @@
// -*- C++ -*-
//
// GeneralTwoBodyDecayer.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_GeneralTwoBodyDecayer_H
-#define HERWIG_GeneralTwoBodyDecayer_H
+#ifndef HERWIG_GeneralTwoBodyDecayer2_H
+#define HERWIG_GeneralTwoBodyDecayer2_H
//
-// This is the declaration of the GeneralTwoBodyDecayer class.
+// This is the declaration of the GeneralTwoBodyDecayer2 class.
//
-#include "Herwig/Decay/PerturbativeDecayer.h"
-#include "Herwig/Decay/DecayPhaseSpaceMode.h"
+#include "Herwig/Decay/PerturbativeDecayer2.h"
+#include "Herwig/Decay/PhaseSpaceMode.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
-#include "GeneralTwoBodyDecayer.fh"
+#include "GeneralTwoBodyDecayer2.fh"
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
/** \ingroup Decay
- * The GeneralTwoBodyDecayer class is designed to be the base class
+ * The GeneralTwoBodyDecayer2 class is designed to be the base class
* for 2 body decays for some general model. It inherits from
- * PerturbativeDecayer and implements the modeNumber() virtual function
+ * PerturbativeDecayer2 and implements the modeNumber() virtual function
* that is the same for all of the decays. A decayer for
* a specific spin configuration should inherit from this and implement
* the me2() and partialWidth() member functions. The colourConnections()
* member should be called from inside me2() in the inheriting decayer
* to set up the colour lines.
*
- * @see \ref GeneralTwoBodyDecayerInterfaces "The interfaces"
- * defined for GeneralTwoBodyDecayer.
- * @see PerturbativeDecayer
+ * @see \ref GeneralTwoBodyDecayer2Interfaces "The interfaces"
+ * defined for GeneralTwoBodyDecayer2.
+ * @see PerturbativeDecayer2
*/
-class GeneralTwoBodyDecayer: public PerturbativeDecayer {
+class GeneralTwoBodyDecayer2: public PerturbativeDecayer2 {
public:
/** A ParticleData ptr and (possible) mass pair.*/
typedef pair<tcPDPtr, Energy> PMPair;
public:
/**
* The default constructor.
*/
- GeneralTwoBodyDecayer() : maxWeight_(1.), colour_(1,DVector(1,1.))
+ GeneralTwoBodyDecayer2() : maxWeight_(1.), colour_(1,DVector(1,1.))
{}
- /** @name Virtual functions required by the Decayer class. */
+ /** @name Virtual functions required by the Decayer2 class. */
//@{
/**
* 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;
/**
* Return the matrix element squared for a given mode and phase-space channel
* @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
* @param decay The particles produced in the decay.
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
virtual double me2(const int , const Particle & part,
const ParticleVector & decay, MEOption meopt) const = 0;
/**
* Function to return partial Width
* @param inpart The decaying particle.
* @param outa One of the decay products.
* @param outb The other decay product.
*/
virtual Energy partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const;
/**
* Specify the \f$1\to2\f$ matrix element to be used in the running width
* calculation.
* @param dm The DecayMode
* @param mecode The code for the matrix element as described
* in the GenericWidthGenerator class.
* @param coupling The coupling for the matrix element.
* @return True if the the order of the particles in the
* decayer is the same as the DecayMode tag.
*/
virtual bool twoBodyMEcode(const DecayMode & dm, int & mecode,
double & coupling) const;
/**
* An overidden member to calculate a branching ratio for a certain
* particle instance.
* @param dm The DecayMode of the particle
* @param p The particle object
* @param oldbrat The branching fraction given in the DecayMode object
*/
virtual double brat(const DecayMode & dm, const Particle & p,
double oldbrat) const;
//@}
/**
* Set the information on the decay
*/
virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr>,
map<ShowerInteraction,VertexBasePtr> &,
const vector<map<ShowerInteraction,VertexBasePtr> > &,
map<ShowerInteraction,VertexBasePtr>) =0;
protected:
/** @name Functions used by inheriting decayers. */
//@{
/**
* Set integration weight
* @param wgt Maximum integration weight
*/
void setWeight(double wgt) { maxWeight_ = wgt; }
/**
* Set colour connections
* @param parent Parent particle
* @param out Particle vector containing particles to
* connect colour lines
*/
void colourConnections(const Particle & parent,
const ParticleVector & out) const;
/**
* Compute the spin and colour factor
*/
double colourFactor(tcPDPtr in, tcPDPtr out1, tcPDPtr out2) const;
/**
* Calculate matrix element ratio R/B
*/
double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
const ParticleVector & decay3, MEOption meopt,
ShowerInteraction inter);
/**
* Set the information on the decay
*/
void decayInfo(PDPtr incoming, PDPair outgoing);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/**
* Member for the generation of additional hard radiation
*/
//@{
/**
* Return the matrix of colour factors
*/
typedef vector<pair<int,double > > CFlowPairVec;
typedef vector<CFlowPairVec> CFlow;
const vector<DVector> & getColourFactors(const Particle & inpart,
const ParticleVector & decay,
unsigned int & nflow);
const CFlow & colourFlows(const Particle & inpart,
const ParticleVector & decay);
/**
* Three-body matrix element including additional QCD radiation
*/
virtual double threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt);
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- GeneralTwoBodyDecayer & operator=(const GeneralTwoBodyDecayer &);
+ GeneralTwoBodyDecayer2 & operator=(const GeneralTwoBodyDecayer2 &);
private:
/**
* Store the incoming particle
*/
PDPtr incoming_;
/**
* Outgoing particles
*/
vector<PDPtr> outgoing_;
/**
* Maximum weight for integration
*/
double maxWeight_;
/**
* Store colour factors for ME calc.
*/
vector<DVector> colour_;
};
/**
* Write a map with ShowerInteraction as the key
*/
template<typename T, typename Cmp, typename A>
inline PersistentOStream & operator<<(PersistentOStream & os,
const map<ShowerInteraction,T,Cmp,A> & m) {
os << m.size();
if(m.find(ShowerInteraction::QCD)!=m.end()) {
os << 0 << m.at(ShowerInteraction::QCD);
}
if(m.find(ShowerInteraction::QED)!=m.end()) {
os << 1 << m.at(ShowerInteraction::QED);
}
return os;
}
/**
* Read a map with ShowerInteraction as the key
*/
template <typename T, typename Cmp, typename A>
inline PersistentIStream & operator>>(PersistentIStream & is, map<ShowerInteraction,T,Cmp,A> & m) {
m.clear();
long size;
int k;
is >> size;
while ( size-- && is ) {
is >> k;
if(k==0)
is >> m[ShowerInteraction::QCD];
else if(k==1)
is >> m[ShowerInteraction::QED];
else
assert(false);
}
return is;
}
}
-#endif /* HERWIG_GeneralTwoBodyDecayer_H */
+#endif /* HERWIG_GeneralTwoBodyDecayer2_H */
diff --git a/Decay/General/Makefile.am b/Decay/General/Makefile.am
--- a/Decay/General/Makefile.am
+++ b/Decay/General/Makefile.am
@@ -1,74 +1,75 @@
noinst_LTLIBRARIES = libHwGeneralDecay.la
libHwGeneralDecay_la_SOURCES = \
GeneralTwoBodyDecayer.h GeneralTwoBodyDecayer.fh GeneralTwoBodyDecayer.cc \
+GeneralTwoBodyDecayer2.h GeneralTwoBodyDecayer2.fh GeneralTwoBodyDecayer2.cc \
GeneralThreeBodyDecayer.h GeneralThreeBodyDecayer.fh \
GeneralThreeBodyDecayer.cc \
GeneralCurrentDecayer.h GeneralCurrentDecayer.fh GeneralCurrentDecayer.cc \
GeneralFourBodyDecayer.h GeneralFourBodyDecayer.fh \
GeneralFourBodyDecayer.cc \
StoFFFFDecayer.h StoFFFFDecayer.cc
# check the gauge invariance of corrections (for testing ONLY)
#libHwGeneralDecay_la_CPPFLAGS= $(AM_CPPFLAGS) -DGAUGE_CHECK
nodist_libHwGeneralDecay_la_SOURCES = \
GeneralDecayer__all.cc
BUILT_SOURCES = GeneralDecayer__all.cc
CLEANFILES = GeneralDecayer__all.cc
GeneralDecayer__all.cc : $(DIR_H_FILES) $(DIR_CC_FILES) Makefile
@echo "Concatenating .cc files into $@"
@$(top_srcdir)/cat_with_cpplines $(DIR_CC_FILES) > $@
EXTRA_DIST = $(ALL_H_FILES) $(ALL_CC_FILES)
DIR_H_FILES = $(addprefix $(srcdir)/,$(ALL_H_FILES))
ALL_H_FILES = \
FFSDecayer.h \
FFVDecayer.h \
SFFDecayer.h \
SSSDecayer.h \
SSVDecayer.h \
SVVDecayer.h \
TFFDecayer.h \
TSSDecayer.h \
TVVDecayer.h \
VFFDecayer.h \
VSSDecayer.h \
VVSDecayer.h \
VVVDecayer.h \
SRFDecayer.h \
FRSDecayer.h \
FRVDecayer.h \
FtoFFFDecayer.h \
StoSFFDecayer.h \
StoFFVDecayer.h \
VtoFFVDecayer.h \
FtoFVVDecayer.h \
FFVCurrentDecayer.h
DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES))
ALL_CC_FILES = \
FFSDecayer.cc \
FFVDecayer.cc \
SFFDecayer.cc \
SSSDecayer.cc \
SSVDecayer.cc \
SVVDecayer.cc \
TFFDecayer.cc \
TSSDecayer.cc \
TVVDecayer.cc \
VFFDecayer.cc \
VSSDecayer.cc \
VVSDecayer.cc \
VVVDecayer.cc \
SRFDecayer.cc \
FRSDecayer.cc \
FRVDecayer.cc \
FtoFFFDecayer.cc \
StoSFFDecayer.cc \
StoFFVDecayer.cc \
VtoFFVDecayer.cc \
FtoFVVDecayer.cc \
FFVCurrentDecayer.cc

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:37 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805858
Default Alt Text
(62 KB)

Event Timeline