diff --git a/src/EvtGenExternal/EvtHmePythia.cpp b/src/EvtGenExternal/EvtHmePythia.cpp
index 55a9013..002d53c 100644
--- a/src/EvtGenExternal/EvtHmePythia.cpp
+++ b/src/EvtGenExternal/EvtHmePythia.cpp
@@ -1,90 +1,91 @@
/***********************************************************************
* 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 . *
***********************************************************************/
#include "EvtGenExternal/EvtHmePythia.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenExternal/EvtExternalGenFactory.hh"
#include "EvtGenExternal/EvtPythiaEngine.hh"
// Name
std::string EvtHmePythia::getName()
{
return "HMEPYTHIA";
}
// 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.
// Allowed values are (all implicitly include the tau neutrino):
// 1521 = HMETau2Meson
// 1531 = HMETau2TwoLeptons
// 1532 = HMETau2TwoMesonsViaVector
// 1533 = HMETau2TwoMesonsViaVectorScalar
// 1541 = HMETau2ThreePions
// 1542 = HMETau2ThreeMesonsWithKaons
// 1543 = HMETau2ThreeMesonsGeneric
// 1544 = HMETau2TwoPionsGamma
// 1551 = HMETau2FourPions
// 1561 = HMETau2FivePions
checkNArg( 1 );
_hmeMode = getArg( 0 );
if ( _pythiaEngine->hmeAmplitudes( _hmeMode ) == 0 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Helicity matrix element mode " << _hmeMode
<< " is not available." << std::endl;
+ ::abort();
}
// Check the number of daughters matches the matrix element
// e.g. tau- -> e- anti-nu_e nu_tau, HMETau2TwoLeptons (1531), 3 daughters
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( _hmeMode, p, this );
}
}
diff --git a/src/EvtGenExternal/EvtPythiaEngine.cpp b/src/EvtGenExternal/EvtPythiaEngine.cpp
index ebe314d..60261b0 100644
--- a/src/EvtGenExternal/EvtPythiaEngine.cpp
+++ b/src/EvtGenExternal/EvtPythiaEngine.cpp
@@ -1,952 +1,972 @@
/***********************************************************************
* 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 . *
***********************************************************************/
#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
#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" << endl;
_genericPythiaGen = std::make_unique( xmlDir );
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Creating alias Pythia generator" << endl;
_aliasPythiaGen = std::make_unique( 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
_hmes[0] = std::make_unique();
_hmes[1521] = std::make_unique();
_hmes[1531] = std::make_unique();
_hmes[1532] = std::make_unique();
_hmes[1533] = std::make_unique();
_hmes[1541] = std::make_unique();
_hmes[1542] = std::make_unique();
_hmes[1543] = std::make_unique();
_hmes[1544] = std::make_unique();
_hmes[1551] = std::make_unique();
_hmes[1561] = std::make_unique();
// Initialize the tau decay model pointers
for ( auto& model : _hmes ) {
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."
<< 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 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"
<< 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 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 "
<< 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( 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 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::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 );
}
}
}
}
double EvtPythiaEngine::hmeAmplitudes( int hmeMode, EvtParticle* mom,
EvtDecayAmp* amp )
{
+ if ( _initialised == false ) {
+ this->initialise();
+ }
+
// Get the helicity matrix element (return 0 if no matrix element)
std::unique_ptr& hme = _hmes[hmeMode];
if ( hme == nullptr ) {
return 0.0;
}
// Return the maximum amplitude if no particle is given
std::vector prts;
if ( mom == nullptr ) {
- prts.push_back( Pythia8::HelicityParticle( mom->getPDGId() ) );
- prts[0].rho[0][0] = 1;
+ // Use tau particle
+ int PDGCode = EvtPDL::getStdHep( EvtPDL::getId( "tau-" ) );
+ const Pythia8::Particle particle( PDGCode );
+ prts.push_back( Pythia8::HelicityParticle(
+ particle, &_genericPythiaGen->particleData ) );
+
+ // Set first phase space rho element to unity
+ prts[0].rho[0][0] = 1.0;
+
+ // Initialise channel constants, including the maximum decay weight hard limit.
+ // This usually requires the tau mass via Particle::m(), which is initialised
+ // in the Particle constructor. However this mass is set to zero since we
+ // only pass the PDG integer to the constructor, and m() is not automatically
+ // updated by the particle data. So we need to set it (to the pole mass):
+ double m0 = prts[0].m0();
+ prts[0].m( m0 );
+ hme->initChannel( prts );
return hme->decayWeightMax( prts );
}
+
int ih( 0 );
std::vector ihs( 3, -1 ), nhs( 3, 0 ), hs( mom->getNDaug() + 1, 0 );
// Translate the mother to HelicityParticle
EvtVector4R p = mom->getP4Lab();
prts.push_back( Pythia8::HelicityParticle(
mom->getPDGId(), 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;
}
// Translate the children to HelicityParticle
for ( int iDtr = 0; iDtr < (int)mom->getNDaug(); ++iDtr ) {
EvtParticle* dtr = mom->getDaug( iDtr );
p = dtr->getP4Lab();
prts.push_back( Pythia8::HelicityParticle(
dtr->getPDGId(), 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 ( 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;
}
}
if ( dtr->getSpinType() == EvtSpinType::NEUTRINO && dtr->getPDGId() < 0 )
hs[iDtr + 1] = 1;
}
// Calculate and set the amplitudes (maximum three helicities)
hme->initChannel( prts );
// The initWaves function in pythia8/include/Pythia8/HelicityMatrixElements.h
// needs to be made public (from protected) before we can use it:
// hme->initWaves( prts );
if ( ihs[0] < 0 ) {
amp->vertex( hme->calculateME( hs ) );
} else {
for ( hs[ihs[0]] = 0; hs[ihs[0]] < nhs[0]; ++hs[ihs[0]] ) {
if ( ihs[1] < 0 ) {
amp->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 ) {
amp->vertex( hs[ihs[0]], hs[ihs[1]],
hme->calculateME( hs ) );
} else {
for ( hs[ihs[2]] = 0; hs[ihs[2]] < nhs[2]; ++hs[ihs[2]] ) {
amp->vertex( hs[ihs[0]], hs[ihs[1]], hs[ihs[2]],
hme->calculateME( hs ) );
}
} // End third helicity
}
} // End second helicity
}
} // End first helicity
// Return the decay weight
return hme->decayWeight( prts );
}
#endif
diff --git a/test/jsonFiles/HmePythiaTau1.json b/test/jsonFiles/HmePythiaTau1.json
new file mode 100644
index 0000000..fb4c71d
--- /dev/null
+++ b/test/jsonFiles/HmePythiaTau1.json
@@ -0,0 +1,15 @@
+{
+ "parent" : "tau-",
+ "daughters" : ["e-", "anti-nu_e", "nu_tau"],
+ "model" : "HMEPYTHIA",
+ "parameters" : [1531],
+ "events" : 10000,
+ "histograms" : [
+ {"variable" : "mass", "title" : "mENuE" , "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+ {"variable" : "psumsq", "title" : "qSq", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+ {"variable" : "coshel", "title" : "cosE", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+ {"variable" : "prob", "title" : "Prob", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0}
+ ],
+ "outfile" : "HmePythiaTau1.root",
+ "reference" : "RefHmePythiaTau1.root"
+}
diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
index 817e8c0..e968c8c 100644
--- a/test/testDecayModel.cc
+++ b/test/testDecayModel.cc
@@ -1,357 +1,348 @@
/***********************************************************************
* 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 . *
***********************************************************************/
#include "testDecayModel.hh"
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtAbsRadCorr.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtMTRandomEngine.hh"
+#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
-#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4R.hh"
#ifdef EVTGEN_EXTERNAL
#include "EvtGenExternal/EvtExternalGenList.hh"
#endif
#include
#include
-#include
#include
+#include
#include
testDecayModel::testDecayModel( const std::string& jsonFile ) :
- m_jsonFile( jsonFile ),
- m_infoVect(),
- m_histVect(),
- m_nHistos( 0 )
+ m_jsonFile( jsonFile ), m_infoVect(), m_histVect(), m_nHistos( 0 )
{
}
-void testDecayModel::run() {
-
+void testDecayModel::run()
+{
Json::Value config;
std::ifstream inputStr( m_jsonFile.c_str() );
inputStr >> config;
inputStr.close();
const std::string parentName = config["parent"].asString();
Json::Value jDaughters = config["daughters"];
std::vector daugNames;
for ( auto daugN : jDaughters ) {
- daugNames.push_back( daugN.asString() );
+ daugNames.push_back( daugN.asString() );
}
const std::string modelName = config["model"].asString();
Json::Value jParameters = config["parameters"];
std::vector modelPars;
for ( auto par : jParameters ) {
- modelPars.push_back( par.asDouble() );
+ modelPars.push_back( par.asDouble() );
}
// Create decay file based on json file input
- const std::string decFile = createDecFile( parentName, daugNames,
- modelName, modelPars );
+ const std::string decFile = createDecFile( parentName, daugNames, modelName,
+ modelPars );
int nEvents = config["events"].asInt();
std::cout << "Number of events = " << nEvents << std::endl;
// ROOT output file
const std::string outFileName = config["outfile"].asString();
TFile* outFile = TFile::Open( outFileName.c_str(), "recreate" );
// Define histograms
defineHistos( config, outFile );
// Generate events and fill histograms
generateEvents( decFile, parentName, nEvents );
// Normalize histograms
for ( auto hist : m_histVect ) {
- double area = hist->Integral();
- if ( area > 0.0 ) { hist->Scale( 1.0/area ); }
+ double area = hist->Integral();
+ if ( area > 0.0 ) {
+ hist->Scale( 1.0 / area );
+ }
}
-
+
// Compare with reference histograms
const std::string refFileName = config["reference"].asString();
compareHistos( refFileName );
-
+
// Write output
outFile->cd();
for ( auto hist : m_histVect ) {
- hist->Write();
+ hist->Write();
}
outFile->Close();
-
}
-std::string testDecayModel::createDecFile( const std::string& parent,
- const std::vector& daugNames,
- const std::string& modelName,
- const std::vector& parameters ) const
+std::string testDecayModel::createDecFile(
+ const std::string& parent, const std::vector& daugNames,
+ const std::string& modelName, const std::vector& parameters ) const
{
-
// Create (or overwrite) the decay file
- std::string decName( modelName ); decName += ".dec";
+ std::string decName( modelName );
+ decName += ".dec";
std::ofstream decFile;
decFile.open( decName.c_str() );
decFile << "Decay " << parent << std::endl;
decFile << "1.0";
for ( auto daug : daugNames ) {
- decFile << " " << daug;
+ decFile << " " << daug;
}
decFile << " " << modelName;
for ( auto par : parameters ) {
- decFile << " " << par;
+ decFile << " " << par;
}
decFile << ";" << std::endl;
decFile << "Enddecay" << std::endl;
decFile << "End" << std::endl;
decFile.close();
-
+
return decName;
-
}
-void testDecayModel::defineHistos( Json::Value& config, TFile* outFile ) {
-
+void testDecayModel::defineHistos( Json::Value& config, TFile* outFile )
+{
// Histogram information
Json::Value jHistos = config["histograms"];
size_t nHistos = jHistos.size();
-
+
m_infoVect.reserve( nHistos );
m_histVect.reserve( nHistos );
-
- for ( auto hInfo : jHistos ) {
- const std::string varName = hInfo["variable"].asString();
- // Integer values that define what particles need to be used
- // for invariant mass combinations or helicity angles etc
- int d1 = hInfo["d1"].asInt();
- int d2 = hInfo["d2"].asInt();
-
- testInfo info( varName, d1, d2 );
- m_infoVect.push_back( info );
-
- const std::string varTitle = hInfo["title"].asString();
- int nBins = hInfo["nbins"].asInt();
- double xmin = hInfo["xmin"].asDouble();
- double xmax = hInfo["xmax"].asDouble();
-
- TString histName( varName.c_str() );
- if ( d1 != 0 ) { histName += "_"; histName += d1; }
- if ( d2 != 0 ) { histName += "_"; histName += d2; }
- TH1D* hist = new TH1D( histName.Data(), varTitle.c_str(), nBins, xmin, xmax );
- hist->SetDirectory( outFile );
- m_histVect.push_back( hist );
-
- }
+ for ( auto hInfo : jHistos ) {
+ const std::string varName = hInfo["variable"].asString();
+ // Integer values that define what particles need to be used
+ // for invariant mass combinations or helicity angles etc
+ int d1 = hInfo["d1"].asInt();
+ int d2 = hInfo["d2"].asInt();
+
+ testInfo info( varName, d1, d2 );
+ m_infoVect.push_back( info );
+
+ const std::string varTitle = hInfo["title"].asString();
+ int nBins = hInfo["nbins"].asInt();
+ double xmin = hInfo["xmin"].asDouble();
+ double xmax = hInfo["xmax"].asDouble();
+
+ TString histName( varName.c_str() );
+ if ( d1 != 0 ) {
+ histName += "_";
+ histName += d1;
+ }
+ if ( d2 != 0 ) {
+ histName += "_";
+ histName += d2;
+ }
+ TH1D* hist = new TH1D( histName.Data(), varTitle.c_str(), nBins, xmin,
+ xmax );
+ hist->SetDirectory( outFile );
+ m_histVect.push_back( hist );
+ }
m_nHistos = m_histVect.size();
-
}
-void testDecayModel::generateEvents( const std::string& decFile, const std::string& parentName,
- int nEvents ) {
-
+void testDecayModel::generateEvents( const std::string& decFile,
+ const std::string& parentName, int nEvents )
+{
// Define the random number generator
auto randomEngine = std::make_unique();
EvtAbsRadCorr* radCorrEngine = nullptr;
std::list extraModels;
#ifdef EVTGEN_EXTERNAL
bool convertPythiaCodes( false );
bool useEvtGenRandom( true );
EvtExternalGenList genList( convertPythiaCodes, "", "gamma", useEvtGenRandom );
radCorrEngine = genList.getPhotosModel();
extraModels = genList.getListOfModels();
#endif
-
+
EvtGen theGen( "../DECAY.DEC", "../evt.pdl", randomEngine.get(),
- radCorrEngine, &extraModels );
+ radCorrEngine, &extraModels );
theGen.readUDecay( decFile.c_str() );
// Generate the decays
EvtId parId = EvtPDL::getId( parentName.c_str() );
for ( int i = 0; i < nEvents; i++ ) {
+ if ( i % 1000 == 0 ) {
+ std::cout << "Event " << nEvents - i << std::endl;
+ }
+
+ // Initial 4-momentum and particle
+ EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 );
+ EvtParticle* parent = EvtParticleFactory::particleFactory( parId, pInit );
+
+ theGen.generateDecay( parent );
- if ( i%1000 == 0 ) { std::cout << "Event " << nEvents - i << std::endl; }
-
- // Initial 4-momentum and particle
- EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 );
- EvtParticle* parent = EvtParticleFactory::particleFactory( parId, pInit );
-
- theGen.generateDecay( parent );
-
- // Store information
- for ( size_t i = 0; i < m_nHistos; i++ ) {
-
- testInfo info = m_infoVect[i];
- double value = getValue( parent, info.getName(), info.getd1(), info.getd2() );
-
- TH1D* hist = m_histVect[i];
- if ( hist ) {
- hist->Fill( value );
- }
-
- }
-
- //parent->printTree();
-
- // Cleanup
- parent->deleteTree();
-
+ // Store information
+ for ( size_t i = 0; i < m_nHistos; i++ ) {
+ testInfo info = m_infoVect[i];
+ double value = getValue( parent, info.getName(), info.getd1(),
+ info.getd2() );
+
+ TH1D* hist = m_histVect[i];
+ if ( hist ) {
+ hist->Fill( value );
+ }
+ }
+
+ //parent->printTree();
+
+ // Cleanup
+ parent->deleteTree();
}
-
}
double testDecayModel::getValue( EvtParticle* parent, const std::string& varName,
- const int d1, const int d2 ) const {
-
+ const int d1, const int d2 ) const
+{
double value( 0.0 );
- if ( !parent ) { return value; }
+ if ( !parent ) {
+ return value;
+ }
// Daughters
- const EvtParticle* par1 = d1 > 0 ? parent->getDaug( d1 - 1 ) : nullptr;
- const EvtParticle* par2 = d2 > 0 ? parent->getDaug( d2 - 1 ) : nullptr;
+ const int NDaugMax( parent->getNDaug() );
+ const EvtParticle* par1 = d1 > 0 && d1 <= NDaugMax ? parent->getDaug( d1 - 1 )
+ : nullptr;
+ const EvtParticle* par2 = d2 > 0 && d2 <= NDaugMax ? parent->getDaug( d2 - 1 )
+ : nullptr;
// 4-momenta (in parent rest frame)
const EvtVector4R p1 = par1 != nullptr ? par1->getP4() : EvtVector4R();
const EvtVector4R p2 = par2 != nullptr ? par2->getP4() : EvtVector4R();
if ( !varName.compare( "mass" ) ) {
+ // Invariant mass
+ if ( d2 != 0 ) {
+ // Invariant 4-mass combination of particles d1 and d2
+ value = ( p1 + p2 ).mass();
+ } else {
+ // Invariant mass of particle d1 only
+ value = p1.mass();
+ }
- // Invariant mass
- if ( d2 != 0 ) {
- // Invariant 4-mass combination of particles d1 and d2
- value = ( p1 + p2 ).mass();
- } else {
- // Invariant mass of particle d1 only
- value = p1.mass();
- }
-
} else if ( !varName.compare( "psumsq" ) ) {
-
- // Invariant momentum sum squared
- value = ( p1 + p2 ).mass2();
+ // Invariant momentum sum squared
+ value = ( p1 + p2 ).mass2();
} else if ( !varName.compare( "pdiffsq" ) ) {
-
- // Invariant momentum difference squared
- value = ( p1 - p2 ).mass2();
+ // Invariant momentum difference squared
+ value = ( p1 - p2 ).mass2();
} else if ( !varName.compare( "mtm" ) ) {
+ // Momentum of particle d1
+ value = p1.d3mag();
- // Momentum of particle d1
- value = p1.d3mag();
-
} else if ( !varName.compare( "coshel" ) ) {
-
- // Cosine of helicity angle
-
- // Resonance center-of-mass system (d1 and d2)
- const EvtVector4R p12 = p1 + p2;
- // Boost vector
- const EvtVector4R boost( p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ), -p12.get( 3 ) );
- // Momentum of particle d1 in resonance frame, p1Res
- const EvtVector4R p1Res = boostTo( p1, boost );
- // Cosine of angle between p1Res and momentum of resonance in parent frame
- double p1ResMag = p1Res.d3mag();
- double p12Mag = p12.d3mag();
- if ( p1ResMag > 0.0 && p12Mag > 0.0 ) {
- value = p1Res.dot( p12 ) / ( p1ResMag * p12Mag );
- }
+ // Cosine of helicity angle
+
+ // Resonance center-of-mass system (d1 and d2)
+ const EvtVector4R p12 = p1 + p2;
+ // Boost vector
+ const EvtVector4R boost( p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ),
+ -p12.get( 3 ) );
+ // Momentum of particle d1 in resonance frame, p1Res
+ const EvtVector4R p1Res = boostTo( p1, boost );
+ // Cosine of angle between p1Res and momentum of resonance in parent frame
+ double p1ResMag = p1Res.d3mag();
+ double p12Mag = p12.d3mag();
+ if ( p1ResMag > 0.0 && p12Mag > 0.0 ) {
+ value = p1Res.dot( p12 ) / ( p1ResMag * p12Mag );
+ }
} else if ( !varName.compare( "prob" ) ) {
-
- // Decay probability
- double* dProb = parent->decayProb();
+ // Decay probability
+ double* dProb = parent->decayProb();
if ( dProb ) {
- value = *dProb;
+ value = *dProb;
}
-
}
-
+
return value;
-
}
-void testDecayModel::compareHistos( const std::string& refFileName ) const {
-
+void testDecayModel::compareHistos( const std::string& refFileName ) const
+{
// Compare histograms with the same name, calculating the chi-squared
TFile* refFile = TFile::Open( refFileName.c_str(), "read" );
if ( !refFile ) {
- std::cout << "Could not open reference file " << refFileName << std::endl;
- return;
+ std::cout << "Could not open reference file " << refFileName << std::endl;
+ return;
}
-
- for ( auto hist : m_histVect ) {
-
- const std::string histName = hist->GetName();
- // Get equivalent reference histogram
- const TH1* refHist = dynamic_cast( refFile->Get ( histName.c_str() ) );
-
- if ( refHist ) {
-
- double chiSq( 0.0 );
- int nDof( 0 );
- int iGood( 0 );
- double pValue = refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" );
- std::cout << "Histogram " << histName << " chiSq/nDof = "
- << chiSq << "/" << nDof << ", pValue = " << pValue << std::endl;
- } else {
-
- std::cout << "Could not find reference histogram " << histName << std::endl;
-
- }
-
+ for ( auto hist : m_histVect ) {
+ const std::string histName = hist->GetName();
+ // Get equivalent reference histogram
+ const TH1* refHist = dynamic_cast(
+ refFile->Get( histName.c_str() ) );
+
+ if ( refHist ) {
+ double chiSq( 0.0 );
+ int nDof( 0 );
+ int iGood( 0 );
+ double pValue = refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" );
+ std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq
+ << "/" << nDof << ", pValue = " << pValue << std::endl;
+
+ } else {
+ std::cout << "Could not find reference histogram " << histName
+ << std::endl;
+ }
}
refFile->Close();
-
}
-int main( int argc, char* argv[] ) {
-
+int main( int argc, char* argv[] )
+{
if ( argc != 2 ) {
- std::cout << "Expecting one argument: json input file" << std::endl;
- return -1;
+ std::cout << "Expecting one argument: json input file" << std::endl;
+ return -1;
}
const std::string jsonFile = argv[1];
testDecayModel test( jsonFile );
test.run();
return 0;
-
}
diff --git a/test/testDecayModel.hh b/test/testDecayModel.hh
index 10b2196..5fddd8c 100644
--- a/test/testDecayModel.hh
+++ b/test/testDecayModel.hh
@@ -1,87 +1,80 @@
/***********************************************************************
* 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 . *
***********************************************************************/
#ifndef TEST_DECAY_MODEL_HH
#define TEST_DECAY_MODEL_HH
-#include "json/json.h"
-
#include "TFile.h"
#include "TH1D.h"
+#include "json/json.h"
+
#include
#include
class EvtParticle;
class testInfo {
-
public:
testInfo( const std::string& name, int d1, int d2 ) :
- m_name( name ),
- m_d1( d1 ),
- m_d2( d2 )
- {;}
-
+ m_name( name ), m_d1( d1 ), m_d2( d2 )
+ {
+ ;
+ }
+
const std::string getName() const { return m_name; }
int getd1() const { return m_d1; }
int getd2() const { return m_d2; }
-
+
private:
std::string m_name;
int m_d1;
int m_d2;
-
};
class testDecayModel {
-
public:
-
testDecayModel( const std::string& jsonFile );
- void run();
+ void run();
private:
-
std::string createDecFile( const std::string& parent,
- const std::vector& daugNames,
- const std::string& model,
- const std::vector& parameters ) const;
+ const std::vector& daugNames,
+ const std::string& model,
+ const std::vector& parameters ) const;
void defineHistos( Json::Value& config, TFile* theFile );
- void generateEvents( const std::string& decFile, const std::string& parentName,
- int nEvents );
+ void generateEvents( const std::string& decFile,
+ const std::string& parentName, int nEvents );
double getValue( EvtParticle* rootPart, const std::string& varName,
- const int d1, const int d2 ) const;
-
+ const int d1, const int d2 ) const;
+
void compareHistos( const std::string& refFileName ) const;
-
+
std::string m_jsonFile;
std::vector m_infoVect;
std::vector m_histVect;
size_t m_nHistos;
-
-
};
#endif