No OneTemporary
Actions

Size

39 KB

Subscribers

None

View Options

	diff --git a/EvtGenExternal/EvtPythiaEngine.hh b/EvtGenExternal/EvtPythiaEngine.hh
	index 04b76a8..cc739e8 100644
	--- a/EvtGenExternal/EvtPythiaEngine.hh
	+++ b/EvtGenExternal/EvtPythiaEngine.hh
	@@ -1,92 +1,92 @@

	/***********************************************************************
	* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
	* *
	* This file is part of EvtGen. *
	* *
	* EvtGen is free software: you can redistribute it and/or modify *
	* it under the terms of the GNU General Public License as published by *
	* the Free Software Foundation, either version 3 of the License, or *
	* (at your option) any later version. *
	* *
	* EvtGen is distributed in the hope that it will be useful, *
	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	* GNU General Public License for more details. *
	* *
	* You should have received a copy of the GNU General Public License *
	* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
	***********************************************************************/

	#ifdef EVTGEN_PYTHIA

	#ifndef EVTPYTHIAENGINE_HH
	#define EVTPYTHIAENGINE_HH

	#include "EvtGenBase/EvtDecayBase.hh"
	#include "EvtGenBase/EvtId.hh"
	#include "EvtGenBase/EvtParticle.hh"
	#include "EvtGenBase/EvtVector4R.hh"

	#include "EvtGenModels/EvtAbsExternalGen.hh"

	#include "EvtGenExternal/EvtPythiaRandom.hh"

	#include "Pythia8/ParticleData.h"
	#include "Pythia8/Pythia.h"

	#include <map>
	#include <memory>
	#include <string>
	#include <vector>

	// Description: Interface to the Pytha 8 external generator

	class EvtPythiaEngine : public EvtAbsExternalGen {
	public:
	EvtPythiaEngine( std::string xmlDir = "./xmldoc",
	bool convertPhysCodes = false, bool useEvtGenRandom = true );

	virtual ~EvtPythiaEngine();

	bool doDecay( EvtParticle* theMother ) override;

	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<Pythia8::Pythia> _genericPythiaGen;
	std::unique_ptr<Pythia8::Pythia> _aliasPythiaGen;
	Pythia8::Pythia* _thePythiaGenerator;

	std::vector<int> _daugPDGVector;
	std::vector<EvtVector4R> _daugP4Vector;

	typedef std::map<int, std::vector<int>> PythiaModeMap;
	PythiaModeMap _pythiaModeMap;

	bool _convertPhysCodes, _initialised, _useEvtGenRandom;

	- std::unique_ptr<EvtPythiaRandom> _evtgenRandom;
	+ std::shared_ptr<EvtPythiaRandom> _evtgenRandom;

	std::map<int, int> _addedPDGCodes;
	};

	#endif

	#endif
	diff --git a/src/EvtGenExternal/EvtPythiaEngine.cpp b/src/EvtGenExternal/EvtPythiaEngine.cpp
	index 72aa46d..55c1ce3 100644
	--- a/src/EvtGenExternal/EvtPythiaEngine.cpp
	+++ b/src/EvtGenExternal/EvtPythiaEngine.cpp
	@@ -1,853 +1,858 @@

	/***********************************************************************
	* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
	* *
	* This file is part of EvtGen. *
	* *
	* EvtGen is free software: you can redistribute it and/or modify *
	* it under the terms of the GNU General Public License as published by *
	* the Free Software Foundation, either version 3 of the License, or *
	* (at your option) any later version. *
	* *
	* EvtGen is distributed in the hope that it will be useful, *
	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	* GNU General Public License for more details. *
	* *
	* You should have received a copy of the GNU General Public License *
	* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
	***********************************************************************/

	#ifdef EVTGEN_PYTHIA

	#include "EvtGenExternal/EvtPythiaEngine.hh"

	#include "EvtGenBase/EvtDecayTable.hh"
	#include "EvtGenBase/EvtExtGeneratorCommandsTable.hh"
	#include "EvtGenBase/EvtPDL.hh"
	#include "EvtGenBase/EvtParticleFactory.hh"
	#include "EvtGenBase/EvtReport.hh"
	#include "EvtGenBase/EvtSpinType.hh"

	#include "EvtGenExternal/EvtPythia6CommandConverter.hh"

	#include "Pythia8/Event.h"
	#include "Pythia8/ParticleData.h"

	#include <cmath>
	#include <iostream>
	#include <sstream>

	#if PYTHIA_VERSION_INTEGER < 8304
	typedef Pythia8::ParticleDataEntry* ParticleDataEntryPtr;
	#else
	typedef Pythia8::ParticleDataEntryPtr ParticleDataEntryPtr;
	#endif

	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" << endl;
	_genericPythiaGen = std::make_unique<Pythia8::Pythia>( xmlDir );

	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Creating alias Pythia generator" << endl;
	_aliasPythiaGen = std::make_unique<Pythia8::Pythia>( 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<EvtPythiaRandom>();
	+ _evtgenRandom = std::make_shared<EvtPythiaRandom>();

	_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<int> 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 ) {
	+#if PYTHIA_VERSION_INTEGER < 8310
	_genericPythiaGen->setRndmEnginePtr( _evtgenRandom.get() );
	_aliasPythiaGen->setRndmEnginePtr( _evtgenRandom.get() );
	+#else
	+ _genericPythiaGen->setRndmEnginePtr( _evtgenRandom );
	+ _aliasPythiaGen->setRndmEnginePtr( _evtgenRandom );
	+#endif
	}

	_genericPythiaGen->init();
	_aliasPythiaGen->init();

	_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."
	<< endl;
	return false;
	}

	// Delete EvtParticle daughters if they already exist
	if ( theParticle->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<int> daugList = theEvent.daughterList( startInt );

	std::vector<int>::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"
	<< endl;
	return;
	}

	// Get the list of Pythia decay modes defined for this particle id alias.
	// It would be easier to just use the decay channel number that Pythia chose to use
	// for the particle decay, but this is not accessible from the Pythia interface at present.

	int nDaughters = _daugPDGVector.size();
	std::vector<EvtId> 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<int> 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<int>::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.

	ParticleDataEntryPtr entry_generic =
	_genericPythiaGen->particleData.particleDataEntryPtr( PDGCode );
	ParticleDataEntryPtr 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"<<endl;
	//_genericPythiaGen->particleData.listChanged();

	//EvtGenReport(EVTGEN_INFO,"EvtGen")<<"Writing out changed alias Pythia decay list"<<endl;
	//_aliasPythiaGen->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<int> 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 "
	<< EvtPDL::name( particleId )
	<< " for a decay channel that has no daughters."
	<< " Keeping Pythia default (if available)." << endl;
	}

	} else {
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Error in EvtPythiaEngine. decayModel is null for particle "
	<< EvtPDL::name( particleId ) << " mode number " << iMode
	<< endl;
	}

	} // Loop over modes

	_pythiaModeMap[aliasInt] = pythiaModes;

	// Now, renormalise the decay branching fractions to sum to 1.0
	std::ostringstream rescaleStr;
	rescaleStr.setf( std::ios::scientific );
	rescaleStr << PDGCode << ":rescaleBR = 1.0";

	_thePythiaGenerator->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<int>( 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 " << tmpModeInt
	<< " is not recognised. Using phase-space model" << endl;
	modeInt = 0; // Use phase-space for anything else
	}

	} else {
	// No need to convert the physics mode integer code
	modeInt = tmpModeInt;
	}

	return modeInt;
	}

	void EvtPythiaEngine::createPythiaParticle( EvtId& particleId, int PDGCode )
	{
	// Use the EvtGen name, PDGId and other variables to define the new Pythia particle.
	EvtId antiPartId = EvtPDL::chargeConj( particleId );

	std::string aliasName = EvtPDL::name(
	particleId ); // If not an alias, aliasName = normal name
	std::string antiName = EvtPDL::name( antiPartId );

	EvtSpinType::spintype spinType = EvtPDL::getSpinType( particleId );
	int spin = EvtSpinType::getSpin2( spinType );

	int charge = EvtPDL::chg3( particleId );

	// Must set the correct colour type manually here, since the evt.pdl file
	// does not store this information. This is required for quarks otherwise
	// Pythia cannot generate the decay properly.
	int PDGId = EvtPDL::getStdHep( particleId );
	int colour( 0 );
	if ( PDGId == 21 ) {
	colour = 2; // gluons
	} else if ( PDGId <= 8 && PDGId > 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<std::string> commandStrings;

	if ( command["VERSION"] == "PYTHIA6" ) {
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Converting Pythia 6 command: " << command["MODULE"] << "("
	<< command["PARAM"] << ")=" << command["VALUE"] << "..." << endl;
	commandStrings = convertPythia6Command( command );
	} else if ( command["VERSION"] == "PYTHIA8" ) {
	commandStrings.push_back( command["MODULE"] + ":" + command["PARAM"] +
	" = " + command["VALUE"] );
	} else {
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Pythia command received by EvtPythiaEngine has bad version:"
	<< endl;
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "Received " << command["VERSION"]
	<< " but expected PYTHIA6 or PYTHIA8." << endl;
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "The error is likely to be in EvtDecayTable.cpp" << endl;
	EvtGenReport( EVTGEN_ERROR, "EvtGen" )
	<< "EvtGen will now abort." << endl;
	::abort();
	}
	std::string generator = command["GENERATOR"];
	if ( generator == "GENERIC" \|\| generator == "Generic" \|\|
	generator == "generic" \|\| generator == "BOTH" \|\|
	generator == "Both" \|\| generator == "both" ) {
	std::vector<std::string>::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<std::string>::iterator it2 = commandStrings.begin();
	for ( ; it2 != commandStrings.end(); it2++ ) {
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Configuring alias Pythia generator: " << ( *it2 )
	<< endl;
	_aliasPythiaGen->readString( *it2 );
	}
	}
	}
	}

	#endif

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions