Page MenuHomeHEPForge
Files F8308627

No OneTemporary
Actions

View Options

diff --git a/Models/General/TwoBodyDecayConstructor.cc b/Models/General/TwoBodyDecayConstructor.cc
--- a/Models/General/TwoBodyDecayConstructor.cc
+++ b/Models/General/TwoBodyDecayConstructor.cc
@@ -1,432 +1,446 @@
// -*- C++ -*-
//
// TwoBodyDecayConstructor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TwoBodyDecayConstructor class.
//
#include "TwoBodyDecayConstructor.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "Herwig/Decay/General/GeneralTwoBodyDecayer.h"
#include "Herwig/Models/StandardModel/StandardModel.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "DecayConstructor.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractSSTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractRFSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractRFVVertex.fh"
+#include "VectorCurrentDecayConstructor.h"
using namespace Herwig;
using ThePEG::Helicity::VertexBasePtr;
IBPtr TwoBodyDecayConstructor::clone() const {
return new_ptr(*this);
}
IBPtr TwoBodyDecayConstructor::fullclone() const {
return new_ptr(*this);
}
void TwoBodyDecayConstructor::persistentOutput(PersistentOStream & os) const {
os << alphaQCD_ << alphaQED_ << oenum(inter_);
}
void TwoBodyDecayConstructor::persistentInput(PersistentIStream & is, int) {
is >> alphaQCD_ >> alphaQED_>> ienum(inter_);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<TwoBodyDecayConstructor,NBodyDecayConstructorBase>
describeHerwigTwoBodyDecayConstructor("Herwig::TwoBodyDecayConstructor", "Herwig.so");
void TwoBodyDecayConstructor::Init() {
static ClassDocumentation<TwoBodyDecayConstructor> documentation
("The TwoBodyDecayConstructor implements to creation of 2 body decaymodes "
"and decayers that do not already exist for the given set of vertices.");
static Reference<TwoBodyDecayConstructor,ShowerAlpha> interfaceShowerAlphaQCD
("AlphaQCD",
"The coupling for QCD corrections",
&TwoBodyDecayConstructor::alphaQCD_, false, false, true, false, false);
static Reference<TwoBodyDecayConstructor,ShowerAlpha> interfaceShowerAlphaQED
("AlphaQED",
"The coupling for QED corrections",
&TwoBodyDecayConstructor::alphaQED_, false, false, true, false, false);
static Switch<TwoBodyDecayConstructor,ShowerInteraction> interfaceInteractions
("Interactions",
"which interactions to include for the hard corrections",
&TwoBodyDecayConstructor::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);
}
void TwoBodyDecayConstructor::DecayList(const set<PDPtr> & particles) {
+ // special for weak decays
+ for(unsigned int ix=0;ix<decayConstructor()->decayConstructors().size();++ix) {
+ Ptr<Herwig::VectorCurrentDecayConstructor>::pointer
+ weak = dynamic_ptr_cast<Ptr<Herwig::VectorCurrentDecayConstructor>::pointer >
+ (decayConstructor()->decayConstructors()[ix]);
+ if(!weak) continue;
+ weakMassCut_ = max(weakMassCut_,weak->massCut());
+ }
if( particles.empty() ) return;
tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
unsigned int nv(model->numberOfVertices());
for(set<PDPtr>::const_iterator ip=particles.begin();
ip!=particles.end();++ip) {
tPDPtr parent = *ip;
if ( Debug::level > 0 )
Repository::cout() << "Constructing 2-body decays for "
<< parent->PDGName() << '\n';
multiset<TwoBodyDecay> decays;
for(unsigned int iv = 0; iv < nv; ++iv) {
if(excluded(model->vertex(iv)) ||
model->vertex(iv)->getNpoint()>3) continue;
for(unsigned int il = 0; il < 3; ++il)
createModes(parent, model->vertex(iv), il,decays);
}
if( !decays.empty() ) createDecayMode(decays);
}
}
void TwoBodyDecayConstructor::
createModes(tPDPtr inpart, VertexBasePtr vertex,
unsigned int list, multiset<TwoBodyDecay> & modes) {
if( !vertex->isIncoming(inpart) || vertex->getNpoint() != 3 )
return;
Energy m1(inpart->mass());
tPDPtr ccpart = inpart->CC() ? inpart->CC() : inpart;
long id = ccpart->id();
tPDVector decaylist = vertex->search(list, ccpart);
tPDVector::size_type nd = decaylist.size();
for( tPDVector::size_type i = 0; i < nd; i += 3 ) {
tPDPtr pa(decaylist[i]), pb(decaylist[i + 1]), pc(decaylist[i + 2]);
if( pb->id() == id ) swap(pa, pb);
if( pc->id() == id ) swap(pa, pc);
//allowed on-shell decay?
if( m1 <= pb->mass() + pc->mass() ) continue;
+ // double counting with current decays?
+ if(inpart->iSpin()==PDT::Spin1 && inpart->iCharge()==0 &&
+ pb->id()==-pc->id() && abs(pb->id())<=3 && inpart->mass() <= weakMassCut_ ) {
+ continue;
+ }
//vertices are defined with all particles incoming
modes.insert( TwoBodyDecay(inpart,pb, pc, vertex) );
}
}
GeneralTwoBodyDecayerPtr
TwoBodyDecayConstructor::createDecayer(TwoBodyDecay decay,
vector<VertexBasePtr> vertices) {
string name;
using namespace Helicity::VertexType;
PDT::Spin in = decay.parent_->iSpin();
// PDT::Spin out1 = decay.children_.first ->iSpin();
PDT::Spin out2 = decay.children_.second->iSpin();
switch(decay.vertex_->getName()) {
case FFV :
if(in == PDT::Spin1Half) {
name = "FFVDecayer";
if(out2==PDT::Spin1Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "VFFDecayer";
}
break;
case FFS :
if(in == PDT::Spin1Half) {
name = "FFSDecayer";
if(out2==PDT::Spin1Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "SFFDecayer";
}
break;
case VVS :
if(in == PDT::Spin1) {
name = "VVSDecayer";
if(out2==PDT::Spin1)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "SVVDecayer";
}
break;
case VSS :
if(in == PDT::Spin1) {
name = "VSSDecayer";
}
else {
name = "SSVDecayer";
if(out2==PDT::Spin0)
swap(decay.children_.first,decay.children_.second);
}
break;
case VVT :
name = in==PDT::Spin2 ? "TVVDecayer" : "Unknown";
break;
case FFT :
name = in==PDT::Spin2 ? "TFFDecayer" : "Unknown";
break;
case SST :
name = in==PDT::Spin2 ? "TSSDecayer" : "Unknown";
break;
case SSS :
name = "SSSDecayer";
break;
case VVV :
name = "VVVDecayer";
break;
case RFS :
if(in==PDT::Spin1Half) {
name = "FRSDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else if(in==PDT::Spin0) {
name = "SRFDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "Unknown";
}
break;
case RFV :
if(in==PDT::Spin1Half) {
name = "FRVDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else
name = "Unknown";
break;
default : Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName() << Exception::runerror;
}
if(name=="Unknown")
Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName() << Exception::runerror;
ostringstream fullname;
fullname << "/Herwig/Decays/" << name << "_" << decay.parent_->PDGName()
<< "_" << decay.children_.first ->PDGName()
<< "_" << decay.children_.second->PDGName();
string classname = "Herwig::" + name;
GeneralTwoBodyDecayerPtr decayer;
decayer = dynamic_ptr_cast<GeneralTwoBodyDecayerPtr>
(generator()->preinitCreate(classname,fullname.str()));
if(!decayer)
Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName() << Exception::runerror;
// set the strong coupling for radiation
generator()->preinitInterface(decayer, "AlphaS" , "set", alphaQCD_->fullName());
// set the EM coupling for radiation
generator()->preinitInterface(decayer, "AlphaEM", "set", alphaQED_->fullName());
// set the type of interactions for the correction
if(inter_==ShowerInteraction::QCD)
generator()->preinitInterface(decayer, "Interactions", "set", "QCD");
else if(inter_==ShowerInteraction::QED)
generator()->preinitInterface(decayer, "Interactions", "set", "QED");
else
generator()->preinitInterface(decayer, "Interactions", "set", "QCDandQED");
// get the vertices for radiation from the external legs
map<ShowerInteraction,VertexBasePtr> inRad,fourRad;
vector<map<ShowerInteraction,VertexBasePtr> > outRad(2);
vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
for(auto & inter : itemp) {
inRad[inter] = radiationVertex(decay.parent_,inter);
outRad[0][inter] = radiationVertex(decay.children_.first ,inter);
outRad[1][inter] = radiationVertex(decay.children_.second,inter);
// get any contributing 4 point vertices
fourRad[inter] = radiationVertex(decay.parent_,inter, decay.children_);
}
// set info on decay
decayer->setDecayInfo(decay.parent_,decay.children_,vertices,
inRad,outRad,fourRad);
// initialised the decayer
setDecayerInterfaces(fullname.str());
decayer->init();
return decayer;
}
void TwoBodyDecayConstructor::
createDecayMode(multiset<TwoBodyDecay> & decays) {
tPDPtr inpart = decays.begin()->parent_;
for( multiset<TwoBodyDecay>::iterator dit = decays.begin();
dit != decays.end(); ) {
TwoBodyDecay mode = *dit;
// get all the moees with the same in and outgoing particles
pair<multiset<TwoBodyDecay>::iterator,
multiset<TwoBodyDecay>::iterator> range = decays.equal_range(mode);
// construct the decay mode
tPDPtr pb((mode).children_.first), pc((mode).children_.second);
string tag = inpart->name() + "->" + pb->name() + "," +
pc->name() + ";";
// Does it exist already ?
tDMPtr dm = generator()->findDecayMode(tag);
// find the vertices
vector<VertexBasePtr> vertices;
for ( multiset<TwoBodyDecay>::iterator dit2 = range.first;
dit2 != range.second; ++dit2) {
vertices.push_back(dit2->vertex_);
}
dit=range.second;
// now create DecayMode objects that do not already exist
if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) {
tDMPtr ndm = generator()->preinitCreateDecayMode(tag);
if(ndm) {
inpart->stable(false);
GeneralTwoBodyDecayerPtr decayer=createDecayer(mode,vertices);
if(!decayer) continue;
generator()->preinitInterface(ndm, "Decayer", "set",
decayer->fullName());
generator()->preinitInterface(ndm, "Active", "set", "Yes");
Energy width =
decayer->partialWidth(make_pair(inpart,inpart->mass()),
make_pair(pb,pb->mass()) ,
make_pair(pc,pc->mass()));
setBranchingRatio(ndm, width);
if(width==ZERO || ndm->brat()<decayConstructor()->minimumBR()) {
generator()->preinitInterface(decayer->fullName(),
"Initialize", "set","0");
}
}
else
Throw<NBodyDecayConstructorError>()
<< "TwoBodyDecayConstructor::createDecayMode - Needed to create "
<< "new decaymode but one could not be created for the tag "
<< tag << Exception::warning;
}
else if( dm ) {
if(dm->brat()<decayConstructor()->minimumBR()) {
continue;
}
if((dm->decayer()->fullName()).find("Mambo") != string::npos) {
inpart->stable(false);
GeneralTwoBodyDecayerPtr decayer=createDecayer(mode,vertices);
if(!decayer) continue;
generator()->preinitInterface(dm, "Decayer", "set",
decayer->fullName());
Energy width =
decayer->partialWidth(make_pair(inpart,inpart->mass()),
make_pair(pb,pb->mass()) ,
make_pair(pc,pc->mass()));
if(width/(dm->brat()*inpart->width())<1e-10) {
string message = "Herwig calculation of the partial width for the decay mode "
+ inpart->PDGName() + " -> " + pb->PDGName() + " " + pc->PDGName()
+ " is zero.\n This will cause problems with the calculation of"
+ " spin correlations.\n It is probably due to inconsistent parameters"
+ " and decay modes being passed to Herwig via the SLHA file.\n"
+ " Zeroing the branching ratio for this mode.";
setBranchingRatio(dm,ZERO);
generator()->logWarning(NBodyDecayConstructorError(message,Exception::warning));
}
}
}
}
// update CC mode if it exists
if( inpart->CC() ) inpart->CC()->synchronize();
}
VertexBasePtr TwoBodyDecayConstructor::radiationVertex(tPDPtr particle,
ShowerInteraction inter,
tPDPair children) {
tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
map<tPDPtr,VertexBasePtr>::iterator rit = radiationVertices_[inter].find(particle);
tPDPtr cc = particle->CC() ? particle->CC() : particle;
if(children==tPDPair() && rit!=radiationVertices_[inter].end()) return rit->second;
unsigned int nv(model->numberOfVertices());
long bosonID = inter==ShowerInteraction::QCD ? ParticleID::g : ParticleID::gamma;
tPDPtr gluon = getParticleData(bosonID);
// look for radiation vertices for incoming and outgoing particles
for(unsigned int iv=0;iv<nv;++iv) {
VertexBasePtr vertex = model->vertex(iv);
// look for 3 point vertices
if (children==tPDPair()){
if( !vertex->isIncoming(particle) || vertex->getNpoint() != 3 ||
!vertex->isOutgoing(particle) || !vertex->isOutgoing(gluon)) continue;
for(unsigned int list=0;list<3;++list) {
tPDVector decaylist = vertex->search(list, particle);
for( tPDVector::size_type i = 0; i < decaylist.size(); i += 3 ) {
tPDPtr pa(decaylist[i]), pb(decaylist[i + 1]), pc(decaylist[i + 2]);
if( pb->id() == bosonID ) swap(pa, pb);
if( pc->id() == bosonID ) swap(pa, pc);
if( pb->id() != particle->id()) swap(pb, pc);
if( pa->id() != bosonID) continue;
if( pb != particle) continue;
if( pc != cc) continue;
radiationVertices_[inter][particle] = vertex;
return vertex;
}
}
}
// look for 4 point vertex including a gluon
else {
if( !vertex->isIncoming(particle) || vertex->getNpoint()!=4 ||
!vertex->isOutgoing(children.first) || !vertex->isOutgoing(children.second) ||
!vertex->isOutgoing(gluon)) continue;
for(unsigned int list=0;list<4;++list) {
tPDVector decaylist = vertex->search(list, particle);
for( tPDVector::size_type i = 0; i < decaylist.size(); i += 4 ) {
tPDPtr pa(decaylist[i]), pb(decaylist[i+1]), pc(decaylist[i+2]), pd(decaylist[i+3]);
// order so that a = g, b = parent
if( pb->id() == bosonID ) swap(pa, pb);
if( pc->id() == bosonID ) swap(pa, pc);
if( pd->id() == bosonID ) swap(pa, pd);
if( pc->id() == particle->id()) swap(pb, pc);
if( pd->id() == particle->id()) swap(pb, pd);
if( pa->id() != bosonID) continue;
if( pb->id() != particle->id()) continue;
if( !((abs(pd->id()) == abs(children. first->id()) &&
abs(pc->id()) == abs(children.second->id())) ||
(abs(pc->id()) == abs(children. first->id()) &&
abs(pd->id()) == abs(children.second->id()))))
continue;
return vertex;
}
}
}
}
return VertexBasePtr();
}
diff --git a/Models/General/TwoBodyDecayConstructor.h b/Models/General/TwoBodyDecayConstructor.h
--- a/Models/General/TwoBodyDecayConstructor.h
+++ b/Models/General/TwoBodyDecayConstructor.h
@@ -1,177 +1,182 @@
// -*- C++ -*-
//
// TwoBodyDecayConstructor.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_TwoBodyDecayConstructor_H
#define HERWIG_TwoBodyDecayConstructor_H
//
// This is the declaration of the TwoBodyDecayConstructor class.
//
#include "NBodyDecayConstructorBase.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "Herwig/Decay/General/GeneralTwoBodyDecayer.fh"
#include "Herwig/Shower/ShowerAlpha.h"
#include "TwoBodyDecay.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
using Helicity::tVertexBasePtr;
/**
* The TwoBodyDecayConstructor class inherits from the dummy base class
* NBodyDecayConstructorBase and implements the necessary functions in
* order to create the 2 body decay modes for a given set of vertices
* stored in a Model class.
*
* @see \ref TwoBodyDecayConstructorInterfaces "The interfaces"
* defined for TwoBodyDecayConstructor.
* @see NBodyDecayConstructor
**/
class TwoBodyDecayConstructor: public NBodyDecayConstructorBase {
public:
/**
* The default constructor.
*/
- TwoBodyDecayConstructor() : inter_(ShowerInteraction::Both) {
+ TwoBodyDecayConstructor() : inter_(ShowerInteraction::Both), weakMassCut_(-GeV) {
radiationVertices_[ShowerInteraction::QCD] = map<tPDPtr,VertexBasePtr>();
radiationVertices_[ShowerInteraction::QED] = map<tPDPtr,VertexBasePtr>();
}
/**
* Function used to determine allowed decaymodes
*@param part vector of ParticleData pointers containing particles in model
*/
virtual void DecayList(const set<PDPtr> & part);
/**
* Number of outgoing lines. Required for correct ordering.
*/
virtual unsigned int numBodies() const { return 2; }
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
TwoBodyDecayConstructor & operator=(const TwoBodyDecayConstructor &) = delete;
private:
/** @name Functions to create decayers and decaymodes. */
//@{
/**
* Function to create decays
* @param inpart Incoming particle
* @param vert The vertex to create decays for
* @param ilist Which list to search
* @param iv Row number in _theExistingDecayers member
* @return A vector a decay modes
*/
void createModes(tPDPtr inpart, VertexBasePtr vert,
unsigned int ilist,
multiset<TwoBodyDecay> & modes);
/**
* Function to create decayer for specific vertex
* @param decay decay mode for this decay
* member variable
*/
GeneralTwoBodyDecayerPtr createDecayer(TwoBodyDecay decay,
vector<VertexBasePtr> );
/**
* Create decay mode(s) from given part and decay modes
* @param decays The vector of decay modes
* @param decayer The decayer responsible for this decay
*/
void createDecayMode(multiset<TwoBodyDecay> & decays);
//@}
/**
* Get the vertex for QED/QCD radiation
*/
VertexBasePtr radiationVertex(tPDPtr particle,ShowerInteraction inter,
tPDPair children = tPDPair ());
private:
/**
* Map of particles and the vertices which generate their QCD
* radiation
*/
map<ShowerInteraction,map<tPDPtr,VertexBasePtr> > radiationVertices_;
/**
* Default choice for the strong coupling object for hard QCD radiation
*/
ShowerAlphaPtr alphaQCD_;
/**
* Default choice for the strong coupling object for hard QED radiation
*/
ShowerAlphaPtr alphaQED_;
/**
* Which type of corrections to the decays to include
*/
ShowerInteraction inter_;
+
+ /**
+ * Cut off or decays via the weak current
+ */
+ Energy weakMassCut_;
};
}
#endif /* HERWIG_TwoBodyDecayConstructor_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 12:40 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4017539
Default Alt Text
(22 KB)

Event Timeline