Page MenuHomeHEPForge
Files F8723962

No OneTemporary
Actions

View Options

diff --git a/Decay/DecayIntegrator2.cc b/Decay/DecayIntegrator2.cc
--- a/Decay/DecayIntegrator2.cc
+++ b/Decay/DecayIntegrator2.cc
@@ -1,269 +1,277 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DecayIntegrator2 class.
//
#include "DecayIntegrator2.h"
#include "PhaseSpaceMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
void DecayIntegrator2::persistentOutput(PersistentOStream & os) const {
os << modes_ << nIter_ << nPoint_ << nTry_
<< photonGen_ << generateInter_ << ounit(eps_,GeV);
}
void DecayIntegrator2::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<DecayIntegrator2,HwDecayerBase>
describeHerwigDecayIntegrator2("Herwig::DecayIntegrator2", "Herwig.so");
void DecayIntegrator2::Init() {
static ClassDocumentation<DecayIntegrator2> documentation
("The DecayIntegrator class is a base decayer class "
"including a multi-channel integrator.");
static Parameter<DecayIntegrator2,unsigned int> interfaceIteration
("Iteration",
"Number of iterations for the initialization of the phase space",
&DecayIntegrator2::nIter_, 10, 0, 100,
false, false, true);
static Parameter<DecayIntegrator2,unsigned int> interfacePoints
("Points",
"number of phase space points to generate in the initialisation.",
&DecayIntegrator2::nPoint_, 10000, 1, 1000000000,
false, false, true);
static Parameter<DecayIntegrator2,unsigned int> interfaceNtry
("Ntry",
"Number of attempts to generate the decay",
&DecayIntegrator2::nTry_, 500, 0, 100000,
false, false, true);
static Reference<DecayIntegrator2,DecayRadiationGenerator> interfacePhotonGenerator
("PhotonGenerator",
"Object responsible for generating photons in the decay.",
&DecayIntegrator2::photonGen_, false, false, true, true, false);
static Switch<DecayIntegrator2,bool> interfaceGenerateIntermediates
("GenerateIntermediates",
"Whether or not to include intermediate particles in the output",
&DecayIntegrator2::generateInter_, false, false, false);
static SwitchOption interfaceGenerateIntermediatesNoIntermediates
(interfaceGenerateIntermediates,
"No",
"Don't include the intermediates",
false);
static SwitchOption interfaceGenerateIntermediatesIncludeIntermediates
(interfaceGenerateIntermediates,
"Yes",
"include the intermediates",
true);
}
double DecayIntegrator2::oneLoopVirtualME(unsigned int ,
const Particle &,
const ParticleVector &) {
throw DecayIntegrator2Error()
<< "DecayIntegrator2::oneLoopVirtualME() called. This should"
<< " have been overidden in an inheriting class if it is used"
<< Exception::runerror;
}
InvEnergy2 DecayIntegrator2::realEmissionME(unsigned int,
const Particle &,
ParticleVector &,
unsigned int,
double, double,
const LorentzRotation &,
const LorentzRotation &) {
throw DecayIntegrator2Error()
<< "DecayIntegrator2::realEmmisionME() called. This should"
<< " have been overidden in an inheriting class if it is used"
<< Exception::runerror;
}
ParticleVector DecayIntegrator2::decay(const Particle & parent,
const tPDVector & children) const {
// return empty vector if products heavier than parent
Energy mout(ZERO);
for(tPDPtr pd : children) mout += pd->massMin();
if(mout>parent.mass()) return ParticleVector();
// generate the decay
bool cc;
iMode_ = modeNumber(cc,parent.dataPtr(),children);
return modes_[iMode_]->generateDecay(parent,this,generateInter_,cc);
}
void DecayIntegrator2::doinitrun() {
HwDecayerBase::doinitrun();
if ( initialize() && Debug::level > 1 )
CurrentGenerator::current().log() << "Start of the initialisation for "
<< name() << "\n";
for(unsigned int ix=0;ix<modes_.size();++ix) {
if(!modes_[ix]) continue;
modes_[ix]->initrun();
iMode_=ix;
modes_[ix]->initializePhaseSpace(initialize(),this);
}
}
// output the information for the database
void DecayIntegrator2::dataBaseOutput(ofstream & output,bool header) const {
// header for MySQL
if(header) output << "update decayers set parameters=\"";
output << "newdef " << name() << ":Iteration " << nIter_ << "\n";
output << "newdef " << name() << ":Ntry " << nTry_ << "\n";
output << "newdef " << name() << ":Points " << nPoint_ << "\n";
//if(_photongen){;}
output << "newdef " << name() << ":GenerateIntermediates " << generateInter_ << " \n";
// footer for MySQL
if(header) {
output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";\n";
}
}
// set the code for the partial width
void DecayIntegrator2::setPartialWidth(const DecayMode & dm, int imode) {
int ifound = findMode(dm);
if(ifound>=0) modes_[ifound]->setPartialWidth(imode);
}
WidthCalculatorBasePtr
DecayIntegrator2::threeBodyMEIntegrator(const DecayMode &) const {
return WidthCalculatorBasePtr();
}
int DecayIntegrator2::findMode(const DecayMode & dm) {
int imode(-1);
vector<int> extid;
bool found(false);
int id;
unsigned int ix(0),iy,N,iz,tmax,nmatched;
if(modes_.size()==0) return -1;
do {
if(!modes_[ix]) {
++ix;
continue;
}
tcPDPtr in = modes_[ix]->incoming().first;
tcPDPtr cc = modes_[ix]->incoming().first->CC();
tmax=1;if(!cc){++tmax;}
for(iz=0;iz<tmax;++iz) {
extid.clear();
// check the parent
if(dm.parent()!=in && dm.parent()!=cc) continue;
if(dm.parent()->id()==in->id()&&iz==0) {
for(iy=0,N=modes_[ix]->numberOfParticles();iy<N;++iy) {
extid.push_back(modes_[ix]->outgoing()[iy]->id());
}
}
else if(dm.parent()->id()==in->id()&&iz==1) {
for(iy=0,N=modes_[ix]->numberOfParticles();iy<N;++iy) {
tcPDPtr cc2=modes_[ix]->outgoing()[iy]->CC();
extid.push_back( cc2 ? cc2->id() : modes_[ix]->outgoing()[iy]->id());
}
}
else if(cc&&dm.parent()->id()==cc->id()) {
for(iy=0,N=modes_[ix]->numberOfParticles();iy<N;++iy) {
tcPDPtr cc2 = modes_[ix]->outgoing()[iy]->CC();
extid.push_back( cc2 ? cc2->id() : modes_[ix]->outgoing()[iy]->id());
}
}
// if the parents match
if(!extid.empty()) {
vector<bool> matched(extid.size(),false);
bool done;
nmatched=0;
ParticleMSet::const_iterator pit = dm.products().begin();
do {
id=(**pit).id();
done=false;
iy=0;
do {
if(id==extid[iy]&&!matched[iy]) {
matched[iy]=true;
++nmatched;
done=true;
}
++iy;
}
while(iy<extid.size()&&!done);
++pit;
}
while(pit!=dm.products().end());
if(nmatched==extid.size()) {
imode=ix;
found=true;
}
}
}
++ix;
}
while(!found&&ix<modes_.size());
return imode;
}
// the matrix element to be integrated for the me
double DecayIntegrator2::threeBodyMatrixElement(const int,const Energy2,
const Energy2,
const Energy2,const Energy2,
const Energy, const Energy,
const Energy) const {
throw DecayIntegratorError()
<< "Calling the virtual DecayIntegrator2::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 DecayIntegrator2::threeBodydGammads(const int, const Energy2,
const Energy2,
const Energy, const Energy,
const Energy) const {
throw DecayIntegratorError()
<< "Calling the virtual DecayIntegrator2::threeBodydGammads()"
<<"method. This must be overwritten in the classes "
<< "inheriting from DecayIntegrator where it is needed"
<< Exception::runerror;
}
+// generate the momenta for the decay
+ParticleVector DecayIntegrator2::generate(bool inter,bool cc,
+ const unsigned int & imode,
+ const Particle & inpart) const {
+ iMode_=imode;
+ return modes_[imode]->generateDecay(inpart,this,inter,cc);
+}
+
void DecayIntegrator2::addMode(PhaseSpaceModePtr mode) const {
modes_.push_back(mode);
if(mode) mode->init();
}
ostream & Herwig::operator<<(ostream & os, const DecayIntegrator2 & decay) {
os << "The integrator has " << decay.modes_.size() << " modes" << endl;
for(unsigned int ix=0;ix<decay.modes_.size();++ix) {
os << "Information on mode " << ix << endl;
os << *(decay.modes_[ix]);
}
return os;
}
// reset the properities of all intermediates
void DecayIntegrator2::resetIntermediate(tcPDPtr part, Energy mass, Energy width) {
if(!part) return;
for(unsigned int ix=0,N=modes_.size();ix<N;++ix) {
modes_[ix]->resetIntermediate(part,mass,width);
}
}
diff --git a/Decay/DecayIntegrator2.h b/Decay/DecayIntegrator2.h
--- a/Decay/DecayIntegrator2.h
+++ b/Decay/DecayIntegrator2.h
@@ -1,477 +1,492 @@
// -*- 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.
- * 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
+ * @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 opt) const=0;
+ MEOption meopt) const = 0;
- virtual void constructSpinInfo(const Particle & part, ParticleVector outgoing) 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);
protected:
/**
* Methods to set variables in inheriting classes
*/
//@{
/**
* Set whether or not the intermediates are included
*/
void generateIntermediates(bool in) {generateInter_=in;}
/**
* 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/Makefile.am b/Decay/Makefile.am
--- a/Decay/Makefile.am
+++ b/Decay/Makefile.am
@@ -1,245 +1,248 @@
SUBDIRS = FormFactors Tau Baryon VectorMeson Perturbative \
WeakCurrents ScalarMeson TensorMeson Partonic General Radiation
if HAVE_EVTGEN
SUBDIRS += EvtGen
endif
noinst_LTLIBRARIES = libHwDecay.la
libHwDecay_la_LIBADD = \
$(top_builddir)/PDT/libHwPDT.la
nodist_libHwDecay_la_SOURCES = \
hwdecay__all.cc
BUILT_SOURCES = hwdecay__all.cc
CLEANFILES = hwdecay__all.cc
hwdecay__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 = \
DecayIntegrator.fh DecayIntegrator.h \
DecayIntegrator2.fh DecayIntegrator2.h \
DecayPhaseSpaceChannel.fh DecayPhaseSpaceChannel.h \
DecayPhaseSpaceMode.fh DecayPhaseSpaceMode.h \
PhaseSpaceMode.fh PhaseSpaceMode.h \
HwDecayerBase.fh HwDecayerBase.h \
HwDecayHandler.h \
DecayVertex.fh DecayVertex.h \
DecayMatrixElement.fh DecayMatrixElement.h \
TwoBodyDecayMatrixElement.h \
GeneralDecayMatrixElement.fh GeneralDecayMatrixElement.h \
BranchingRatioReweighter.h\
-PerturbativeDecayer.h ResonanceHelpers.h PhaseSpaceChannel.h
+PerturbativeDecayer.h \
+PerturbativeDecayer2.h \
+ResonanceHelpers.h PhaseSpaceChannel.h
DIR_CC_FILES = $(addprefix $(srcdir)/,$(ALL_CC_FILES))
ALL_CC_FILES = \
DecayIntegrator.cc \
DecayIntegrator2.cc \
DecayPhaseSpaceChannel.cc \
DecayPhaseSpaceMode.cc \
PhaseSpaceChannel.cc \
PhaseSpaceMode.cc \
HwDecayerBase.cc \
HwDecayHandler.cc \
DecayVertex.cc \
DecayMatrixElement.cc \
TwoBodyDecayMatrixElement.cc \
GeneralDecayMatrixElement.cc \
BranchingRatioReweighter.cc\
-PerturbativeDecayer.cc
+PerturbativeDecayer.cc\
+PerturbativeDecayer2.cc
##################
pkglib_LTLIBRARIES = Hw64Decay.la
Hw64Decay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 9:0:0
Hw64Decay_la_SOURCES = \
Hw64Decayer.h Hw64Decayer.cc
##################
pkglib_LTLIBRARIES += HwMamboDecay.la
HwMamboDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwMamboDecay_la_SOURCES = \
MamboDecayer.h MamboDecayer.cc
##################
noinst_LTLIBRARIES += libHwFormFactor.la
libHwFormFactor_la_SOURCES = \
FormFactors/BaryonFormFactor.cc \
FormFactors/ScalarFormFactor.cc \
FormFactors/BtoSGammaHadronicMass.cc \
FormFactors/BaryonFormFactor.fh \
FormFactors/BaryonFormFactor.h \
FormFactors/ScalarFormFactor.fh \
FormFactors/ScalarFormFactor.h \
FormFactors/BtoSGammaHadronicMass.h \
FormFactors/BtoSGammaHadronicMass.fh
pkglib_LTLIBRARIES += HwFormFactors.la
HwFormFactors_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwFormFactors_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/FormFactors
nodist_HwFormFactors_la_SOURCES = \
FormFactors/Formfactor__all.cc
##################
pkglib_LTLIBRARIES += HwTauDecay.la
HwTauDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwTauDecay_la_SOURCES = \
Tau/TauDecayer.cc\
Tau/TauDecayer2.cc
##################
pkglib_LTLIBRARIES += HwBaryonDecay.la
HwBaryonDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 9:0:0
HwBaryonDecay_la_LIBADD = \
$(top_builddir)/PDT/libHwBaryonWidth.la
HwBaryonDecay_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Baryon
nodist_HwBaryonDecay_la_SOURCES = \
Baryon/BaryonDecayer__all.cc
##################
pkglib_LTLIBRARIES += HwVMDecay.la
HwVMDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 9:0:0
HwVMDecay_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/VectorMeson
nodist_HwVMDecay_la_SOURCES = \
VectorMeson/VMDecayer__all.cc
##################
pkglib_LTLIBRARIES += HwPerturbativeDecay.la
HwPerturbativeDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwPerturbativeDecay_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Perturbative
nodist_HwPerturbativeDecay_la_SOURCES = \
Perturbative/Perturbative__all.cc
##################
noinst_LTLIBRARIES += libHwWeakCurrent.la
libHwWeakCurrent_la_SOURCES = \
WeakCurrents/WeakDecayCurrent.cc \
WeakCurrents/WeakCurrent.cc \
WeakCurrents/LeptonNeutrinoCurrent.cc \
WeakCurrents/WeakDecayCurrent.fh \
WeakCurrents/WeakDecayCurrent.h \
WeakCurrents/WeakCurrent.fh \
WeakCurrents/WeakCurrent.h \
WeakCurrents/LeptonNeutrinoCurrent.fh \
WeakCurrents/LeptonNeutrinoCurrent.h
pkglib_LTLIBRARIES += HwWeakCurrents.la
HwWeakCurrents_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwWeakCurrents_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/WeakCurrents
nodist_HwWeakCurrents_la_SOURCES = \
WeakCurrents/WeakCurrents__all.cc
##################
pkglib_LTLIBRARIES += HwSMDecay.la
HwSMDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 11:0:0
HwSMDecay_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/ScalarMeson
nodist_HwSMDecay_la_SOURCES = \
ScalarMeson/SMDecayer__all.cc
##################
pkglib_LTLIBRARIES += HwTMDecay.la
HwTMDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 9:0:0
HwTMDecay_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/TensorMeson
nodist_HwTMDecay_la_SOURCES = \
TensorMeson/TMDecayer__all.cc
##################
pkglib_LTLIBRARIES += HwPartonicDecay.la
HwPartonicDecay_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 10:0:0
HwPartonicDecay_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Partonic
nodist_HwPartonicDecay_la_SOURCES = \
Partonic/Partonic__all.cc
##################
noinst_LTLIBRARIES += libHwDecRad.la
libHwDecRad_la_SOURCES = \
Radiation/DecayRadiationGenerator.cc \
Radiation/QEDRadiationHandler.cc \
Radiation/DecayRadiationGenerator.h \
Radiation/DecayRadiationGenerator.fh \
Radiation/QEDRadiationHandler.fh \
Radiation/QEDRadiationHandler.h
pkglib_LTLIBRARIES += HwSOPHTY.la
HwSOPHTY_la_LDFLAGS = \
$(AM_LDFLAGS) -module -version-info 4:0:0
HwSOPHTY_la_CPPFLAGS = \
$(AM_CPPFLAGS) -I$(srcdir)/Radiation
nodist_HwSOPHTY_la_SOURCES = \
Radiation/Sophty__all.cc
##################
diff --git a/Decay/Perturbative/SMHiggsFermionsDecayer.cc b/Decay/Perturbative/SMHiggsFermionsDecayer.cc
--- a/Decay/Perturbative/SMHiggsFermionsDecayer.cc
+++ b/Decay/Perturbative/SMHiggsFermionsDecayer.cc
@@ -1,395 +1,398 @@
// -*- C++ -*-
//
// SMHiggsFermionsDecayer.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 SMHiggsFermionsDecayer class.
//
#include "SMHiggsFermionsDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Decay/DecayVertex.h"
#include "ThePEG/Helicity/ScalarSpinInfo.h"
#include "ThePEG/Helicity/FermionSpinInfo.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include "Herwig/Utilities/Maths.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Shower/ShowerAlpha.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
SMHiggsFermionsDecayer::SMHiggsFermionsDecayer() :
CF_(4./3.), NLO_(false) {
_maxwgt.resize(9);
_maxwgt[0]=0.;
_maxwgt[1]=0;
_maxwgt[2]=0;
_maxwgt[3]=0.0194397;
_maxwgt[4]=0.463542;
_maxwgt[5]=0.;
_maxwgt[6]=6.7048e-09;
_maxwgt[7]=0.00028665;
_maxwgt[8]=0.0809643;
}
void SMHiggsFermionsDecayer::doinit() {
- PerturbativeDecayer::doinit();
+ PerturbativeDecayer2::doinit();
// get the vertices from the Standard Model object
tcHwSMPtr hwsm=dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm)
throw InitException() << "SMHiggsFermionsDecayer needs the StandardModel class"
<< " to be either the Herwig one or a class inheriting"
<< " from it";
_hvertex = hwsm->vertexFFH();
// make sure they are initialized
_hvertex->init();
// get the width generator for the higgs
tPDPtr higgs = getParticleData(ParticleID::h0);
// set up the decay modes
- vector<double> wgt(0);
unsigned int imode=0;
- tPDVector extpart(3);
- DecayPhaseSpaceModePtr mode;
- int iy;
- extpart[0]=higgs;
for(unsigned int istep=0;istep<11;istep+=10) {
for(unsigned ix=1;ix<7;++ix) {
if(istep<10||ix%2!=0) {
- iy = ix+istep;
- extpart[1]=getParticleData( iy);
- extpart[2]=getParticleData(-iy);
- mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
- addMode(mode,_maxwgt[imode],wgt);
+ int iy = ix+istep;
+ tPDVector out = {getParticleData( iy),
+ getParticleData(-iy)};
+ PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(higgs,out,_maxwgt[imode]));
++imode;
}
}
}
// Energy quarkMass = getParticleData(ParticleID::b )->mass();
// Energy higgsMass = getParticleData(ParticleID::h0)->mass();
// double mu = quarkMass/higgsMass;
// double beta = sqrt(1.-4.*sqr(mu));
// double beta2 = sqr(beta);
// double aS = SM().alphaS(sqr(higgsMass));
// double L = log((1.+beta)/(1.-beta));
// cerr << "testing " << beta << " " << mu << "\n";
// cerr << "testing " << aS << " " << L << "\n";
// double fact =
// 6.-0.75*(1.+beta2)/beta2+12.*log(mu)-8.*log(beta)
// +(5./beta-2.*beta+0.375*sqr(1.-beta2)/beta2/beta)*L
// +(1.+beta2)/beta*(4.*L*log(0.5*(1.+beta)/beta)
// -2.*log(0.5*(1.+beta))*log(0.5*(1.-beta))
// +8.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))
// -4.*Herwig::Math::ReLi2(0.5*(1.-beta)));
// cerr << "testing correction "
// << 1.+4./3.*aS/Constants::twopi*fact
// << "\n";
// double real = 4./3.*aS/Constants::twopi*
// (8.-0.75*(1.+beta2)/beta2+8.*log(mu)-8.*log(beta)
// +(3./beta+0.375*sqr(1.-beta2)/pow(beta,3))*L
// +(1.+beta2)/beta*(-0.5*sqr(L)+4.*L*log(0.5*(1.+beta))
// -2.*L*log(beta)-2.*log(0.5*(1.+beta))*log(0.5*(1.-beta))
// +6.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))
// -4.*Herwig::Math::ReLi2(0.5*(1.-beta))
// -2./3.*sqr(Constants::pi)));
// double virt = 4./3.*aS/Constants::twopi*
// (-2.+4.*log(mu)+(2./beta-2.*beta)*L
// +(1.+beta2)/beta*(0.5*sqr(L)-2.*L*log(beta)+2.*sqr(Constants::pi)/3.
// +2.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))));
// cerr << "testing real " << real << "\n";
// cerr << "testing virtual " << virt << "\n";
// cerr << "testing total no mb corr " << 1.+real+virt << "\n";
// cerr << "testing total mb corr " << 1.+real+virt +(8./3. - 2.*log(sqr(mu)))*aS/Constants::pi << "\n";
// InvEnergy2 Gf = 1.166371e-5/GeV2;
// Gf = sqrt(2.)*4*Constants::pi*SM().alphaEM(sqr(higgsMass))/8./SM().sin2ThetaW()/
// sqr(getParticleData(ParticleID::Wplus)->mass());
// cerr << "testing GF " << Gf*GeV2 << "\n";
// Energy LO = (3./8./Constants::pi)*sqrt(2)*sqr(quarkMass)*Gf*higgsMass*beta*beta*beta;
// cerr << "testing LO " << LO/GeV << "\n";
// cerr << "testing quark mass " << quarkMass/GeV << "\n";
// cerr << "testing gamma " << (1.+real+virt)*LO/MeV << "\n";
}
bool SMHiggsFermionsDecayer::accept(tcPDPtr parent, const tPDVector & children) const {
if(parent->id()!=ParticleID::h0||children.size()!=2) return false;
tPDVector::const_iterator pit = children.begin();
int id1=(**pit).id();
++pit;
int id2=(**pit).id();
if(id1==-id2&&(abs(id1)<=6||(abs(id1)>=11&&abs(id1)<=16)))
return true;
else
return false;
}
ParticleVector SMHiggsFermionsDecayer::decay(const Particle & parent,
const tPDVector & children) const {
// id's of the decaying particles
tPDVector::const_iterator pit(children.begin());
int id1((**pit).id());
int imode=-1;
if(abs(id1)<=6) imode = abs(id1)-1;
else if(abs(id1)>=11&&abs(id1)<=16) imode = (abs(id1)-11)/2+6;
ParticleVector output(generate(false,false,imode,parent));
// set up the colour flow
if(output[0]->hasColour()) output[0]->antiColourNeighbour(output[1]);
else if(output[1]->hasColour()) output[1]->antiColourNeighbour(output[0]);
return output;
}
void SMHiggsFermionsDecayer::persistentOutput(PersistentOStream & os) const {
os << _maxwgt << _hvertex << NLO_;
}
void SMHiggsFermionsDecayer::persistentInput(PersistentIStream & is, int) {
is >> _maxwgt >> _hvertex >> NLO_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<SMHiggsFermionsDecayer,PerturbativeDecayer>
+DescribeClass<SMHiggsFermionsDecayer,PerturbativeDecayer2>
describeHerwigSMHiggsFermionsDecayer("Herwig::SMHiggsFermionsDecayer", "HwPerturbativeHiggsDecay.so");
void SMHiggsFermionsDecayer::Init() {
static ClassDocumentation<SMHiggsFermionsDecayer> documentation
("The SMHiggsFermionsDecayer class implements the decat of the Standard Model"
" Higgs boson to the Standard Model fermions");
static ParVector<SMHiggsFermionsDecayer,double> interfaceMaxWeights
("MaxWeights",
"Maximum weights for the various decays",
&SMHiggsFermionsDecayer::_maxwgt, 9, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static Switch<SMHiggsFermionsDecayer,bool> interfaceNLO
("NLO",
"Whether to return the LO or NLO result",
&SMHiggsFermionsDecayer::NLO_, false, false, false);
static SwitchOption interfaceNLOLO
(interfaceNLO,
"No",
"Leading-order result",
false);
static SwitchOption interfaceNLONLO
(interfaceNLO,
"Yes",
"NLO result",
true);
}
+void SMHiggsFermionsDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ int iferm(1),ianti(0);
+ if(decay[0]->id()>0) swap(iferm,ianti);
+ ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
+ incoming,true);
+ SpinorBarWaveFunction::
+ constructSpinInfo(_wavebar,decay[iferm],outgoing,true);
+ SpinorWaveFunction::
+ constructSpinInfo(_wave ,decay[ianti],outgoing,true);
+}
// return the matrix element squared
-double SMHiggsFermionsDecayer::me2(const int, const Particle & part,
- const ParticleVector & decay,
+double SMHiggsFermionsDecayer::me2(const int,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half)));
int iferm(1),ianti(0);
- if(decay[0]->id()>0) swap(iferm,ianti);
+ if(outgoing[0]->id()>0) swap(iferm,ianti);
if(meopt==Initialize) {
ScalarWaveFunction::
calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
_swave = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming);
// fix rho if no correlations
fixRho(_rho);
}
- if(meopt==Terminate) {
- ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
- incoming,true);
- SpinorBarWaveFunction::
- constructSpinInfo(_wavebar,decay[iferm],outgoing,true);
- SpinorWaveFunction::
- constructSpinInfo(_wave ,decay[ianti],outgoing,true);
- return 0.;
+ SpinorBarWaveFunction wbar(momenta[iferm],outgoing[iferm],Helicity::outgoing);
+ SpinorWaveFunction w (momenta[ianti],outgoing[ianti],Helicity::outgoing);
+ _wavebar.resize(2);
+ _wave.resize(2);
+ for(unsigned int ihel=0;ihel<2;++ihel) {
+ wbar.reset(ihel);
+ _wavebar[ihel] = wbar;
+ w.reset(ihel);
+ _wave[ihel] = w;
}
- SpinorBarWaveFunction::
- calculateWaveFunctions(_wavebar,decay[iferm],outgoing);
- SpinorWaveFunction::
- calculateWaveFunctions(_wave ,decay[ianti],outgoing);
Energy2 scale(sqr(part.mass()));
unsigned int ifm,ia;
for(ifm=0;ifm<2;++ifm) {
for(ia=0;ia<2;++ia) {
if(iferm>ianti)
(*ME())(0,ia,ifm)=_hvertex->evaluate(scale,_wave[ia],
_wavebar[ifm],_swave);
else
(*ME())(0,ifm,ia)=_hvertex->evaluate(scale,_wave[ia],
_wavebar[ifm],_swave);
}
}
- int id = abs(decay[0]->id());
+ int id = abs(outgoing[0]->id());
double output=(ME()->contract(_rho)).real()*UnitRemoval::E2/scale;
if(id <=6) output*=3.;
// test of the partial width
-// Ptr<Herwig::StandardModel>::transient_const_pointer
-// hwsm=dynamic_ptr_cast<Ptr<Herwig::StandardModel>::transient_const_pointer>(standardModel());
-// double g2(hwsm->alphaEM(scale)*4.*Constants::pi/hwsm->sin2ThetaW());
-// Energy mass(hwsm->mass(scale,decay[0]->dataPtr())),
-// mw(getParticleData(ParticleID::Wplus)->mass());
-// double beta(sqrt(1.-4.*decay[0]->mass()*decay[0]->mass()/scale));
-// cerr << "testing alpha " << hwsm->alphaEM(scale) << "\n";
-// Energy test(g2*mass*mass*beta*beta*beta*part.mass()/32./Constants::pi/mw/mw);
-// if(abs(decay[0]->id())<=6){test *=3.;}
-// cout << "testing the answer " << output << " "
-// << test/GeV
-// << endl;
+ // Ptr<Herwig::StandardModel>::transient_const_pointer
+ // hwsm=dynamic_ptr_cast<Ptr<Herwig::StandardModel>::transient_const_pointer>(standardModel());
+ // double g2(hwsm->alphaEM(scale)*4.*Constants::pi/hwsm->sin2ThetaW());
+ // Energy mass(hwsm->mass(scale,outgoing[0])),
+ // mw(getParticleData(ParticleID::Wplus)->mass());
+ // double beta(sqrt(1.-4.*decay[0]->mass()*decay[0]->mass()/scale));
+ // cerr << "testing alpha " << hwsm->alphaEM(scale) << "\n";
+ // Energy test(g2*mass*mass*beta*beta*beta*part.mass()/32./Constants::pi/mw/mw);
+ // if(abs(decay[0]->id())<=6){test *=3.;}
+ // cout << "testing the answer " << output << " "
+ // << test/GeV
+ // << endl;
// leading-order result
if(!NLO_) return output;
// fermion mass
- Energy particleMass = decay[0]->dataPtr()->mass();
+ Energy particleMass = outgoing[0]->mass();
// check decay products coloured, otherwise return
- if(!decay[0]->dataPtr()->coloured()||
+ if(!outgoing[0]->coloured()||
particleMass==ZERO) return output;
// inital masses, couplings etc
// higgs mass
mHiggs_ = part.mass();
// strong coupling
aS_ = SM().alphaS(sqr(mHiggs_));
// reduced mass
mu_ = particleMass/mHiggs_;
mu2_ = sqr(mu_);
// generate y
double yminus = 0.;
double yplus = 1.-2.*mu_*(1.-mu_)/(1.-2*mu2_);
double y = yminus + UseRandom::rnd()*(yplus-yminus);
//generate z for D31,2
double v = sqrt(sqr(2.*mu2_+(1.-2.*mu2_)*(1.-y))-4.*mu2_)/(1.-2.*mu2_)/(1.-y);
double zplus = (1.+v)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
double zminus = (1.-v)*(1.-2.*mu2_)*y/2./(mu2_ +(1.-2.*mu2_)*y);
double z = zminus + UseRandom::rnd()*(zplus-zminus);
// map y,z to x1,x2 for both possible emissions
double x2 = 1. - y*(1.-2.*mu2_);
double x1 = 1. - z*(x2-2.*mu2_);
//get the dipoles
InvEnergy2 D1 = dipoleSubtractionTerm( x1, x2);
InvEnergy2 D2 = dipoleSubtractionTerm( x2, x1);
InvEnergy2 dipoleSum = abs(D1) + abs(D2);
//jacobian
double jac = (1.-y)*(yplus-yminus)*(zplus-zminus);
//calculate real
Energy2 realPrefactor = 0.25*sqr(mHiggs_)*sqr(1.-2.*mu2_)
/sqrt(calculateLambda(1,mu2_,mu2_))/sqr(Constants::twopi);
InvEnergy2 realEmission = 4.*Constants::pi*aS_*CF_*calculateRealEmission( x1, x2);
// calculate the virtual
double virtualTerm = calculateVirtualTerm();
// running mass correction
virtualTerm += (8./3. - 2.*log(mu2_))*aS_/Constants::pi;
//answer = (born + virtual + real)/born * LO
output *= 1. + virtualTerm + 2.*jac*realPrefactor*(realEmission*abs(D1)/dipoleSum - D1);
// return the answer
return output;
}
void SMHiggsFermionsDecayer::dataBaseOutput(ofstream & os,bool header) const {
if(header) os << "update decayers set parameters=\"";
- // parameters for the PerturbativeDecayer base class
+ // parameters for the PerturbativeDecayer2 base class
for(unsigned int ix=0;ix<_maxwgt.size();++ix) {
os << "newdef " << name() << ":MaxWeights " << ix << " "
<< _maxwgt[ix] << "\n";
}
- PerturbativeDecayer::dataBaseOutput(os,false);
+ PerturbativeDecayer2::dataBaseOutput(os,false);
if(header) os << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
void SMHiggsFermionsDecayer::doinitrun() {
- PerturbativeDecayer::doinitrun();
+ PerturbativeDecayer2::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<numberModes();++ix) {
_maxwgt[ix] = mode(ix)->maxWeight();
}
}
}
//calculate lambda
double SMHiggsFermionsDecayer::calculateLambda(double x, double y, double z) const{
return sqr(x)+sqr(y)+sqr(z)-2.*x*y-2.*x*z-2.*y*z;
}
//calculates the dipole subtraction term for x1, D31,2 (Dij,k),
// 2 is the spectator anti-fermion and 3 is the gluon
InvEnergy2 SMHiggsFermionsDecayer::
dipoleSubtractionTerm(double x1, double x2) const{
InvEnergy2 commonPrefactor = CF_*8.*Constants::pi*aS_/sqr(mHiggs_);
return commonPrefactor/(1.-x2)*
(2.*(1.-2.*mu2_)/(2.-x1-x2)-
sqrt((1.-4.*mu2_)/(sqr(x2)-4.*mu2_))*
(x2-2.*mu2_)*(2.+(x1-1.)/(x2-2.*mu2_)+2.*mu2_/(1.-x2))/(1.-2.*mu2_));
}
//return ME for real emission
InvEnergy2 SMHiggsFermionsDecayer::
calculateRealEmission(double x1, double x2) const {
InvEnergy2 prefactor = 2./sqr(mHiggs_)/(1.-4.*mu2_);
return prefactor*(2. + (1.-x1)/(1.-x2) + (1.-x2)/(1.-x1)
+ 2.*(1.-2.*mu2_)*(1.-4.*mu2_)/(1.-x1)/(1.-x2)
- 2.*(1.-4.*mu2_)*(1./(1.-x2)+1./(1.-x1))
- 2.*mu2_*(1.-4.*mu2_)*(1./sqr(1.-x2)+1./sqr(1.-x1)));
}
double SMHiggsFermionsDecayer::
calculateVirtualTerm() const {
// logs and prefactors
double beta = sqrt(1.-4.*mu2_);
double L = log((1.+beta)/(1.-beta));
double prefactor = CF_*aS_/Constants::twopi;
// non-singlet piece
double nonSingletTerm = calculateNonSingletTerm(beta, L);
double virtualTerm =
-2.+4.*log(mu_)+(2./beta - 2.*beta)*L
+ (2.-4.*mu2_)/beta*(0.5*sqr(L) - 2.*L*log(beta)
+ 2.*Herwig::Math::ReLi2((1.-beta)/(1.+beta))
+ 2.*sqr(Constants::pi)/3.);
double iEpsilonTerm =
2.*(3.-sqr(Constants::pi)/2. + 0.5*log(mu2_) - 1.5*log(1.-2.*mu2_)
-(1.-2.*mu2_)/beta*(0.5*sqr(L)+sqr(Constants::pi)/6.
-2.*L*log(1.-2.*mu2_))
+ nonSingletTerm);
return prefactor*(virtualTerm+iEpsilonTerm);
}
//non-singlet piece of I(epsilon) insertion operator
double SMHiggsFermionsDecayer::
calculateNonSingletTerm(double beta, double L) const {
return 1.5*log(1.-2.*mu2_)
+ (1.-2.*mu2_)/beta*(- 2.*L*log(4.*(1.-2.*mu2_)/sqr(1.+beta))+
+ 2.*Herwig::Math::ReLi2(sqr((1.-beta)/(1.+beta)))
- 2.*Herwig::Math::ReLi2(2.*beta/(1.+beta))
- sqr(Constants::pi)/6.)
+ log(1.-mu_)
- 2.*log(1.-2.*mu_)
- 2.*mu2_/(1.-2.*mu2_)*log(mu_/(1.-mu_))
- mu_/(1.-mu_)
+ 2.*mu_*(2*mu_-1.)/(1.-2.*mu2_)
+ 0.5*sqr(Constants::pi);
}
double SMHiggsFermionsDecayer::matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
const ParticleVector & decay3, MEOption,
ShowerInteraction inter) {
mHiggs_ = inpart.mass();
mu_ = decay2[0]->mass()/mHiggs_;
mu2_ = sqr(mu_);
double x1 = 2.*decay3[0]->momentum().t()/mHiggs_;
double x2 = 2.*decay3[1]->momentum().t()/mHiggs_;
double pre = inter==ShowerInteraction::QCD ? CF_ : sqr(double(decay2[0]->dataPtr()->iCharge())/3.);
return pre*calculateRealEmission(x1,x2)*4.*Constants::pi*sqr(mHiggs_);
}
diff --git a/Decay/Perturbative/SMHiggsFermionsDecayer.h b/Decay/Perturbative/SMHiggsFermionsDecayer.h
--- a/Decay/Perturbative/SMHiggsFermionsDecayer.h
+++ b/Decay/Perturbative/SMHiggsFermionsDecayer.h
@@ -1,279 +1,289 @@
// -*- C++ -*-
//
// SMHiggsFermionsDecayer.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_SMHiggsFermionsDecayer_H
#define HERWIG_SMHiggsFermionsDecayer_H
//
// This is the declaration of the SMHiggsFermionsDecayer class.
//
-#include "Herwig/Decay/PerturbativeDecayer.h"
+#include "Herwig/Decay/PerturbativeDecayer2.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
-#include "Herwig/Decay/DecayPhaseSpaceMode.h"
+#include "Herwig/Decay/PhaseSpaceMode.h"
namespace Herwig {
using namespace ThePEG;
/**
* The SMHiggsFermionsDecayer class is designed to decay the Standard Model Higgs
* to the Standard Model fermions.
*
- * @see PerturbativeDecayer
+ * @see PerturbativeDecayer2
*/
-class SMHiggsFermionsDecayer: public PerturbativeDecayer {
+class SMHiggsFermionsDecayer: public PerturbativeDecayer2 {
public:
/**
* The default constructor.
*/
SMHiggsFermionsDecayer();
/**
* Which of the possible decays is required
*/
virtual int modeNumber(bool & , tcPDPtr , const tPDVector & ) const {return -1;}
/**
* 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;
/**
* 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;
+
/**
* 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 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 ParticleVector & decay, MEOption meopt) const;
+ double me2(const int ichan,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ MEOption meopt) const;
+
+ /**
+ * Construct the SpinInfos for the particles produced in the decay
+ */
+ virtual void constructSpinInfo(const Particle & part,
+ ParticleVector outgoing) const;
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
*/
virtual void dataBaseOutput(ofstream & os,bool header) const;
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {return FSR;}
/**
* Calculate matrix element ratio R/B
*/
virtual double matrixElementRatio(const Particle & inpart, const ParticleVector & decay2,
const ParticleVector & decay3, MEOption meopt,
ShowerInteraction inter);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {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:
/**
* Calcluate the Kallen function
*/
double calculateLambda(double x, double y, double z) const;
/**
* Dipole subtraction term
*/
InvEnergy2 dipoleSubtractionTerm(double x1, double x2) const;
/**
* Real emission term
*/
InvEnergy2 calculateRealEmission(double x1, double x2) const;
/**
* Virtual term
*/
double calculateVirtualTerm() const;
/**
* Non-singlet term
*/
double calculateNonSingletTerm(double beta, double L) const;
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SMHiggsFermionsDecayer & operator=(const SMHiggsFermionsDecayer &);
private:
/**
* Pointer to the Higgs vertex
*/
AbstractFFSVertexPtr _hvertex;
/**
* maximum weights for the different decay modes
*/
vector<double> _maxwgt;
/**
* Spin density matrix
*/
mutable RhoDMatrix _rho;
/**
* Scalar wavefunction
*/
mutable ScalarWaveFunction _swave;
/**
* Spinor wavefunction
*/
mutable vector<SpinorWaveFunction> _wave;
/**
* Barred spinor wavefunction
*/
mutable vector<SpinorBarWaveFunction> _wavebar;
private:
/**
* The colour factor
*/
double CF_;
/**
* The Higgs mass
*/
mutable Energy mHiggs_;
/**
* The reduced mass
*/
mutable double mu_;
/**
* The square of the reduced mass
*/
mutable double mu2_;
/**
* The strong coupling
*/
mutable double aS_;
/**
* Stuff for the POWHEG correction
*/
//@{
/**
* The ParticleData objects for the fermions
*/
vector<tcPDPtr> partons_;
/**
* The fermion momenta
*/
vector<Lorentz5Momentum> quark_;
/**
* The momentum of the radiated gauge boson
*/
Lorentz5Momentum gauge_;
/**
* The Higgs boson
*/
PPtr higgs_;
/**
* Higgs mass squared
*/
Energy2 mh2_;
//@}
/**
* LO or NLO ?
*/
bool NLO_;
};
}
#endif /* HERWIG_SMHiggsFermionsDecayer_H */
diff --git a/Decay/PerturbativeDecayer.cc b/Decay/PerturbativeDecayer.cc
--- a/Decay/PerturbativeDecayer.cc
+++ b/Decay/PerturbativeDecayer.cc
@@ -1,1172 +1,1172 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PerturbativeDecayer class.
//
#include "PerturbativeDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/EnumIO.h"
using namespace Herwig;
void PerturbativeDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(pTmin_,GeV) << oenum(inter_) << alphaS_ << alphaEM_
<< useMEforT2_ << C_ << ymax_ << phaseOpt_;
}
void PerturbativeDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(pTmin_,GeV) >> ienum(inter_) >> alphaS_ >> alphaEM_
>> useMEforT2_ >> C_ >> ymax_ >> phaseOpt_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeAbstractClass<PerturbativeDecayer,DecayIntegrator>
describeHerwigPerturbativeDecayer("Herwig::PerturbativeDecayer",
"Herwig.so HwPerturbativeDecay.so");
void PerturbativeDecayer::Init() {
static ClassDocumentation<PerturbativeDecayer> documentation
("The PerturbativeDecayer class is the mase class for "
"perturbative decays in Herwig");
static Parameter<PerturbativeDecayer,Energy> interfacepTmin
("pTmin",
"Minimum transverse momentum from gluon radiation",
&PerturbativeDecayer::pTmin_, GeV, 1.0*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Switch<PerturbativeDecayer,ShowerInteraction> interfaceInteractions
("Interactions",
"which interactions to include for the hard corrections",
&PerturbativeDecayer::inter_, ShowerInteraction::QCD, false, false);
static SwitchOption interfaceInteractionsQCD
(interfaceInteractions,
"QCD",
"QCD Only",
ShowerInteraction::QCD);
static SwitchOption interfaceInteractionsQED
(interfaceInteractions,
"QED",
"QED only",
ShowerInteraction::QED);
static SwitchOption interfaceInteractionsQCDandQED
(interfaceInteractions,
"QCDandQED",
"Both QCD and QED",
ShowerInteraction::Both);
static Reference<PerturbativeDecayer,ShowerAlpha> interfaceAlphaS
("AlphaS",
"Object for the coupling in the generation of hard QCD radiation",
&PerturbativeDecayer::alphaS_, false, false, true, true, false);
static Reference<PerturbativeDecayer,ShowerAlpha> interfaceAlphaEM
("AlphaEM",
"Object for the coupling in the generation of hard QED radiation",
&PerturbativeDecayer::alphaEM_, false, false, true, true, false);
static Switch<PerturbativeDecayer,bool> interfaceUseMEForT2
("UseMEForT2",
"Use the matrix element correction, if available to fill the T2"
" region for the decay shower and don't fill using the shower",
&PerturbativeDecayer::useMEforT2_, true, false, false);
static SwitchOption interfaceUseMEForT2Shower
(interfaceUseMEForT2,
"Shower",
"Use the shower to fill the T2 region",
false);
static SwitchOption interfaceUseMEForT2ME
(interfaceUseMEForT2,
"ME",
"Use the Matrix element to fill the T2 region",
true);
static Parameter<PerturbativeDecayer,double> interfacePrefactor
("Prefactor",
"The prefactor for the sampling of the powheg Sudakov",
&PerturbativeDecayer::C_, 6.3, 0.0, 1e10,
false, false, Interface::limited);
static Parameter<PerturbativeDecayer,double> interfaceYMax
("YMax",
"The maximum value for the rapidity",
&PerturbativeDecayer::ymax_, 10., 0.0, 100.,
false, false, Interface::limited);
static Switch<PerturbativeDecayer,unsigned int> interfacePhaseSpaceOption
("PhaseSpaceOption",
"Option for the phase-space sampling",
&PerturbativeDecayer::phaseOpt_, 0, false, false);
static SwitchOption interfacePhaseSpaceOptionFixedYLimits
(interfacePhaseSpaceOption,
"FixedYLimits",
"Use a fixed limit for the rapidity",
0);
static SwitchOption interfacePhaseSpaceOptionVariableYLimits
(interfacePhaseSpaceOption,
"VariableYLimits",
"Change limit for the rapidity with pT",
1);
}
double PerturbativeDecayer::matrixElementRatio(const Particle & ,
const ParticleVector & ,
const ParticleVector & ,
MEOption ,
ShowerInteraction ) {
throw Exception() << "Base class PerturbativeDecayer::matrixElementRatio() "
<< "called, should have an implementation in the inheriting class"
<< Exception::runerror;
return 0.;
}
RealEmissionProcessPtr PerturbativeDecayer::generateHardest(RealEmissionProcessPtr born) {
return getHardEvent(born,false,inter_);
}
RealEmissionProcessPtr PerturbativeDecayer::applyHardMatrixElementCorrection(RealEmissionProcessPtr born) {
return getHardEvent(born,true,ShowerInteraction::QCD);
}
RealEmissionProcessPtr PerturbativeDecayer::getHardEvent(RealEmissionProcessPtr born,
bool inDeadZone,
ShowerInteraction inter) {
// check one incoming
assert(born->bornIncoming().size()==1);
// check exactly two outgoing particles
assert(born->bornOutgoing().size()==2); // search for coloured particles
bool colouredParticles=born->bornIncoming()[0]->dataPtr()->coloured();
bool chargedParticles=born->bornIncoming()[0]->dataPtr()->charged();
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
if(born->bornOutgoing()[ix]->dataPtr()->coloured())
colouredParticles=true;
if(born->bornOutgoing()[ix]->dataPtr()->charged())
chargedParticles=true;
}
// if no coloured/charged particles return
if ( !colouredParticles && !chargedParticles ) return RealEmissionProcessPtr();
if ( !colouredParticles && inter==ShowerInteraction::QCD ) return RealEmissionProcessPtr();
if ( ! chargedParticles && inter==ShowerInteraction::QED ) return RealEmissionProcessPtr();
// for decay b -> a c
// set progenitors
PPtr cProgenitor = born->bornOutgoing()[0];
PPtr aProgenitor = born->bornOutgoing()[1];
// get the decaying particle
PPtr bProgenitor = born->bornIncoming()[0];
// identify which dipoles are required
vector<DipoleType> dipoles;
if(!identifyDipoles(dipoles,aProgenitor,bProgenitor,cProgenitor,inter)) {
return RealEmissionProcessPtr();
}
Energy trialpT = pTmin_;
LorentzRotation eventFrame;
vector<Lorentz5Momentum> momenta;
vector<Lorentz5Momentum> trialMomenta(4);
PPtr finalEmitter, finalSpectator;
PPtr trialEmitter, trialSpectator;
DipoleType finalType(FFa,ShowerInteraction::QCD);
for (int i=0; i<int(dipoles.size()); ++i) {
if(dipoles[i].type==FFg) continue;
// assign emitter and spectator based on current dipole
if (dipoles[i].type==FFc || dipoles[i].type==IFc || dipoles[i].type==IFbc) {
trialEmitter = cProgenitor;
trialSpectator = aProgenitor;
}
else if (dipoles[i].type==FFa || dipoles[i].type==IFa || dipoles[i].type==IFba) {
trialEmitter = aProgenitor;
trialSpectator = cProgenitor;
}
// find rotation from lab to frame with the spectator along -z
LorentzRotation trialEventFrame(bProgenitor->momentum().findBoostToCM());
Lorentz5Momentum pspectator = (trialEventFrame*trialSpectator->momentum());
trialEventFrame.rotateZ( -pspectator.phi() );
trialEventFrame.rotateY( -pspectator.theta() - Constants::pi );
// invert it
trialEventFrame.invert();
// try to generate an emission
pT_ = pTmin_;
vector<Lorentz5Momentum> trialMomenta
= hardMomenta(bProgenitor, trialEmitter, trialSpectator,
dipoles, i, inDeadZone);
// select dipole which gives highest pT emission
if(pT_>trialpT) {
trialpT = pT_;
momenta = trialMomenta;
eventFrame = trialEventFrame;
finalEmitter = trialEmitter;
finalSpectator = trialSpectator;
finalType = dipoles[i];
if (dipoles[i].type==FFc || dipoles[i].type==FFa ) {
if((momenta[3]+momenta[1]).m2()-momenta[1].m2()>
(momenta[3]+momenta[2]).m2()-momenta[2].m2()) {
swap(finalEmitter,finalSpectator);
swap(momenta[1],momenta[2]);
}
}
}
}
pT_ = trialpT;
// if no emission return
if(momenta.empty()) {
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QCD)
born->pT()[ShowerInteraction::QCD] = pTmin_;
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QED)
born->pT()[ShowerInteraction::QED] = pTmin_;
return born;
}
// rotate momenta back to the lab
for(unsigned int ix=0;ix<momenta.size();++ix) {
momenta[ix] *= eventFrame;
}
// set maximum pT for subsequent branchings
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QCD)
born->pT()[ShowerInteraction::QCD] = pT_;
if(inter==ShowerInteraction::Both || inter==ShowerInteraction::QED)
born->pT()[ShowerInteraction::QED] = pT_;
// get ParticleData objects
tcPDPtr b = bProgenitor ->dataPtr();
tcPDPtr e = finalEmitter ->dataPtr();
tcPDPtr s = finalSpectator->dataPtr();
tcPDPtr boson = getParticleData(finalType.interaction==ShowerInteraction::QCD ?
ParticleID::g : ParticleID::gamma);
// create new ShowerParticles
PPtr emitter = e ->produceParticle(momenta[1]);
PPtr spectator = s ->produceParticle(momenta[2]);
PPtr gauge = boson->produceParticle(momenta[3]);
PPtr incoming = b ->produceParticle(bProgenitor->momentum());
// insert the particles
born->incoming().push_back(incoming);
unsigned int iemit(0),ispect(0);
for(unsigned int ix=0;ix<born->bornOutgoing().size();++ix) {
if(born->bornOutgoing()[ix]==finalEmitter) {
born->outgoing().push_back(emitter);
iemit = born->outgoing().size();
}
else if(born->bornOutgoing()[ix]==finalSpectator) {
born->outgoing().push_back(spectator);
ispect = born->outgoing().size();
}
}
born->outgoing().push_back(gauge);
if(!spectator->dataPtr()->coloured() ||
(finalType.type != FFa && finalType.type!=FFc) ) ispect = 0;
born->emitter(iemit);
born->spectator(ispect);
born->emitted(3);
// boost if being use as ME correction
if(inDeadZone) {
if(finalType.type==IFa || finalType.type==IFba) {
LorentzRotation trans(cProgenitor->momentum().findBoostToCM());
trans.boost(spectator->momentum().boostVector());
born->transformation(trans);
}
else if(finalType.type==IFc || finalType.type==IFbc) {
LorentzRotation trans(bProgenitor->momentum().findBoostToCM());
trans.boost(spectator->momentum().boostVector());
born->transformation(trans);
}
}
// set the interaction
born->interaction(finalType.interaction);
// set up colour lines
getColourLines(born);
// return the tree
return born;
}
bool PerturbativeDecayer::identifyDipoles(vector<DipoleType> & dipoles,
PPtr & aProgenitor,
PPtr & bProgenitor,
PPtr & cProgenitor,
ShowerInteraction inter) const {
enhance_ = 1.;
// identify any QCD dipoles
if(inter==ShowerInteraction::QCD ||
inter==ShowerInteraction::Both) {
PDT::Colour bColour = bProgenitor->dataPtr()->iColour();
PDT::Colour cColour = cProgenitor->dataPtr()->iColour();
PDT::Colour aColour = aProgenitor->dataPtr()->iColour();
// decaying colour singlet
if (bColour==PDT::Colour0 ) {
if ((cColour==PDT::Colour3 && aColour==PDT::Colour3bar) ||
(cColour==PDT::Colour3bar && aColour==PDT::Colour3) ||
(cColour==PDT::Colour8 && aColour==PDT::Colour8)){
dipoles.push_back(DipoleType(FFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc,ShowerInteraction::QCD));
if(aProgenitor->id()==ParticleID::g &&
cProgenitor->id()==ParticleID::g ) {
enhance_ = 1.5;
dipoles.push_back(DipoleType(FFg,ShowerInteraction::QCD));
}
}
}
// decaying colour triplet
else if (bColour==PDT::Colour3 ) {
if (cColour==PDT::Colour3 && aColour==PDT::Colour0){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour0 && aColour==PDT::Colour3){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour3){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour3 && aColour==PDT::Colour8){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
else if(cColour==PDT::Colour3bar && aColour==PDT::Colour3bar) {
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
}
// decaying colour anti-triplet
else if (bColour==PDT::Colour3bar) {
if ((cColour==PDT::Colour3bar && aColour==PDT::Colour0)){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
}
else if ((cColour==PDT::Colour0 && aColour==PDT::Colour3bar)){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour3bar){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour3bar && aColour==PDT::Colour8){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
else if(cColour==PDT::Colour3 && aColour==PDT::Colour3) {
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFc ,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(FFa ,ShowerInteraction::QCD));
}
}
// decaying colour octet
else if (bColour==PDT::Colour8){
if ((cColour==PDT::Colour3 && aColour==PDT::Colour3bar) ||
(cColour==PDT::Colour3bar && aColour==PDT::Colour3)){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour0){
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
}
else if (cColour==PDT::Colour0 && aColour==PDT::Colour8){
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
}
}
// decaying colour sextet
else if(bColour==PDT::Colour6) {
if (cColour==PDT::Colour3 && aColour==PDT::Colour3) {
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
}
}
// decaying colour antisextet
else if(bColour==PDT::Colour6bar) {
if (cColour==PDT::Colour3bar && aColour==PDT::Colour3bar) {
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFa,ShowerInteraction::QCD));
dipoles.push_back(DipoleType(IFc,ShowerInteraction::QCD));
}
}
}
// QED dipoles
if(inter==ShowerInteraction::Both ||
inter==ShowerInteraction::QED) {
const bool & bCharged = bProgenitor->dataPtr()->charged();
const bool & cCharged = cProgenitor->dataPtr()->charged();
const bool & aCharged = aProgenitor->dataPtr()->charged();
// initial-final
if(bCharged && aCharged) {
dipoles.push_back(DipoleType(IFba,ShowerInteraction::QED));
dipoles.push_back(DipoleType(IFa ,ShowerInteraction::QED));
}
if(bCharged && cCharged) {
dipoles.push_back(DipoleType(IFbc,ShowerInteraction::QED));
dipoles.push_back(DipoleType(IFc ,ShowerInteraction::QED));
}
// final-state
if(aCharged && cCharged) {
dipoles.push_back(DipoleType(FFa,ShowerInteraction::QED));
dipoles.push_back(DipoleType(FFc,ShowerInteraction::QED));
}
}
// check colour structure is allowed
return !dipoles.empty();
}
vector<Lorentz5Momentum> PerturbativeDecayer::hardMomenta(PPtr in, PPtr emitter,
PPtr spectator,
const vector<DipoleType> &dipoles,
int i, bool inDeadZone) {
// get masses of the particles
mb_ = in ->momentum().mass();
e_ = emitter ->momentum().mass()/mb_;
s_ = spectator->momentum().mass()/mb_;
e2_ = sqr(e_);
s2_ = sqr(s_);
vector<Lorentz5Momentum> particleMomenta;
Energy2 lambda = sqr(mb_)*sqrt(1.+sqr(s2_)+sqr(e2_)-2.*s2_-2.*e2_-2.*s2_*e2_);
// calculate A
double pre = C_;
// multiply by the colour factor of the dipole
// ISR
if (dipoles[i].type==IFba || dipoles[i].type==IFbc) {
pre *= colourCoeff(in->dataPtr(),emitter->dataPtr(),spectator->dataPtr(),dipoles[i]);
}
// radiation from a/c with initial-final connection
else if (dipoles[i].type==IFa || dipoles[i].type==IFc) {
pre *= colourCoeff(emitter->dataPtr(),in->dataPtr(),spectator->dataPtr(),dipoles[i]);
}
// radiation from a/c with final-final connection
else if (dipoles[i].type==FFa || dipoles[i].type==FFc) {
pre *= colourCoeff(emitter->dataPtr(),spectator->dataPtr(),in->dataPtr(),dipoles[i]);
}
double A = 2.*abs(pre)/Constants::twopi;
// factor due sampling choice
if(phaseOpt_==0) A *= ymax_;
// coupling factor
if(dipoles[i].interaction==ShowerInteraction::QCD)
A *= alphaS() ->overestimateValue();
else
A *= alphaEM()->overestimateValue();
Energy pTmax = 0.5*mb_*(1.-sqr(s_+e_));
// if no possible branching return
if ( pTmax < pTmin_ ) return particleMomenta;
// loop over the two regions
for(unsigned int j=0;j<2;++j) {
Energy pT=pTmax;
vector<Lorentz5Momentum> momenta(4);
while (pT >= pTmin_) {
double ymax;
// overestimate with flat y limit
if(phaseOpt_==0) {
pT *= pow(UseRandom::rnd(),(1./A));
ymax=ymax_;
}
// pT sampling including tighter pT dependent y limit
else {
pT = 2.*pTmax*exp(-sqrt(-2.*log(UseRandom::rnd())/A+sqr(log(2.*pTmax/pT))));
// choice of limit overestimate ln(2*pTmax/pT) (true limit acosh(pTmax/pT))
ymax = log(2.*pTmax/pT);
}
if (pT < pTmin_) break;
double phi = UseRandom::rnd()*Constants::twopi;
double y = ymax*(2.*UseRandom::rnd()-1.);
double xs, xe, xe_z, xg;
// check if the momenta are physical
if (!calcMomenta(j, pT, y, phi, xg, xs, xe,
xe_z, momenta))
continue;
// check if point lies within phase space
if (!psCheck(xg, xs)) continue;
// check if point lies within the dead-zone (if required)
if(inDeadZone && !inTotalDeadZone(xg,xs,dipoles,i)) continue;
// decay products for 3 body decay
PPtr inpart = in ->dataPtr()->produceParticle(momenta[0]);
ParticleVector decay3;
decay3.push_back(emitter ->dataPtr()->produceParticle(momenta[1]));
decay3.push_back(spectator->dataPtr()->produceParticle(momenta[2]));
if(dipoles[i].interaction==ShowerInteraction::QCD)
decay3.push_back(getParticleData(ParticleID::g )->produceParticle(momenta[3]));
else
decay3.push_back(getParticleData(ParticleID::gamma)->produceParticle(momenta[3]));
// decay products for 2 body decay
Lorentz5Momentum p1(ZERO,ZERO, lambda/2./mb_,(mb_/2.)*(1.+e2_-s2_),mb_*e_);
Lorentz5Momentum p2(ZERO,ZERO,-lambda/2./mb_,(mb_/2.)*(1.+s2_-e2_),mb_*s_);
ParticleVector decay2;
decay2.push_back(emitter ->dataPtr()->produceParticle(p1));
decay2.push_back(spectator->dataPtr()->produceParticle(p2));
if (dipoles[i].type==FFa || dipoles[i].type==IFa || dipoles[i].type==IFba) {
swap(decay2[0],decay2[1]);
swap(decay3[0],decay3[1]);
}
// calculate matrix element ratio R/B
double meRatio = matrixElementRatio(*inpart,decay2,decay3,Initialize,dipoles[i].interaction);
// calculate dipole factor
double dipoleSum(0.),numerator(0.);
for (int k=0; k<int(dipoles.size()); ++k) {
// skip dipoles which are not of the interaction being considered
if(dipoles[k].interaction!=dipoles[i].interaction) continue;
pair<double,double> dipole = calculateDipole(dipoles[k],*inpart,decay3);
dipoleSum += abs(dipole.first);
if (k==i) numerator = abs(dipole.second);
}
meRatio *= numerator/dipoleSum;
// calculate jacobian
Energy2 denom = (mb_-momenta[3].e())*momenta[2].vect().mag() -
momenta[2].e()*momenta[3].z();
InvEnergy2 J = (momenta[2].vect().mag2())/(lambda*denom);
// calculate weight
double weight = enhance_*meRatio*fabs(sqr(pT)*J)/pre/Constants::twopi;
if(dipoles[i].interaction==ShowerInteraction::QCD)
weight *= alphaS() ->ratio(pT*pT);
else
weight *= alphaEM()->ratio(pT*pT);
// accept point if weight > R
if (pT > pT_ && weight > UseRandom::rnd()) {
particleMomenta=momenta;
if (weight > 1.) {
generator()->log() << "WEIGHT PROBLEM " << fullName() << " " << weight << "\n";
generator()->log() << xe << " " << xs << " " << xg << "\n";
for(unsigned int ix=0;ix<particleMomenta.size();++ix)
generator()->log() << particleMomenta[ix]/GeV << "\n";
}
pT_ = pT;
break;
}
}
}
return particleMomenta;
}
bool PerturbativeDecayer::calcMomenta(int j, Energy pT, double y, double phi,
double& xg, double& xs, double& xe, double& xe_z,
vector<Lorentz5Momentum>& particleMomenta) {
// calculate xg
xg = 2.*pT*cosh(y) / mb_;
if (xg>(1. - sqr(e_ + s_)) || xg<0.) return false;
// calculate the two values of zs
double xT = 2.*pT / mb_;
double zg = 2.*pT*sinh(y) / mb_;
double A = (sqr(xT) - 4. * xg + 4.);
double B = 2. * zg * (s2_ - e2_ - xg + 1.);
double det = -4. * (-sqr(s2_) + (2. * e2_ + sqr(xT) - 2. * xg + 2.) * s2_ - sqr(e2_ + xg - 1.)) * sqr(xg - 2.);
if (det<0.) return false;
double zs= j==0 ? (-B+sqrt(det))/A : (-B-sqrt(det))/A;
// zs must be negative
if(zs>0.) return false;
xs = sqrt(sqr(zs)+4.*s2_);
// check value of xs is physical
if (xs>(1.+s2_-e2_) || xs<2.*s_) return false;
// calculate xe
xe = 2.-xs-xg;
// check value of xe is physical
if (xe>(1.+e2_-s2_) || xe<2.*e_) return false;
// calculate xe_z
xe_z = -zg-zs;
// calculate 4 momenta
particleMomenta[0].setE ( mb_);
particleMomenta[0].setX ( ZERO);
particleMomenta[0].setY ( ZERO);
particleMomenta[0].setZ ( ZERO);
particleMomenta[0].setMass( mb_);
particleMomenta[1].setE ( mb_*xe/2.);
particleMomenta[1].setX (-pT*cos(phi));
particleMomenta[1].setY (-pT*sin(phi));
particleMomenta[1].setZ ( mb_*xe_z/2.);
particleMomenta[1].setMass( mb_*e_);
particleMomenta[2].setE ( mb_*xs/2.);
particleMomenta[2].setX ( ZERO);
particleMomenta[2].setY ( ZERO);
particleMomenta[2].setZ ( mb_*zs/2.);
particleMomenta[2].setMass( mb_*s_);
particleMomenta[3].setE ( pT*cosh(y));
particleMomenta[3].setX ( pT*cos(phi));
particleMomenta[3].setY ( pT*sin(phi));
particleMomenta[3].setZ ( pT*sinh(y));
particleMomenta[3].setMass( ZERO);
return true;
}
bool PerturbativeDecayer::psCheck(const double xg, const double xs) {
// check is point is in allowed region of phase space
double xe_star = (1.-s2_+e2_-xg)/sqrt(1.-xg);
double xg_star = xg/sqrt(1.-xg);
if ((sqr(xe_star)-4.*e2_) < 1e-10) return false;
double xs_max = (4.+4.*s2_-sqr(xe_star+xg_star)+
sqr(sqrt(sqr(xe_star)-4.*e2_)+xg_star))/ 4.;
double xs_min = (4.+4.*s2_-sqr(xe_star+xg_star)+
sqr(sqrt(sqr(xe_star)-4.*e2_)-xg_star))/ 4.;
if (xs < xs_min || xs > xs_max) return false;
return true;
}
pair<double,double> PerturbativeDecayer::calculateDipole(const DipoleType & dipoleId,
const Particle & inpart,
const ParticleVector & decay3) {
// calculate dipole for decay b->ac
pair<double,double> dipole = make_pair(0.,0.);
double x1 = 2.*decay3[0]->momentum().e()/mb_;
double x2 = 2.*decay3[1]->momentum().e()/mb_;
double xg = 2.*decay3[2]->momentum().e()/mb_;
double mu12 = sqr(decay3[0]->mass()/mb_);
double mu22 = sqr(decay3[1]->mass()/mb_);
tcPDPtr part[3] = {inpart.dataPtr(),decay3[0]->dataPtr(),decay3[1]->dataPtr()};
if(dipoleId.type==FFa || dipoleId.type == IFa || dipoleId.type == IFba) {
swap(part[1],part[2]);
swap(x1,x2);
swap(mu12,mu22);
}
// radiation from b with initial-final connection
if (dipoleId.type==IFba || dipoleId.type==IFbc) {
dipole.first = -2./sqr(xg);
dipole.first *= colourCoeff(part[0],part[1],part[2],dipoleId);
}
// radiation from a/c with initial-final connection
else if (dipoleId.type==IFa || dipoleId.type==IFc) {
double z = 1. - xg/(1.-mu22+mu12);
dipole.first = (-2.*mu12/sqr(1.-x2+mu22-mu12) + (1./(1.-x2+mu22-mu12))*
(2./(1.-z)-dipoleSpinFactor(part[1],z)));
dipole.first *= colourCoeff(part[1],part[0],part[2],dipoleId);
}
// radiation from a/c with final-final connection
else if (dipoleId.type==FFa || dipoleId.type==FFc) {
double z = 1. + ((x1-1.+mu22-mu12)/(x2-2.*mu22));
double y = (1.-x2-mu12+mu22)/(1.-mu12-mu22);
double vt = sqrt((1.-sqr(e_+s_))*(1.-sqr(e_-s_)))/(1.-mu12-mu22);
double v = sqrt(sqr(2.*mu22+(1.-mu12-mu22)*(1.-y))-4.*mu22)
/(1.-y)/(1.-mu12-mu22);
if(part[1]->iSpin()!=PDT::Spin1) {
dipole.first = (1./(1.-x2+mu22-mu12))*
((2./(1.-z*(1.-y)))-vt/v*(dipoleSpinFactor(part[1],z)+(2.*mu12/(1.+mu22-mu12-x2))));
}
else {
dipole.first = (1./(1.-x2+mu22-mu12))*
(1./(1.-z*(1.-y))+1./(1.-(1.-z)*(1.-y))+(z*(1.-z)-2.)/v-vt/v*(2.*mu12/(1.+mu22-mu12-x2)));
dipole.second = (1./(1.-x2+mu22-mu12))*
(2./(1.-z*(1.-y))+(z*(1.-z)-2.)/v-vt/v*(2.*mu12/(1.+mu22-mu12-x2)));
dipole.second *= colourCoeff(part[1],part[2],part[0],dipoleId);
}
dipole.first *= colourCoeff(part[1],part[2],part[0],dipoleId);
}
// special for the case that all particles are gluons
else if(dipoleId.type==FFg) {
double z = (1.-x2)/xg;
double y = 1.-xg;
dipole.first = 1./(1.-xg)*(1./(1.-z*(1.-y))+1./(1.-(1.-z)*(1.-y))+(z*(1.-z)-2.));
dipole.first *= colourCoeff(part[1],part[2],part[0],dipoleId);
}
else
assert(false);
// coupling prefactors
if(dipole.second==0.) dipole.second=dipole.first;
dipole.first *= 8.*Constants::pi;
dipole.second *= 8.*Constants::pi;
// return the answer
return dipole;
}
double PerturbativeDecayer::dipoleSpinFactor(tcPDPtr part, double z){
// calculate the spin dependent component of the dipole
if (part->iSpin()==PDT::Spin0)
return 2.;
else if (part->iSpin()==PDT::Spin1Half)
return (1. + z);
else if (part->iSpin()==PDT::Spin1)
return -(z*(1.-z) - 1./(1.-z) + 1./z -2.);
return 0.;
}
namespace {
-double colourCharge(PDT::Colour icol) {
+double colourChargeOld(PDT::Colour icol) {
switch(icol) {
case PDT::Colour0 :
return 0.;
case PDT::Colour3 : case PDT::Colour3bar :
return 4./3.;
case PDT::Colour8:
return 3.;
case PDT::Colour6 : case PDT::Colour6bar :
return 10./3.;
default :
assert(false);
}
}
}
double PerturbativeDecayer::colourCoeff(tcPDPtr emitter,
tcPDPtr spectator,
tcPDPtr other,
DipoleType dipole) {
if(dipole.interaction==ShowerInteraction::QCD) {
- double emitterColour = colourCharge(emitter ->iColour());
- double spectatorColour = colourCharge(spectator->iColour());
- double otherColour = colourCharge(other ->iColour());
+ double emitterColour = colourChargeOld(emitter ->iColour());
+ double spectatorColour = colourChargeOld(spectator->iColour());
+ double otherColour = colourChargeOld(other ->iColour());
double val = 0.5*(sqr(emitterColour)+sqr(spectatorColour)-sqr(otherColour))/emitterColour;
return val;
}
else {
double val = double(emitter->iCharge()*spectator->iCharge())/9.;
// FF dipoles
if(dipole.type==FFa || dipole.type == FFc) return -val;
// IF dipoles
else return val;
}
}
void PerturbativeDecayer::getColourLines(RealEmissionProcessPtr real) {
// extract the particles
vector<tPPtr> branchingPart;
branchingPart.push_back(real->incoming()[0]);
for(unsigned int ix=0;ix<real->outgoing().size();++ix) {
branchingPart.push_back(real->outgoing()[ix]);
}
vector<unsigned int> sing,trip,atrip,oct,sex,asex;
for (size_t ib=0;ib<branchingPart.size()-1;++ib) {
if (branchingPart[ib]->dataPtr()->iColour()==PDT::Colour0 ) sing. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour3 ) trip. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour3bar) atrip.push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour8 ) oct. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour6 ) sex. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour6bar) asex. push_back(ib);
}
// decaying colour singlet
if (branchingPart[0]->dataPtr()->iColour()==PDT::Colour0) {
// 0 -> 3 3bar
if (trip.size()==1 && atrip.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[atrip[0]]->colourConnect(branchingPart[ 3 ]);
branchingPart[ 3 ]->colourConnect(branchingPart[trip[0]]);
}
else {
branchingPart[atrip[0]]->colourConnect(branchingPart[trip[0]]);
}
}
// 0 -> 8 8
else if (oct.size()==2 ) {
if(real->interaction()==ShowerInteraction::QCD) {
bool col = UseRandom::rndbool();
branchingPart[oct[0]]->colourConnect(branchingPart[ 3 ],col);
branchingPart[ 3 ]->colourConnect(branchingPart[oct[1]],col);
branchingPart[oct[1]]->colourConnect(branchingPart[oct[0]],col);
}
else {
branchingPart[oct[0]]->colourConnect(branchingPart[oct[1]]);
branchingPart[oct[1]]->colourConnect(branchingPart[oct[0]]);
}
}
else
assert(real->interaction()==ShowerInteraction::QED);
}
// decaying colour triplet
else if (branchingPart[0]->dataPtr()->iColour()==PDT::Colour3 ) {
// 3 -> 3 0
if (trip.size()==2 && sing.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[3]->incomingColour(branchingPart[trip[0]]);
branchingPart[3]-> colourConnect(branchingPart[trip[1]]);
}
else {
branchingPart[trip[1]]->incomingColour(branchingPart[trip[0]]);
}
}
// 3 -> 3 8
else if (trip.size()==2 && oct.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
// 8 emit incoming partner
if(real->emitter()==oct[0]&&real->spectator()==0) {
branchingPart[ 3 ]->incomingColour(branchingPart[trip[0]]);
branchingPart[ 3 ]-> colourConnect(branchingPart[oct[0] ]);
branchingPart[oct[0]]-> colourConnect(branchingPart[trip[1]]);
}
// 8 emit final spectator or vice veras
else {
branchingPart[oct[0]]->incomingColour(branchingPart[trip[0]]);
branchingPart[oct[0]]-> colourConnect(branchingPart[ 3 ]);
branchingPart[ 3 ]-> colourConnect(branchingPart[trip[1]]);
}
}
else {
branchingPart[oct[0]]->incomingColour(branchingPart[trip[0]]);
branchingPart[oct[0]]-> colourConnect(branchingPart[trip[1]]);
}
}
// 3 -> 3bar 3bar
else if(trip.size() ==1 && atrip.size()==2) {
if(real->interaction()==ShowerInteraction::QCD) {
if(real->emitter()==atrip[0]) {
branchingPart[3]->colourConnect(branchingPart[atrip[0]],true);
tColinePtr col[3] = {ColourLine::create(branchingPart[ trip[0]],false),
ColourLine::create(branchingPart[ 3],true ),
ColourLine::create(branchingPart[atrip[1]],true)};
col[0]->setSinkNeighbours(col[1],col[2]);
}
else {
branchingPart[3]->colourConnect(branchingPart[atrip[1]],true);
tColinePtr col[3] = {ColourLine::create(branchingPart[ trip[0]],false),
ColourLine::create(branchingPart[atrip[0]],true ),
ColourLine::create(branchingPart[ 3],true)};
col[0]->setSinkNeighbours(col[1],col[2]);
}
}
else {
tColinePtr col[3] = {ColourLine::create(branchingPart[ trip[0]],false),
ColourLine::create(branchingPart[atrip[0]],true ),
ColourLine::create(branchingPart[atrip[1]],true)};
col[0]->setSinkNeighbours(col[1],col[2]);
}
}
else
assert(false);
}
// decaying colour anti-triplet
else if (branchingPart[0]->dataPtr()->iColour()==PDT::Colour3bar) {
// 3bar -> 3bar 0
if (atrip.size()==2 && sing.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[3]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[3]-> colourConnect(branchingPart[atrip[1]],true);
}
else {
branchingPart[atrip[1]]->incomingColour(branchingPart[atrip[0]],true);
}
}
// 3 -> 3 8
else if (atrip.size()==2 && oct.size()==1){
if(real->interaction()==ShowerInteraction::QCD) {
// 8 emit incoming partner
if(real->emitter()==oct[0]&&real->spectator()==0) {
branchingPart[ 3 ]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[ 3 ]-> colourConnect(branchingPart[oct[0] ],true);
branchingPart[oct[0]]-> colourConnect(branchingPart[atrip[1]],true);
}
// 8 emit final spectator or vice veras
else {
if(real->interaction()==ShowerInteraction::QCD) {
branchingPart[oct[0]]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[oct[0]]-> colourConnect(branchingPart[ 3 ],true);
branchingPart[3]-> colourConnect(branchingPart[atrip[1]] ,true);
}
}
}
else {
branchingPart[oct[0]]->incomingColour(branchingPart[atrip[0]],true);
branchingPart[oct[0]]-> colourConnect(branchingPart[atrip[1]],true);
}
}
// 3bar -> 3 3
else if(atrip.size() ==1 && trip.size()==2) {
if(real->interaction()==ShowerInteraction::QCD) {
if(real->emitter()==trip[0]) {
branchingPart[3]->colourConnect(branchingPart[trip[0]],false);
tColinePtr col[3] = {ColourLine::create(branchingPart[atrip[0]],true ),
ColourLine::create(branchingPart[ 3],false),
ColourLine::create(branchingPart[ trip[1]],false)};
col[0]->setSourceNeighbours(col[1],col[2]);
}
else {
branchingPart[3]->colourConnect(branchingPart[trip[1]],false);
tColinePtr col[3] = {ColourLine::create(branchingPart[atrip[0]],true ),
ColourLine::create(branchingPart[ trip[0]],false),
ColourLine::create(branchingPart[ 3],false)};
col[0]->setSourceNeighbours(col[1],col[2]);
}
}
else {
tColinePtr col[3] = {ColourLine::create(branchingPart[atrip[0]],true ),
ColourLine::create(branchingPart[ trip[0]],false),
ColourLine::create(branchingPart[ trip[1]],false)};
col[0]->setSourceNeighbours(col[1],col[2]);
}
}
else
assert(false);
}
// decaying colour octet
else if(branchingPart[0]->dataPtr()->iColour()==PDT::Colour8 ) {
// 8 -> 3 3bar
if (trip.size()==1 && atrip.size()==1) {
if(real->interaction()==ShowerInteraction::QCD) {
// 3 emits
if(trip[0]==real->emitter()) {
branchingPart[3] ->incomingColour(branchingPart[oct[0]] );
branchingPart[3] -> colourConnect(branchingPart[trip[0]]);
branchingPart[atrip[0]]->incomingColour(branchingPart[oct[0]],true);
}
// 3bar emits
else {
branchingPart[3] ->incomingColour(branchingPart[oct[0]] ,true);
branchingPart[3] -> colourConnect(branchingPart[atrip[0]],true);
branchingPart[trip[0]]->incomingColour(branchingPart[oct[0]] );
}
}
else {
branchingPart[trip[0]]->incomingColour(branchingPart[oct[0]] );
branchingPart[atrip[0]]->incomingColour(branchingPart[oct[0]],true);
}
}
// 8 -> 8 0
else if (sing.size()==1 && oct.size()==2) {
if(real->interaction()==ShowerInteraction::QCD) {
bool col = UseRandom::rndbool();
branchingPart[ 3 ]->colourConnect (branchingPart[oct[1]], col);
branchingPart[ 3 ]->incomingColour(branchingPart[oct[0]], col);
branchingPart[oct[1]]->incomingColour(branchingPart[oct[0]],!col);
}
else {
branchingPart[oct[1]]->incomingColour(branchingPart[oct[0]]);
branchingPart[oct[1]]->incomingColour(branchingPart[oct[0]],true);
}
}
else
assert(false);
}
// sextet
else if(branchingPart[0]->dataPtr()->iColour() == PDT::Colour6) {
if(trip.size()==2) {
if(real->interaction()==ShowerInteraction::QCD) {
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(branchingPart[0]->colourInfo());
if(trip[0]==real->emitter()) {
ColinePtr cline = new_ptr(ColourLine());
parentColour->colourLine(cline);
cline->addColoured(branchingPart[3]);
branchingPart[3] -> colourConnect(branchingPart[trip[0]]);
cline = new_ptr(ColourLine());
parentColour->colourLine(cline);
cline->addColoured(branchingPart[trip[1]]);
}
else {
ColinePtr cline = new_ptr(ColourLine());
parentColour->colourLine(cline);
cline->addColoured(branchingPart[3]);
branchingPart[3] -> colourConnect(branchingPart[trip[1]]);
cline = new_ptr(ColourLine());
parentColour->colourLine(cline);
cline->addColoured(branchingPart[trip[0]]);
}
}
else {
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(branchingPart[0]->colourInfo());
for(unsigned int ix=0;ix<2;++ix) {
ColinePtr cline = new_ptr(ColourLine());
parentColour->colourLine(cline);
cline->addColoured(branchingPart[trip[ix]]);
}
}
}
else
assert(false);
}
// antisextet
else if(branchingPart[0]->dataPtr()->iColour() == PDT::Colour6bar) {
if(atrip.size()==2) {
if(real->interaction()==ShowerInteraction::QCD) {
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(branchingPart[0]->colourInfo());
if(atrip[0]==real->emitter()) {
ColinePtr cline = new_ptr(ColourLine());
parentColour->antiColourLine(cline);
cline->addAntiColoured(branchingPart[3]);
branchingPart[3]->antiColourConnect(branchingPart[atrip[0]]);
cline = new_ptr(ColourLine());
parentColour->antiColourLine(cline);
cline->addAntiColoured(branchingPart[atrip[1]]);
}
else {
ColinePtr cline = new_ptr(ColourLine());
parentColour->antiColourLine(cline);
cline->addAntiColoured(branchingPart[3]);
branchingPart[3]->antiColourConnect(branchingPart[atrip[1]]);
cline = new_ptr(ColourLine());
parentColour->antiColourLine(cline);
cline->addAntiColoured(branchingPart[trip[0]]);
}
}
else {
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(branchingPart[0]->colourInfo());
for(unsigned int ix=0;ix<2;++ix) {
ColinePtr cline = new_ptr(ColourLine());
parentColour->antiColourLine(cline);
cline->addColoured(branchingPart[atrip[ix]],true);
}
}
}
else
assert(false);
}
else
assert(false);
}
PerturbativeDecayer::phaseSpaceRegion
PerturbativeDecayer::inInitialFinalDeadZone(double xg, double xa,
double a, double c) const {
double lam = sqrt(1.+a*a+c*c-2.*a-2.*c-2.*a*c);
double kappab = 1.+0.5*(1.-a+c+lam);
double kappac = kappab-1.+c;
double kappa(0.);
// check whether or not in the region for emission from c
double r = 0.5;
if(c!=0.) r += 0.5*c/(1.+a-xa);
double pa = sqrt(sqr(xa)-4.*a);
double z = ((2.-xa)*(1.-r)+r*pa-xg)/pa;
if(z<1. && z>0.) {
kappa = (1.+a-c-xa)/(z*(1.-z));
if(kappa<kappac)
return emissionFromC;
}
// check in region for emission from b (T1)
double cq = sqr(1.+a-c)-4*a;
double bq = -2.*kappab*(1.-a-c);
double aq = sqr(kappab)-4.*a*(kappab-1);
double dis = sqr(bq)-4.*aq*cq;
z=1.-(-bq-sqrt(dis))/2./aq;
double w = 1.-(1.-z)*(kappab-1.);
double xgmax = (1.-z)*kappab;
// possibly in T1 region
if(xg<xgmax) {
z = 1.-xg/kappab;
kappa=kappab;
}
// possibly in T2 region
else {
aq = 4.*a;
bq = -4.*a*(2.-xg);
cq = sqr(1.+a-c-xg);
dis = sqr(bq)-4.*aq*cq;
z = (-bq-sqrt(dis))/2./aq;
kappa = xg/(1.-z);
}
// compute limit on xa
double u = 1.+a-c-(1.-z)*kappa;
w = 1.-(1.-z)*(kappa-1.);
double v = sqr(u)-4.*z*a*w;
if(v<0. && v>-1e-10) v= 0.;
v = sqrt(v);
if(xa<0.5*((u+v)/w+(u-v)/z)) {
if(xg<xgmax)
return emissionFromA1;
else if(useMEforT2_)
return deadZone;
else
return emissionFromA2;
}
else
return deadZone;
}
PerturbativeDecayer::phaseSpaceRegion
PerturbativeDecayer::inFinalFinalDeadZone(double xb, double xc,
double b, double c) const {
// basic kinematics
double lam = sqrt(1.+b*b+c*c-2.*b-2.*c-2.*b*c);
// check whether or not in the region for emission from b
double r = 0.5;
if(b!=0.) r+=0.5*b/(1.+c-xc);
double pc = sqrt(sqr(xc)-4.*c);
double z = -((2.-xc)*r-r*pc-xb)/pc;
if(z<1. and z>0.) {
if((1.-b+c-xc)/(z*(1.-z))<0.5*(1.+b-c+lam)) return emissionFromB;
}
// check whether or not in the region for emission from c
r = 0.5;
if(c!=0.) r+=0.5*c/(1.+b-xb);
double pb = sqrt(sqr(xb)-4.*b);
z = -((2.-xb)*r-r*pb-xc)/pb;
if(z<1. and z>0.) {
if((1.-c+b-xb)/(z*(1.-z))<0.5*(1.-b+c+lam)) return emissionFromC;
}
return deadZone;
}
bool PerturbativeDecayer::inTotalDeadZone(double xg, double xs,
const vector<DipoleType> & dipoles,
int i) {
double xb,xc,b,c;
if(dipoles[i].type==FFa || dipoles[i].type == IFa || dipoles[i].type == IFba) {
xc = xs;
xb = 2.-xg-xs;
b = e2_;
c = s2_;
}
else {
xb = xs;
xc = 2.-xg-xs;
b = s2_;
c = e2_;
}
for(unsigned int ix=0;ix<dipoles.size();++ix) {
if(dipoles[ix].interaction!=dipoles[i].interaction)
continue;
// should also remove negative QED dipoles but shouldn't be an issue unless we
// support QED ME corrections
switch (dipoles[ix].type) {
case FFa :
if(inFinalFinalDeadZone(xb,xc,b,c)!=deadZone) return false;
break;
case FFc :
if(inFinalFinalDeadZone(xc,xb,c,b)!=deadZone) return false;
break;
case IFa : case IFba:
if(inInitialFinalDeadZone(xg,xc,c,b)!=deadZone) return false;
break;
case IFc : case IFbc:
if(inInitialFinalDeadZone(xg,xb,b,c)!=deadZone) return false;
break;
case FFg:
break;
}
}
return true;
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 10:03 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242584
Default Alt Text
(99 KB)

Event Timeline