diff --git a/EvtGenExternal/EvtHmePythia.hh b/EvtGenExternal/EvtHmePythia.hh new file mode 100644 index 0000000..7460fed --- /dev/null +++ b/EvtGenExternal/EvtHmePythia.hh @@ -0,0 +1,28 @@ +#ifndef EVT_HMEPYTHIA_HH +#define EVT_HMEPYTHIA_HH + +#include "EvtGenBase/EvtDecayAmp.hh" + +class EvtPythiaEngine; +class EvtParticle; + +class EvtHmePythia: public EvtDecayAmp { + +public: + + EvtHmePythia(); + std::string getName() override; + EvtDecayBase* clone() override; + + void decay(EvtParticle *p) override; + void init() override; + void initProbMax() override; + +private: + + EvtPythiaEngine* _pythiaEngine; + int _hmeMode; + +}; + +#endif diff --git a/EvtGenExternal/EvtPythiaEngine.hh b/EvtGenExternal/EvtPythiaEngine.hh index 512872a..d9cd38e 100644 --- a/EvtGenExternal/EvtPythiaEngine.hh +++ b/EvtGenExternal/EvtPythiaEngine.hh @@ -1,98 +1,98 @@ #ifdef EVTGEN_PYTHIA //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 2011 University of Warwick, UK // // Module: EvtPythiaEngine // // Description: Interface to the Pytha 8 external generator // // Modification history: // // John Back April 2011 Module created // //------------------------------------------------------------------------ #ifndef EVTPYTHIAENGINE_HH #define EVTPYTHIAENGINE_HH #include "EvtGenModels/EvtAbsExternalGen.hh" #include "EvtGenExternal/EvtPythiaRandom.hh" #include "EvtGenBase/EvtId.hh" #include "EvtGenBase/EvtDecayBase.hh" #include "EvtGenBase/EvtParticle.hh" #include "EvtGenBase/EvtVector4R.hh" #include "Pythia8/Pythia.h" #include "Pythia8/ParticleData.h" #include #include #include #include class EvtPythiaEngine : public EvtAbsExternalGen { public: EvtPythiaEngine(std::string xmlDir = "./xmldoc", bool convertPhysCodes = false, bool useEvtGenRandom = true); virtual ~EvtPythiaEngine(); bool doDecay(EvtParticle* theMother) override; - void doTauDecay(EvtParticle* theParticle, int &meMode, - EvtComplex &1, EvtComplex &2); + double hmeAmplitudes(int hmeMode, EvtParticle* mom = nullptr, + EvtDecayAmp* amp = nullptr); void initialise() override; protected: private: void updateParticleLists(); void updatePhysicsParameters(); void createPythiaParticle(EvtId& particleId, int PDGCode); bool validPDGCode(int PDGCode); void updatePythiaDecayTable(EvtId& particleId, int aliasInt, int PDGCode); void storeDaughterInfo(EvtParticle* theParticle, int startInt); void clearDaughterVectors(); void clearPythiaModeMap(); void createDaughterEvtParticles(EvtParticle* theParent); int getModeInt(EvtDecayBase* decayModel); std::unique_ptr _genericPythiaGen; std::unique_ptr _aliasPythiaGen; Pythia8::Pythia* _thePythiaGenerator; std::vector _daugPDGVector; std::vector _daugP4Vector; typedef std::map > PythiaModeMap; PythiaModeMap _pythiaModeMap; bool _convertPhysCodes, _initialised, _useEvtGenRandom; std::unique_ptr _evtgenRandom; std::map _addedPDGCodes; - std::map > _tauMEs; + std::map > _hmes; }; #endif #endif diff --git a/EvtGenExternal/EvtTauPythia.hh b/EvtGenExternal/EvtTauPythia.hh deleted file mode 100644 index e421113..0000000 --- a/EvtGenExternal/EvtTauPythia.hh +++ /dev/null @@ -1,30 +0,0 @@ -#ifndef EVT_TAUPYTHIA_HH -#define EVT_TAUPYTHIA_HH - -#include "EvtGenBase/EvtDecayIncoherent.hh" -#include - -class EvtPythiaEngine; -class EvtParticle; - -class EvtTauPythia: public EvtDecayIncoherent { - -public: - - EvtTauPythia(); - std::string getName() override; - EvtDecayBase* clone() override; - - void decay(EvtParticle *p) override; - void init() override; - void initProbMax() override; - -protected: - - EvtPythiaEngine* _pythiaEngine; - -private: - -}; - -#endif diff --git a/doc/doxygen.cfg b/doc/doxygen.cfg new file mode 100644 index 0000000..50eba61 --- /dev/null +++ b/doc/doxygen.cfg @@ -0,0 +1,209 @@ +# Copyright (C) 2019 Torbjorn Sjostrand. +# PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. +# Please respect the MCnet Guidelines, see GUIDELINES for details. +# Author: Philip Ilten, January 2015. + +# Configuration file to create Doxygen documentation for Pythia. Based +# on the default Doxyfile 1.6.1. + +DOXYFILE_ENCODING = UTF-8 +PROJECT_NAME = EvtGen +PROJECT_NUMBER = tauDecay +OUTPUT_DIRECTORY = +CREATE_SUBDIRS = NO +OUTPUT_LANGUAGE = English +BRIEF_MEMBER_DESC = YES +REPEAT_BRIEF = YES +ABBREVIATE_BRIEF = +ALWAYS_DETAILED_SEC = NO +INLINE_INHERITED_MEMB = NO +FULL_PATH_NAMES = YES +STRIP_FROM_PATH = doxygen/ +STRIP_FROM_INC_PATH = +SHORT_NAMES = NO +JAVADOC_AUTOBRIEF = NO +QT_AUTOBRIEF = NO +MULTILINE_CPP_IS_BRIEF = NO +INHERIT_DOCS = YES +SEPARATE_MEMBER_PAGES = NO +TAB_SIZE = 8 +ALIASES = +OPTIMIZE_OUTPUT_FOR_C = NO +OPTIMIZE_OUTPUT_JAVA = NO +OPTIMIZE_FOR_FORTRAN = NO +OPTIMIZE_OUTPUT_VHDL = NO +EXTENSION_MAPPING = cc=C++ h=C++ +BUILTIN_STL_SUPPORT = NO +CPP_CLI_SUPPORT = NO +SIP_SUPPORT = NO +IDL_PROPERTY_SUPPORT = YES +DISTRIBUTE_GROUP_DOC = NO +SUBGROUPING = YES +TYPEDEF_HIDES_STRUCT = NO +EXTRACT_ALL = NO +EXTRACT_PRIVATE = NO +EXTRACT_STATIC = NO +EXTRACT_LOCAL_CLASSES = YES +EXTRACT_LOCAL_METHODS = NO +EXTRACT_ANON_NSPACES = NO +HIDE_UNDOC_MEMBERS = NO +HIDE_UNDOC_CLASSES = NO +HIDE_FRIEND_COMPOUNDS = NO +HIDE_IN_BODY_DOCS = NO +INTERNAL_DOCS = NO +CASE_SENSE_NAMES = YES +HIDE_SCOPE_NAMES = YES +SHOW_INCLUDE_FILES = YES +INLINE_INFO = YES +SORT_MEMBER_DOCS = YES +SORT_BRIEF_DOCS = NO +SORT_MEMBERS_CTORS_1ST = NO +SORT_GROUP_NAMES = NO +SORT_BY_SCOPE_NAME = NO +GENERATE_TODOLIST = YES +GENERATE_TESTLIST = YES +GENERATE_BUGLIST = YES +GENERATE_DEPRECATEDLIST= YES +ENABLED_SECTIONS = +MAX_INITIALIZER_LINES = 30 +SHOW_USED_FILES = YES +SHOW_FILES = YES +SHOW_NAMESPACES = NO +FILE_VERSION_FILTER = +LAYOUT_FILE = +QUIET = NO +WARNINGS = NO +WARN_IF_UNDOCUMENTED = NO +WARN_IF_DOC_ERROR = NO +WARN_NO_PARAMDOC = NO +WARN_FORMAT = "$file:$line: $text" +WARN_LOGFILE = +INPUT = doxygen/EvtGen doxygen/EvtGenBase doxygen/EvtGenExternal doxygen/EvtGenModels doxygen/src doxygen/src/EvtGenBase doxygen/src/EvtGenModels doxygen/src/EvtGenExternal +INPUT_ENCODING = UTF-8 +FILE_PATTERNS = *.cpp *.hh +RECURSIVE = NO +EXCLUDE_SYMLINKS = NO +EXCLUDE_PATTERNS = +EXCLUDE_SYMBOLS = +EXAMPLE_PATH = +EXAMPLE_PATTERNS = +EXAMPLE_RECURSIVE = NO +IMAGE_PATH = +INPUT_FILTER = +FILTER_PATTERNS = +FILTER_SOURCE_FILES = NO +SOURCE_BROWSER = NO +INLINE_SOURCES = NO +STRIP_CODE_COMMENTS = NO +REFERENCED_BY_RELATION = NO +REFERENCES_RELATION = NO +REFERENCES_LINK_SOURCE = YES +USE_HTAGS = NO +VERBATIM_HEADERS = YES +ALPHABETICAL_INDEX = NO +COLS_IN_ALPHA_INDEX = 5 +IGNORE_PREFIX = +GENERATE_HTML = YES +HTML_OUTPUT = doxygen +HTML_FILE_EXTENSION = .html +HTML_HEADER = +HTML_FOOTER = +HTML_TIMESTAMP = NO +HTML_STYLESHEET = +HTML_DYNAMIC_SECTIONS = NO +GENERATE_DOCSET = NO +DOCSET_FEEDNAME = "Doxygen generated docs" +DOCSET_BUNDLE_ID = org.doxygen.Project +GENERATE_HTMLHELP = NO +CHM_FILE = +HHC_LOCATION = +GENERATE_CHI = NO +CHM_INDEX_ENCODING = +BINARY_TOC = NO +TOC_EXPAND = NO +GENERATE_QHP = NO +QCH_FILE = +QHP_NAMESPACE = +QHP_VIRTUAL_FOLDER = doc +QHP_CUST_FILTER_NAME = +QHP_CUST_FILTER_ATTRS = +QHP_SECT_FILTER_ATTRS = +QHG_LOCATION = +DISABLE_INDEX = NO +ENUM_VALUES_PER_LINE = 4 +GENERATE_TREEVIEW = NO +TREEVIEW_WIDTH = 250 +FORMULA_FONTSIZE = 10 +SEARCHENGINE = YES +GENERATE_LATEX = NO +LATEX_OUTPUT = latex +LATEX_CMD_NAME = latex +MAKEINDEX_CMD_NAME = makeindex +COMPACT_LATEX = NO +PAPER_TYPE = a4wide +EXTRA_PACKAGES = +LATEX_HEADER = +PDF_HYPERLINKS = YES +USE_PDFLATEX = YES +LATEX_BATCHMODE = NO +LATEX_HIDE_INDICES = NO +LATEX_SOURCE_CODE = NO +GENERATE_RTF = NO +RTF_OUTPUT = rtf +COMPACT_RTF = NO +RTF_HYPERLINKS = NO +RTF_STYLESHEET_FILE = +RTF_EXTENSIONS_FILE = +GENERATE_MAN = NO +MAN_OUTPUT = man +MAN_EXTENSION = .3 +MAN_LINKS = NO +GENERATE_XML = NO +XML_OUTPUT = xml +XML_PROGRAMLISTING = YES +GENERATE_AUTOGEN_DEF = NO +GENERATE_PERLMOD = NO +PERLMOD_LATEX = NO +PERLMOD_PRETTY = YES +PERLMOD_MAKEVAR_PREFIX = +ENABLE_PREPROCESSING = YES +MACRO_EXPANSION = NO +EXPAND_ONLY_PREDEF = NO +SEARCH_INCLUDES = YES +INCLUDE_PATH = +INCLUDE_FILE_PATTERNS = +PREDEFINED = +EXPAND_AS_DEFINED = +SKIP_FUNCTION_MACROS = YES +TAGFILES = +GENERATE_TAGFILE = +ALLEXTERNALS = NO +EXTERNAL_GROUPS = YES +PERL_PATH = /usr/bin/perl +CLASS_DIAGRAMS = YES +MSCGEN_PATH = +HIDE_UNDOC_RELATIONS = YES +HAVE_DOT = NO +DOT_FONTNAME = +DOT_FONTSIZE = 10 +DOT_FONTPATH = +CLASS_GRAPH = YES +COLLABORATION_GRAPH = YES +GROUP_GRAPHS = YES +UML_LOOK = NO +TEMPLATE_RELATIONS = NO +INCLUDE_GRAPH = YES +INCLUDED_BY_GRAPH = YES +CALL_GRAPH = NO +CALLER_GRAPH = NO +GRAPHICAL_HIERARCHY = YES +DIRECTORY_GRAPH = YES +DOT_IMAGE_FORMAT = png +DOT_PATH = +DOTFILE_DIRS = +DOT_GRAPH_MAX_NODES = 50 +MAX_DOT_GRAPH_DEPTH = 0 +DOT_TRANSPARENT = NO +DOT_MULTI_TARGETS = NO +GENERATE_LEGEND = YES +DOT_CLEANUP = YES diff --git a/doc/doxygen.sh b/doc/doxygen.sh new file mode 100755 index 0000000..bbd0d6e --- /dev/null +++ b/doc/doxygen.sh @@ -0,0 +1,59 @@ +#!/usr/bin/env bash +# Copyright (C) 2019 Torbjorn Sjostrand. +# PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. +# Please respect the MCnet Guidelines, see GUIDELINES for details. +# Author: Philip Ilten, January 2015. + +# This is a simple script to create Doxygen documentation for +# PYTHIA. It requires a single argument which provides the source +# directory. Modified to generate EvtGen Doxygen. + +# Grab the arguments. +if [ $# -ne 2 ]; then + echo "usage: doxygen.sh " + exit +fi +VER=$1 +SRC=$2 + +# Check Doxygen and SED exist. +if ! type "sed" &> /dev/null; then + echo "Error: SED not found."; exit; fi +if ! type "doxygen" &> /dev/null; then + echo "Error: DOXYGEN not found."; exit; fi + +# Handle in-place SED for OSX. +ARCH=$(uname | grep -i -o -e Linux -e Darwin) +ARCH=$(echo $ARCH | awk '{print toupper($0)}') +if [ "$ARCH" == "DARWIN" ]; then SED="sed -i .sed"; else SED="sed -i"; fi + +# Determine the directory structure. +DIRS="EvtGen EvtGenBase EvtGenExternal EvtGenModels" +DIRS+=" src/EvtGenBase src/EvtGenExternal src/EvtGenModels" + +# Remove and create new Doxygen structure. +rm -rf doxygen +FILES="" +for DIR in $DIRS; do + INPUT+=" doxygen/$DIR" + mkdir -p doxygen/$DIR + cp $SRC/$DIR/* doxygen/$DIR/ + FILES+=" "`find doxygen/$DIR/*` +done +INPUT=`find doxygen/* -type d | tr '\r\n' ' '` + +# Convert comments to Doxygen style comments. +$SED "s|//\+|//|g" $FILES # Ensure no more than // (i.e. ///). +$SED "s|//|///|g" $FILES # Change all // to ///. +$SED "s|\( *\)///--|\1//--|g" $FILES # Change ///- lines to //-. +$SED "s|\( *\)///==|\1//==|g" $FILES # Change ///= lines to //=. + +# Modify doxygen.cfg. +$SED "s|\(PROJECT_NUMBER *= *\)[^ ]*|\1$VER|g" doxygen.cfg # Version. +$SED "s|\(INPUT *= *\).*|\1$INPUT|g" doxygen.cfg # Inputs. + +# Create the documentation and clean up. +doxygen doxygen.cfg +$SED "s|///|//|g" doxygen/*source.html # Convert source comments back to //. +for DIR in $DIRS; do rm -rf doxygen/$DIR; done +rm -rf doxygen.cfg.sed doxygen/include diff --git a/src/EvtGenExternal/EvtHmePythia.cpp b/src/EvtGenExternal/EvtHmePythia.cpp new file mode 100644 index 0000000..621610c --- /dev/null +++ b/src/EvtGenExternal/EvtHmePythia.cpp @@ -0,0 +1,78 @@ +//-------------------------------------------------------------------------- +// +// Environment: +// This software is part of the EvtGen package. If you use all or part +// of it, please give an appropriate acknowledgement. +// +// Copyright Information: See EvtGen/COPYRIGHT +// +// Module: EvtHmePythia.cpp +// +// Description: Helicity matrix element decays using Pythia8. +// +// Modification history: +// +// PI May 30, 2019 Module created +// +//------------------------------------------------------------------------ +// +#include "EvtGenBase/EvtParticle.hh" +#include "EvtGenExternal/EvtExternalGenFactory.hh" +#include "EvtGenExternal/EvtPythiaEngine.hh" +#include "EvtGenExternal/EvtHmePythia.hh" + +// Constructor. +EvtHmePythia::EvtHmePythia() : _pythiaEngine(0) {} + +// Name. +std::string EvtHmePythia::getName() { + + return "TAUPYTHIA"; + +} + +// Clone the model. +EvtDecayBase* EvtHmePythia::clone() { + + return new EvtHmePythia(); + +} + +// Initialize the model. +void EvtHmePythia::init() { + + // Retrieve the Pythia engine. + if (!_pythiaEngine) { + _pythiaEngine = dynamic_cast + (EvtExternalGenFactory::getInstance()->getGenerator + (EvtExternalGenFactory::PythiaGenId)); + } + + // Check the matrix element is specified and is available. + checkNArg(1); + _hmeMode = getArg(0); + if (_pythiaEngine->hmeAmplitudes(_hmeMode) == 0) + EvtGenReport(EVTGEN_ERROR, "EvtGen") + << "Helicity matrix element mode " << _hmeMode + << " is not available." << endl; + + // Check the number of daughters matches the matrix element. + checkNDaug((_hmeMode % 100)/10); + +} + +void EvtHmePythia::initProbMax() { + + setProbMax(_pythiaEngine->hmeAmplitudes(_hmeMode)); + +} + +void EvtHmePythia::decay(EvtParticle* p) { + + // Initialize the phase-space. + p->initializePhaseSpace(getNDaug(), getDaugs()); + + // Call the Pythia helicity matrix element. + if (_pythiaEngine) _pythiaEngine->hmeAmplitudes(p, _hmeMode, this); + +} diff --git a/src/EvtGenExternal/EvtPythiaEngine.cpp b/src/EvtGenExternal/EvtPythiaEngine.cpp index fec0eec..580dbd4 100644 --- a/src/EvtGenExternal/EvtPythiaEngine.cpp +++ b/src/EvtGenExternal/EvtPythiaEngine.cpp @@ -1,915 +1,935 @@ #ifdef EVTGEN_PYTHIA //-------------------------------------------------------------------------- // // Environment: // This software is part of the EvtGen package. If you use all or part // of it, please give an appropriate acknowledgement. // // Copyright Information: See EvtGen/COPYRIGHT // Copyright (C) 2011 University of Warwick, UK // // Module: EvtPythiaEngine // // Description: Interface to the Pytha 8 external generator // // Modification history: // // John Back April 2011 Module created -// +// PI May 2019 Added helicity matrix elements //------------------------------------------------------------------------ #include "EvtGenExternal/EvtPythiaEngine.hh" #include "EvtGenBase/EvtPDL.hh" #include "EvtGenBase/EvtDecayTable.hh" #include "EvtGenBase/EvtSpinType.hh" #include "EvtGenBase/EvtParticleFactory.hh" #include "EvtGenBase/EvtReport.hh" #include "EvtGenBase/EvtExtGeneratorCommandsTable.hh" #include "EvtGenExternal/EvtPythia6CommandConverter.hh" #include "Pythia8/Event.h" #include "Pythia8/ParticleData.h" #include #include #include using std::endl; EvtPythiaEngine::EvtPythiaEngine(std::string xmlDir, bool convertPhysCodes, bool useEvtGenRandom) { // Create two Pythia generators. One will be for generic // Pythia decays in the decay.dec file. The other one will be to // only decay aliased particles, which are in general "signal" // decays different from those in the decay.dec file. // Even though it is not possible to have two different particle // versions in one Pythia generator, we can use two generators to // get the required behaviour. EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Creating generic Pythia generator"<(xmlDir); EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Creating alias Pythia generator"<(xmlDir, false); _thePythiaGenerator = 0; _daugPDGVector.clear(); _daugP4Vector.clear(); _convertPhysCodes = convertPhysCodes; // Specify if we are going to use the random number generator (engine) // from EvtGen for Pythia 8. _useEvtGenRandom = useEvtGenRandom; _evtgenRandom = std::make_unique(); _initialised = false; } EvtPythiaEngine::~EvtPythiaEngine() { _thePythiaGenerator = nullptr; this->clearDaughterVectors(); this->clearPythiaModeMap(); } void EvtPythiaEngine::clearDaughterVectors() { _daugPDGVector.clear(); _daugP4Vector.clear(); } void EvtPythiaEngine::clearPythiaModeMap() { PythiaModeMap::iterator iter; for (iter = _pythiaModeMap.begin(); iter != _pythiaModeMap.end(); ++iter) { std::vector modeVector = iter->second; modeVector.clear(); } _pythiaModeMap.clear(); } void EvtPythiaEngine::initialise() { if (_initialised) {return;} this->clearPythiaModeMap(); this->updateParticleLists(); // Hadron-level processes only (hadronized, string fragmentation and secondary decays). // We do not want to generate the full pp or e+e- event structure etc.. _genericPythiaGen->readString("ProcessLevel:all = off"); _aliasPythiaGen->readString("ProcessLevel:all = off"); // Turn off Pythia warnings, e.g. changes to particle properties _genericPythiaGen->readString("Print:quiet = on"); _aliasPythiaGen->readString("Print:quiet = on"); // Apply any other physics (or special particle) requirements/cuts etc.. this->updatePhysicsParameters(); // Set the random number generator if (_useEvtGenRandom == true) { _genericPythiaGen->setRndmEnginePtr(_evtgenRandom.get()); _aliasPythiaGen->setRndmEnginePtr(_evtgenRandom.get()); } _genericPythiaGen->init(); _aliasPythiaGen->init(); // Create the tau decay model map. _tauMEs[0] = std::make_unique(); _tauMEs[1521] = std::make_unique(); _tauMEs[1531] = std::make_unique(); _tauMEs[1532] = std::make_unique(); _tauMEs[1533] = std::make_unique(); _tauMEs[1541] = std::make_unique(); _tauMEs[1542] = std::make_unique(); _tauMEs[1543] = std::make_unique(); _tauMEs[1544] = std::make_unique(); _tauMEs[1551] = std::make_unique(); _tauMEs[1561] = std::make_unique(); // Initialize the tau decay model pointers. for (auto &model: _tauMEs) { model.second->initPointers(&_genericPythiaGen->particleData, &_genericPythiaGen->couplings); } _initialised = true; } bool EvtPythiaEngine::doDecay(EvtParticle* theParticle) { // Store the mother particle within a Pythia8 Event object. // Then do the hadron level decays. // The EvtParticle must be a colour singlet (meson/baryon/lepton), i.e. not a gluon or quark // We delete any daughters the particle may have, since we are asking Pythia // to generate the decay anew. Also note that _any_ Pythia decay allowed for the particle // will be generated and not the specific Pythia decay mode that EvtGen has already // specified. This is necessary since we only want to initialise the Pythia decay table // once; all Pythia branching fractions for a given mother particle are renormalised to sum to 1.0. // In EvtGen decay.dec files, it may be the case that Pythia decays are only used // for some of the particle decays (i.e. Pythia BF sum < 1.0). As we loop over many events, // the total frequency for each Pythia decay mode will normalise correctly to what // we wanted via the specifications made to the decay.dec file, even though event-by-event // the EvtGen decay channel and the Pythia decay channel may be different. if (_initialised == false) {this->initialise();} if (theParticle == 0) { EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Error in EvtPythiaEngine::doDecay. The mother particle is null. Not doing any Pythia decay."<getNDaug() != 0) { bool keepChannel(false); theParticle->deleteDaughters(keepChannel); } EvtId particleId = theParticle->getId(); int isAlias = particleId.isAlias(); // Choose the generator depending if we have an aliased (parent) particle or not _thePythiaGenerator = (isAlias == 1 ? _aliasPythiaGen.get(): _genericPythiaGen.get() ); // Need to use the reference to the Pythia8::Event object, // otherwise it will just return a new empty, default event object. Pythia8::Event& theEvent = _thePythiaGenerator->event; theEvent.reset(); // Initialise the event to be the particle rest frame int PDGCode = EvtPDL::getStdHep(particleId); int status(1); int colour(0), anticolour(0); double px(0.0), py(0.0), pz(0.0); double m0 = theParticle->mass(); double E = m0; theEvent.append(PDGCode, status, colour, anticolour, px, py, pz, E, m0); // Generate the Pythia event int iTrial(0); bool generatedEvent(false); for (iTrial = 0; iTrial < 10; iTrial++) { generatedEvent = _thePythiaGenerator->next(); if (generatedEvent) {break;} } bool success(false); if (generatedEvent) { // Store the daughters for this particle from the Pythia decay tree. // This is a recursive function that will continue looping through // all available daughters until the first set of non-quark and non-gluon // particles are encountered in the Pythia Event structure. // First, clear up the internal vectors storing the daughter // EvtId types and 4-momenta. this->clearDaughterVectors(); // Now store the daughter info. Since this is a recursive function // to loop through the full Pythia decay tree, we do not want to create // EvtParticles here but in the next step. this->storeDaughterInfo(theParticle, 1); // Now create the EvtParticle daughters of the (parent) particle. // We need to use the EvtParticle::makeDaughters function // owing to the way EvtParticle stores parent-daughter information. this->createDaughterEvtParticles(theParticle); //theParticle->printTree(); //theEvent.list(true, true); success = true; } return success; } void EvtPythiaEngine::storeDaughterInfo(EvtParticle* theParticle, int startInt) { Pythia8::Event& theEvent = _thePythiaGenerator->event; std::vector daugList = theEvent.daughterList(startInt); std::vector::iterator daugIter; for (daugIter = daugList.begin(); daugIter != daugList.end(); ++daugIter) { int daugInt = *daugIter; // Ask if the daughter is a quark or gluon. If so, recursively search again. Pythia8::Particle& daugParticle = theEvent[daugInt]; if (daugParticle.isQuark() || daugParticle.isGluon()) { // Recursively search for correct daughter type this->storeDaughterInfo(theParticle, daugInt); } else { // We have a daughter that is not a quark nor gluon particle. // Make sure we are not double counting particles, since several quarks // and gluons make one particle. // Set the status flag for any "new" particle to say that we have stored it. // Use status flag = 1000 (within the user allowed range for Pythia codes). // Check that the status flag for the particle is not equal to 1000 int statusCode = daugParticle.status(); if (statusCode != 1000) { int daugPDGInt = daugParticle.id(); double px = daugParticle.px(); double py = daugParticle.py(); double pz = daugParticle.pz(); double E = daugParticle.e(); EvtVector4R daughterP4(E, px, py, pz); // Now store the EvtId and 4-momentum in the internal vectors _daugPDGVector.push_back(daugPDGInt); _daugP4Vector.push_back(daughterP4); // Set the status flag for the Pythia particle to let us know // that we have already considered it to avoid double counting. daugParticle.status(1000); } // Status code != 1000 } } } void EvtPythiaEngine::createDaughterEvtParticles(EvtParticle* theParent) { if (theParent == 0) { EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Error in EvtPythiaEngine::createDaughterEvtParticles. The parent is null"< daugAliasIdVect(0); EvtId particleId = theParent->getId(); // Check to see if we have an anti-particle. If we do, charge conjugate the particle id to get the // Pythia "alias" we can compare with the defined (particle) Pythia modes. int PDGId = EvtPDL::getStdHep(particleId); int aliasInt = particleId.getAlias(); int pythiaAliasInt(aliasInt); if (PDGId < 0) { // We have an anti-particle. EvtId conjPartId = EvtPDL::chargeConj(particleId); pythiaAliasInt = conjPartId.getAlias(); } std::vector pythiaModes = _pythiaModeMap[pythiaAliasInt]; // Loop over all available Pythia decay modes and find the channel that matches // the daughter ids. Set each daughter id to also use the alias integer. // This will then convert the Pythia generated channel to the EvtGen alias defined one. std::vector::iterator modeIter; bool gotMode(false); for (modeIter = pythiaModes.begin(); modeIter != pythiaModes.end(); ++modeIter) { // Stop the loop if we have the right decay mode channel if (gotMode) {break;} int pythiaModeInt = *modeIter; EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, pythiaModeInt); if (decayModel != 0) { int nModeDaug = decayModel->getNDaug(); // We need to make sure that the number of daughters match if (nDaughters == nModeDaug) { int iModeDaug(0); for (iModeDaug = 0; iModeDaug < nModeDaug; iModeDaug++) { EvtId daugId = decayModel->getDaug(iModeDaug); int daugPDGId = EvtPDL::getStdHep(daugId); // Pythia has used the right PDG codes for this decay mode, even for conjugate modes int pythiaPDGId = _daugPDGVector[iModeDaug]; if (daugPDGId == pythiaPDGId) { daugAliasIdVect.push_back(daugId); } } // Loop over EvtGen mode daughters int daugAliasSize = daugAliasIdVect.size(); if (daugAliasSize == nDaughters) { // All daughter Id codes are accounted for. Set the flag to stop the loop. gotMode = true; } else { // We do not have the correct daughter ordering. Clear the id vector // and try another mode. daugAliasIdVect.clear(); } } // Same number of daughters } // decayModel != 0 } // Loop over available Pythia modes if (gotMode == false) { // We did not find a match for the daughter aliases. Just use the normal PDG codes // from the Pythia decay result int iPyDaug(0); for (iPyDaug = 0; iPyDaug < nDaughters; iPyDaug++) { int daugPDGCode = _daugPDGVector[iPyDaug]; EvtId daugPyId = EvtPDL::evtIdFromStdHep(daugPDGCode); daugAliasIdVect.push_back(daugPyId); } } // Make the EvtParticle daughters (with correct alias id's). Their 4-momenta are uninitialised. theParent->makeDaughters(nDaughters, daugAliasIdVect); // Now set the 4-momenta of the daughters. int iDaug(0); // Can use an iterator here, but we already had to use the vector size... for (iDaug = 0; iDaug < nDaughters; iDaug++) { EvtParticle* theDaughter = theParent->getDaug(iDaug); // Set the correct 4-momentum for each daughter particle. if (theDaughter != 0) { EvtId theDaugId = daugAliasIdVect[iDaug]; const EvtVector4R theDaugP4 = _daugP4Vector[iDaug]; theDaughter->init(theDaugId, theDaugP4); } } } void EvtPythiaEngine::updateParticleLists() { // Use the EvtGen decay table (decay/user.dec) to update Pythia particle // decay modes. Also, make sure the generic and alias Pythia generators // use the same particle data entries as defined by EvtGen (evt.pdl). // Loop over all entries in the EvtPDL particle data table. // Aliases are added at the end with id numbers equal to the // original particle, but with alias integer = lastPDLEntry+1 etc.. int iPDL; int nPDL = EvtPDL::entries(); // Reset the _addedPDGCodes map that keeps track // of any new particles added to the Pythia input data stream _addedPDGCodes.clear(); for (iPDL = 0; iPDL < nPDL; iPDL++) { EvtId particleId = EvtPDL::getEntry(iPDL); int aliasInt = particleId.getAlias(); // Pythia and EvtGen should use the same PDG codes for particles int PDGCode = EvtPDL::getStdHep(particleId); // Update pole mass, width, lifetime and mass range double mass = EvtPDL::getMeanMass(particleId); double width = EvtPDL::getWidth(particleId); double lifetime = EvtPDL::getctau(particleId); double mmin = EvtPDL::getMinMass(particleId); double mmax = EvtPDL::getMaxMass(particleId); // Particle data pointers. The generic and alias Pythia generator pointers have // their own Pythia8::ParticleData particleData public data members which contain // the particle properties table. We can extract and change individual particle // entries using the particleDataEntryPtr() function within ParticleData. // However, we must be careful when accessing the particleData table. We must do // this directly, since assigning it to another Pythia8::ParticleData object via copying // or assignment will give it a different memory address and it will no longer refer to // the original particleData information from the generator pointer. Pythia8::ParticleDataEntry* entry_generic = _genericPythiaGen->particleData.particleDataEntryPtr(PDGCode); Pythia8::ParticleDataEntry* entry_alias = _aliasPythiaGen->particleData.particleDataEntryPtr(PDGCode); // Check that the PDG code is not zero/null and exclude other // special cases, e.g. those reserved for internal generator use if (entry_generic != 0 && this->validPDGCode(PDGCode)) { entry_generic->setM0(mass); entry_generic->setMWidth(width); entry_generic->setTau0(lifetime); if (std::fabs(width) > 0.0) { entry_generic->setMMin(mmin); entry_generic->setMMax(mmax); } } // Check that the PDG code is not zero/null and exclude other // special cases, e.g. those reserved for internal generator use if (entry_alias != 0 && this->validPDGCode(PDGCode)) { entry_alias->setM0(mass); entry_alias->setMWidth(width); entry_alias->setTau0(lifetime); if (std::fabs(width) > 0.0) { entry_alias->setMMin(mmin); entry_alias->setMMax(mmax); } } // Check which particles have a Pythia decay defined. // Get the list of all possible decays for the particle, using the alias integer. // If the particle is not actually an alias, aliasInt = idInt. bool hasPythiaDecays = EvtDecayTable::getInstance()->hasPythia(aliasInt); if (hasPythiaDecays) { int isAlias = particleId.isAlias(); // Decide what generator to use depending on whether we have // an aliased particle or not _thePythiaGenerator = ( isAlias==1 ? _aliasPythiaGen.get() : _genericPythiaGen.get() ); // Find the Pythia particle name given the standard PDG code integer std::string dataName = _thePythiaGenerator->particleData.name(PDGCode); bool alreadyStored = (_addedPDGCodes.find(abs(PDGCode)) != _addedPDGCodes.end()); if (dataName == " " && !alreadyStored ) { // Particle and its antiparticle do not exist in the Pythia database. // Create a new particle, then create the new decay modes. this->createPythiaParticle(particleId, PDGCode); } // For the particle, create the Pythia decay modes. // Update Pythia data tables. this->updatePythiaDecayTable(particleId, aliasInt, PDGCode); } // Loop over Pythia decays } // Loop over EvtPDL entries //EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed generic Pythia decay list"<particleData.listChanged(); //EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed alias Pythia decay list"<particleData.listChanged(); } bool EvtPythiaEngine::validPDGCode(int PDGCode) { // Exclude certain PDG codes: void = 0 and special values = 81 to 100, which are reserved // for internal generator use (pseudoparticles) according to PDG guidelines. Also exclude // nu'_tau (nu_L) = 18, which has different masses: Pythia8 = 400 GeV, EvtGen = 0 GeV. bool isValid(true); int absPDGCode = abs(PDGCode); if (absPDGCode == 0 || absPDGCode == 18) { // Void and nu_L or nu'_tau isValid = false; } else if (absPDGCode >= 81 && absPDGCode <= 100) { // Pseudoparticles isValid = false; } return isValid; } void EvtPythiaEngine::updatePythiaDecayTable(EvtId& particleId, int aliasInt, int PDGCode) { // Update the particle data table in Pythia. // The tables store information about the allowed decay modes // where the PDGId for all particles must be positive; anti-particles are stored // with the corresponding particle entry. // Since we do not want to implement CP violation here, just use the same branching // fractions for particle and anti-particle modes. int nModes = EvtDecayTable::getInstance()->getNModes(aliasInt); int iMode(0); bool firstMode(true); // Only process positive PDG codes. if (PDGCode < 0) {return;} // Keep track of which decay modes are Pythia decays for each aliasInt std::vector pythiaModes(0); // Loop over the decay modes for this particle for (iMode = 0; iMode < nModes; iMode++) { EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, iMode); if (decayModel != 0) { int nDaug = decayModel->getNDaug(); // If the decay mode has no daughters, then that means that there will be // no entries for any submode re-definitions for Pythia. // This sometimes occurs for any mode using non-standard Pythia 6 codes. // Do not refine the decay mode, i.e. accept the Pythia 8 default (if it exists). if (nDaug > 0) { // Check to see if we have a Pythia decay mode std::string modelName = decayModel->getModelName(); if (modelName == "PYTHIA") { // Keep track which decay mode is a Pythia one. We need this in order to // reassign alias Id values for particles generated in the decay. pythiaModes.push_back(iMode); std::ostringstream oss; oss.setf(std::ios::scientific); // Write out the absolute value of the PDG code, since // particles and anti-particles occupy the same part of the Pythia table. oss << PDGCode; if (firstMode) { // Create a new channel oss <<":oneChannel = "; firstMode = false; } else { // Add the channel oss <<":addChannel = "; } // Select all channels (particle and anti-particle). // For CP violation, or different BFs for particle and anti-particle, // use options 2 or 3 (not here). int onMode(1); oss << onMode << " "; double BF = decayModel->getBranchingFraction(); oss << BF << " "; // Need to convert the old Pythia physics mode integers with the new ones // To do this, get the model argument and write a conversion method. int modeInt = this->getModeInt(decayModel); oss << modeInt; int iDaug(0); for (iDaug = 0; iDaug < nDaug; iDaug++) { EvtId daugId = decayModel->getDaug(iDaug); int daugPDG = EvtPDL::getStdHep(daugId); oss << " " << daugPDG; } // Daughter list _thePythiaGenerator->readString(oss.str()); } // is Pythia } else { EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Warning in EvtPythiaEngine. Trying to redefine Pythia table for " <readString(rescaleStr.str()); } int EvtPythiaEngine::getModeInt(EvtDecayBase* decayModel) { int tmpModeInt(0), modeInt(0); if (decayModel != 0) { int nVars = decayModel->getNArg(); // Just read the first integer, which specifies the Pythia decay model. // Ignore any other values. if (nVars > 0) { tmpModeInt = static_cast(decayModel->getArg(0)); } } if (_convertPhysCodes) { // Extra code to convert the old Pythia decay model integer MDME(ICC,2) to the new one. // This should be removed eventually after updating decay.dec files to use // the new convention. if (tmpModeInt == 0) { modeInt = 0; // phase-space } else if (tmpModeInt == 1) { modeInt = 1; // omega or phi -> 3pi } else if (tmpModeInt == 2) { modeInt = 11; // Dalitz decay } else if (tmpModeInt == 3) { modeInt = 2; // V -> PS PS } else if (tmpModeInt == 4) { modeInt = 92; // onium -> ggg or gg gamma } else if (tmpModeInt == 11) { modeInt = 42; // phase-space of hadrons from available quarks } else if (tmpModeInt == 12) { modeInt = 42; // phase-space for onia resonances } else if (tmpModeInt == 13) { modeInt = 43; // phase-space of at least 3 hadrons } else if (tmpModeInt == 14) { modeInt = 44; // phase-space of at least 4 hadrons } else if (tmpModeInt == 15) { modeInt = 45; // phase-space of at least 5 hadrons } else if (tmpModeInt >= 22 && tmpModeInt <= 30) { modeInt = tmpModeInt + 40; // phase space of hadrons with fixed multiplicity (modeInt - 60) } else if (tmpModeInt == 31) { modeInt = 42; // two or more quarks phase-space; one spectactor quark } else if (tmpModeInt == 32) { modeInt = 91; // qqbar or gg pair } else if (tmpModeInt == 33) { modeInt = 0; // triplet q X qbar, where X = gluon or colour singlet (superfluous, since g's are created anyway) } else if (tmpModeInt == 41) { modeInt = 21; // weak decay phase space, weighting nu_tau spectrum } else if (tmpModeInt == 42) { modeInt = 22; // weak decay V-A matrix element } else if (tmpModeInt == 43) { modeInt = 22; // weak decay V-A matrix element, quarks as jets (superfluous) } else if (tmpModeInt == 44) { modeInt = 22; // weak decay V-A matrix element, parton showers (superfluous) } else if (tmpModeInt == 48) { modeInt = 23; // weak decay V-A matrix element, at least 3 decay products } else if (tmpModeInt == 50) { modeInt = 0; // default behaviour } else if (tmpModeInt == 51) { modeInt = 0; // step threshold (channel switched off when mass daughters > mother mass } else if (tmpModeInt == 52 || tmpModeInt == 53) { modeInt = 0; // beta-factor threshold } else if (tmpModeInt == 84) { modeInt = 42; // unknown physics process - just use phase-space } else if (tmpModeInt == 101) { modeInt = 0; // continuation line } else if (tmpModeInt == 102) { modeInt = 0; // off mass shell particles. } else { EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Pythia mode integer "< 0) { colour = 1; // single quarks } double m0 = EvtPDL::getMeanMass(particleId); double mWidth = EvtPDL::getWidth(particleId); double mMin = EvtPDL::getMinMass(particleId); double mMax = EvtPDL::getMaxMass(particleId); double tau0 = EvtPDL::getctau(particleId); std::ostringstream oss; oss.setf(std::ios::scientific); int absPDGCode = abs(PDGCode); oss << absPDGCode << ":new = " << aliasName << " " << antiName << " " << spin << " " << charge << " " << colour << " " << m0 << " " << mWidth << " " << mMin << " " << mMax << " " << tau0; // Pass this information to Pythia _thePythiaGenerator->readString(oss.str()); // Also store the absolute value of the PDG entry // to keep track of which new particles have been added, // which also automatically includes the anti-particle. // We need to avoid creating new anti-particles when // they already exist when the particle was added. _addedPDGCodes[absPDGCode] = 1; } void EvtPythiaEngine::updatePhysicsParameters() { // Update any more Pythia physics (or special particle) requirements/cuts etc.. // This should be used if any of the Pythia 6 parameters like JetSetPar MSTJ(i) = x // are needed. Such commands will need to be implemented using the new interface // pythiaGenerator->readString(cmd); Here cmd is a string telling Pythia 8 // what physics parameters to change. This will need to be done for the generic and // alias generator pointers, as appropriate. // Set the multiplicity level for hadronic weak decays std::string multiWeakCut("ParticleDecays:multIncreaseWeak = 2.0"); _genericPythiaGen->readString(multiWeakCut); _aliasPythiaGen->readString(multiWeakCut); // Set the multiplicity level for all other decays std::string multiCut("ParticleDecays:multIncrease = 4.5"); _genericPythiaGen->readString(multiCut); _aliasPythiaGen->readString(multiCut); //Now read in any custom configuration entered in the XML GeneratorCommands commands = EvtExtGeneratorCommandsTable::getInstance()->getCommands("PYTHIA"); GeneratorCommands::iterator it = commands.begin(); for( ; it!=commands.end(); it++) { Command command = *it; std::vector commandStrings; if(command["VERSION"] == "PYTHIA6") { EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Converting Pythia 6 command: "<::iterator it2 = commandStrings.begin(); for( ; it2!=commandStrings.end(); it2++) { EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Configuring generic Pythia generator: " << (*it2) << endl; _genericPythiaGen->readString(*it2); } } if(generator == "ALIAS" || generator == "Alias" || generator == "alias" || generator == "BOTH" || generator == "Both" || generator == "both") { std::vector::iterator it2 = commandStrings.begin(); for( ; it2!=commandStrings.end(); it2++) { EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Configuring alias Pythia generator: " << (*it2) << endl; _aliasPythiaGen->readString(*it2); } } } } -void EvtPythiaEngine::doTauDecay(EvtParticle* theParticle, int &meMode, - EvtComplex &1, EvtComplex &2) { +double EvtPythiaEngine::hmeAmplitudes(int hmeMode, EvtParticle* mom, + EvtDecayAmp *amp) { - // Get the helicity matrix element. - std::unique_ptr &me = _tauMEs[meMode]; - if (me == nullptr) return; + // Get the helicity matrix element (return 0 if no matrix element). + std::unique_ptr &hme = _hmes[hmeMode]; + if (hme == nullptr) return 0; + + // Return the maximum amplitude if no particle is given. + vector prts; + if (mom == nullptr) { + prts.push_back(HelicityParticle(mom->PDGId())); + prts[0].rho[0][0] = 1; + return hme->decayWeightMax(); + } - /* // Create the children. - vector children; - - // Decay the tau. - bool accepted = false; - int tries = 0; - while (!accepted) { - //isotropicDecay(children); - double decayWeight = me->decayWeight(children); - double decayWeightMax = me->decayWeightMax(children); - accepted = (rndmPtr->flat() < decayWeight / decayWeightMax); - if (decayWeight > decayWeightMax) - infoPtr->errorMsg("Warning in TauDecays::decay: maximum " - "decay weight exceeded in tau decay"); - if (tries > NTRYDECAY) { - infoPtr->errorMsg("Warning in TauDecays::decay: maximum " - "number of decay attempts exceeded"); - break; + int mult(mom->getNDaug() + 1), ih(0); + vector ihs(3, -1), nhs(3, 0), hs(mult, 0); + EvtVector4R p = mom->getP4Lab(); + prts.push_back(HelicityParticle(mom->PDGId(), 0, 0, 0, 0, 0, 0, 0, + p.get(1), p.get(2), p.get(3), p.get(0), + p.mass(), p.mass(), + &thePythiaGenerator->particleData)); + prts.back().direction = -1; + if (mom->GetSpinStates > 1) { + ihs[ih] = 0; nhs[ih] = mom->GetSpinStates; ++ih;} + for (int iDtr = 0; iDtr < (int)mom->getNDaug(); ++iDtr) { + EvetParticle *dtr = mom->getDaug(iDtr); + p = dtr->getP4Lab(); + prts.push_back(HelicityParticle(dtr->PDGId(), 0, 0, 0, 0, 0, 0, 0, + p.get(1), p.get(2), p.get(3), p.get(0), + p.mass(), p.mass(), + &thePythiaGenerator->particleData)); + prts.back().direction = 1; + states[iDtr + 1] = dtr->getSpinStates(); + if (dtr->GetSpinStates > 1) { + if (ih < 3) {ihs[ih] = 0; nhs[ih] = dtr->GetSpinStates; ++ih;} + else {EvtGenReport(EVTGEN_ERROR, "EvtGen") + << "More than three particles with multiple spins." << endl;} } - ++tries; + if (dtr->spinType() == EvtSpinType::NEUTRINO && dtr->PDGId() < 0) + hs[iDtr + 1] = 1; } - */ - - EvtVector4R p4 = theParticle->getP4(); - - // EvtParticle* parent = theParticle->getParent(); - int PDGId = theParticle->getPDGId(); - - // Get the (forward) spin density matrix - // EvtSpinDensity sd = theParticle->getSpinDensityForward(); - - EvtGenReport(EVTGEN_INFO,"EvtGen")<<"doTauDecay "<< endl; + // Calculate and set the amplitudes (maximum three helicities). + if (ihs[0] < 0) {vertex(hme->calculateME(hs)); + } else { + for (hs[ihs[0]] = 0; hs[ihs[0]] < nhs[0]; ++hs[ihs[0]]) { + if (ihs[1] < 0) {vertex(hs[ihs[0]], hme->calculateME(hs)); + } else { + for (hs[ihs[1]] = 0; hs[ihs[1]] < nhs[1]; ++hs[ihs[1]]) { + if (ihs[2] < 0) {vertex(hs[ihs[0]], hs[ihs[1]], hme->calculateME(hs)); + } else { + for (hs[ihs[2]] = 0; hs[ihs[2]] < nhs[2]; ++hs[ihs[2]]) { + vertex(hs[ihs[0]], hs[ihs[1]], hs[ihs[2]], hme->calculateME(hs)); + } + } // End third helicity. + } + } // End second helicity. + } + } // End first helicity. } #endif diff --git a/src/EvtGenExternal/EvtTauPythia.cpp b/src/EvtGenExternal/EvtTauPythia.cpp deleted file mode 100644 index 2f2cfee..0000000 --- a/src/EvtGenExternal/EvtTauPythia.cpp +++ /dev/null @@ -1,88 +0,0 @@ -//-------------------------------------------------------------------------- -// -// Environment: -// This software is part of the EvtGen package. If you use all or part -// of it, please give an appropriate acknowledgement. -// -// Copyright Information: See EvtGen/COPYRIGHT -// -// Module: EvtTauPythia.cpp -// -// Description: Tau decays using Pythia8 -// -// Modification history: -// -// -//------------------------------------------------------------------------ -// - -#include "EvtGenExternal/EvtTauPythia.hh" - -#include "EvtGenBase/EvtDecayBase.hh" -#include "EvtGenBase/EvtParticle.hh" - -#include "EvtGenExternal/EvtExternalGenFactory.hh" -#include "EvtGenExternal/EvtPythiaEngine.hh" - -EvtTauPythia::EvtTauPythia() : - _pythiaEngine(0) -{} - -std::string EvtTauPythia::getName() { - - return "TAUPYTHIA"; - -} - -EvtDecayBase* EvtTauPythia::clone() { - - return new EvtTauPythia(); - -} - -void EvtTauPythia::init() { - - // Check that there are 0 decay file parameter arguments - checkNArg(0); - - // Check number of daughters: 2 or 3? - // checkNDaug(2, 3); - - // Check that the parent is a tau - checkSpinParent(EvtSpinType::DIRAC); - // Daughter checks? - // checkSpinDaughter(0, EvtSpinType::DIRAC); - // checkSpinDaughter(2,EvtSpinType::NEUTRINO); - -} - -void EvtTauPythia::initProbMax() { - - setProbMax(1.0); - -} - -void EvtTauPythia::decay(EvtParticle* p) { - - // Retrieve the Pythia engine. - if (!_pythiaEngine) { - _pythiaEngine = dynamic_cast - (EvtExternalGenFactory::getInstance()->getGenerator - (EvtExternalGenFactory::PythiaGenId)); - } - - // Initialize the daughters. - p->initializePhaseSpace(getNDaug(), getDaugs()); - - // Call the Pythia tau decay model. - int meMode(0); - EvtComplex amp1(1.0, 0.0), amp2(0.0, 0.0); - if (_pythiaEngine) { - _pythiaEngine->doTauDecay(p, meMode, amp1, amp2); - } - - // Set the helicity amplitude. - //vertex(0, amp1); - //vertex(1, amp2); - -}