No OneTemporary
Actions

Size

192 KB

Subscribers

None

View Options

	diff --git a/.clang-format b/.clang-format
	index 5add6cd..3d32fc3 100644
	--- a/.clang-format
	+++ b/.clang-format
	@@ -1,170 +1,170 @@
	Language: Cpp
	AccessModifierOffset: -2
	AlignAfterOpenBracket: Align
	AlignArrayOfStructures: Right
	AlignConsecutiveAssignments: None
	AlignConsecutiveBitFields: Consecutive
	AlignConsecutiveDeclarations: None
	AlignEscapedNewlines: Left
	AlignOperands: AlignAfterOperator
	AlignTrailingComments:
	Kind: Always
	OverEmptyLines: 0
	AllowAllArgumentsOnNextLine: false
	AllowAllParametersOfDeclarationOnNextLine: false
	AllowShortBlocksOnASingleLine: Never
	AllowShortCaseLabelsOnASingleLine: false
	AllowShortEnumsOnASingleLine: false
	AllowShortFunctionsOnASingleLine: InlineOnly
	AllowShortIfStatementsOnASingleLine: Never
	AllowShortLambdasOnASingleLine: All
	AllowShortLoopsOnASingleLine: false
	AlwaysBreakAfterDefinitionReturnType: None
	AlwaysBreakAfterReturnType: None
	AlwaysBreakBeforeMultilineStrings: false
	AlwaysBreakTemplateDeclarations: Yes
	BinPackArguments: true
	BinPackParameters: true
	BitFieldColonSpacing: Both
	BraceWrapping:
	AfterCaseLabel: false
	AfterClass: false
	AfterControlStatement: Never
	AfterEnum: true
	AfterFunction: true
	AfterNamespace: false
	AfterStruct: false
	AfterUnion: false
	AfterExternBlock: false
	BeforeCatch: false
	BeforeElse: false
	BeforeLambdaBody: true
	BeforeWhile: false
	IndentBraces: false
	SplitEmptyFunction: true
	SplitEmptyRecord: true
	SplitEmptyNamespace: true
	BreakAfterAttributes: Leave
	BreakArrays: false
	BreakBeforeBinaryOperators: None
	BreakBeforeBraces: Custom
	BreakBeforeConceptDeclarations: Always
	BreakBeforeTernaryOperators: true
	BreakConstructorInitializers: AfterColon
	BreakInheritanceList: AfterColon
	BreakStringLiterals: false
	ColumnLimit: 80
	CommentPragmas: '^ IWYU pragma:'
	CompactNamespaces: false
	ConstructorInitializerIndentWidth: 4
	ContinuationIndentWidth: 4
	Cpp11BracedListStyle: true
	DerivePointerAlignment: false
	DisableFormat: false
	EmptyLineAfterAccessModifier: Never
	EmptyLineBeforeAccessModifier: Always
	ExperimentalAutoDetectBinPacking: false
	FixNamespaceComments: true
	ForEachMacros:
	- foreach
	- BOOST_FOREACH
	IncludeBlocks: Regroup
	IncludeCategories:
	- Regex: '^"EvtGen/'
	Priority: 1
	- Regex: '^"EvtGenBase/'
	Priority: 2
	- Regex: '^"EvtGenModels/'
	Priority: 3
	- Regex: '^"EvtGenExternal/'
	Priority: 4
	- Regex: '^"(HepMC\|HepMC3)/'
	Priority: 5
	- - Regex: '^"(Pythia8\|Photos\|Tauola\|SHERPA\|ATOOLS)/'
	+ - Regex: '^"(Pythia8\|Photos\|Tauola\|SHERPA\|ATOOLS\|PHOTONS\+\+)/'
	Priority: 6
	- Regex: '^"T[[:alnum:]]+\.h"'
	Priority: 7
	- Regex: '^<[[:alnum:]_]+>'
	Priority: 9
	- Regex: '^<[[:alnum:]_]+\.h>'
	Priority: 10
	- Regex: '.*'
	Priority: 8
	IncludeIsMainRegex: '(Test)?$'
	IncludeIsMainSourceRegex: '(Impl\.hh)$'
	IndentAccessModifiers: false
	IndentCaseBlocks: false
	IndentCaseLabels: true
	IndentExternBlock: Indent
	IndentGotoLabels: true
	IndentPPDirectives: None
	IndentRequiresClause: false
	IndentWidth: 4
	IndentWrappedFunctionNames: false
	InsertBraces: false
	InsertNewlineAtEOF: true
	IntegerLiteralSeparator:
	Binary: 0
	Decimal: 0
	Hex: 0
	KeepEmptyLinesAtTheStartOfBlocks: false
	LambdaBodyIndentation: Signature
	LineEnding: LF
	MacroBlockBegin: ''
	MacroBlockEnd: ''
	MaxEmptyLinesToKeep: 1
	NamespaceIndentation: All
	PPIndentWidth: 0
	PackConstructorInitializers: NextLine
	PenaltyBreakAssignment: 50
	PenaltyBreakBeforeFirstCallParameter: 50
	PenaltyBreakComment: 250
	PenaltyBreakFirstLessLess: 100
	PenaltyBreakOpenParenthesis: 0
	PenaltyBreakString: 1000
	PenaltyBreakTemplateDeclaration: 1
	PenaltyExcessCharacter: 10
	PenaltyIndentedWhitespace: 0
	PenaltyReturnTypeOnItsOwnLine: 250
	PointerAlignment: Left
	QualifierAlignment: Left
	ReferenceAlignment: Pointer
	ReflowComments: false
	RemoveBracesLLVM: false
	RemoveSemicolon: false
	RequiresClausePosition: OwnLine
	RequiresExpressionIndentation: OuterScope
	SeparateDefinitionBlocks: Leave
	ShortNamespaceLines: 0
	SortIncludes: CaseSensitive
	SortUsingDeclarations: Lexicographic
	SpaceAfterCStyleCast: false
	SpaceAfterLogicalNot: false
	SpaceAfterTemplateKeyword: true
	SpaceAroundPointerQualifiers: Default
	SpaceBeforeAssignmentOperators: true
	SpaceBeforeCaseColon: false
	SpaceBeforeCpp11BracedList: false
	SpaceBeforeCtorInitializerColon: true
	SpaceBeforeInheritanceColon: true
	SpaceBeforeParens: ControlStatements
	SpaceBeforeRangeBasedForLoopColon: true
	SpaceBeforeSquareBrackets: false
	SpaceInEmptyBlock: false
	SpaceInEmptyParentheses: false
	SpacesBeforeTrailingComments: 4
	SpacesInAngles: Never
	SpacesInCStyleCastParentheses: false
	SpacesInConditionalStatement: true
	SpacesInContainerLiterals: true
	SpacesInLineCommentPrefix:
	Minimum: 1
	Maximum: 1
	SpacesInParentheses: true
	SpacesInSquareBrackets: false
	Standard: Latest
	StatementMacros:
	- ClassDef
	- ClassImp
	TabWidth: 4
	UseTab: Never
	diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
	index 75bc784..8dfeb14 100644
	--- a/.gitlab-ci.yml
	+++ b/.gitlab-ci.yml
	@@ -1,191 +1,193 @@
	stages:
	- checks
	- build
	- test

	variables:
	- LCG_VERSION: "LCG_105"
	- GCC_VERSION: "gcc13"
	+ LCG_VERSION: "LCG_106a"
	+ GCC_VERSION: "gcc14"
	CLANG_VERSION: "clang16"
	BUILD_TYPE: "Release"
	BUILD_TESTS: "OFF"
	BUILD_DOXYGEN: "OFF"
	USE_HEPMC3: "ON"
	USE_PYTHIA: "ON"
	USE_PHOTOS: "ON"
	USE_TAUOLA: "ON"
	USE_SHERPA: "ON"
	+ USE_SHERPA3: "ON"

	.production_os_image:
	variables:
	LCG_OS: x86_64-el9
	image: gitlab-registry.cern.ch/lhcb-docker/os-base/alma9-devel:latest
	tags:
	- k8s-cvmfs

	.legacy_os_image:
	variables:
	LCG_OS: x86_64-centos7
	LCG_VERSION: "LCG_103"
	GCC_VERSION: "gcc12"
	CLANG_VERSION: "clang12"
	+ USE_SHERPA3: "OFF"
	image: gitlab-registry.cern.ch/lhcb-docker/os-base/centos7-devel:latest
	tags:
	- k8s-cvmfs

	.lcg_setup:
	before_script:
	- set +e && source /cvmfs/sft.cern.ch/lcg/views/setupViews.sh $LCG_VERSION $LCG_OS-$LCG_COMPILER; set -e

	.check_template:
	stage: checks
	variables:
	LCG_COMPILER: "$CLANG_VERSION-opt"
	extends:
	- .production_os_image
	- .lcg_setup

	.build_template:
	stage: build
	extends:
	- .lcg_setup
	script:
	- mkdir install
	- mkdir build && cd build
	- - cmake -DCMAKE_BUILD_TYPE:STRING=$BUILD_TYPE -DCMAKE_INSTALL_PREFIX:PATH=$CI_PROJECT_DIR/install -DEVTGEN_HEPMC3:BOOL=$USE_HEPMC3 -DEVTGEN_PYTHIA:BOOL=$USE_PYTHIA -DEVTGEN_PHOTOS:BOOL=$USE_PHOTOS -DEVTGEN_TAUOLA:BOOL=$USE_TAUOLA -DEVTGEN_SHERPA:BOOL=$USE_SHERPA -DEVTGEN_BUILD_TESTS:BOOL=$BUILD_TESTS -DEVTGEN_BUILD_VALIDATIONS:BOOL=$BUILD_TESTS -DEVTGEN_BUILD_DOXYGEN:BOOL=$BUILD_DOXYGEN $CI_PROJECT_DIR
	+ - cmake -DCMAKE_BUILD_TYPE:STRING=$BUILD_TYPE -DCMAKE_INSTALL_PREFIX:PATH=$CI_PROJECT_DIR/install -DEVTGEN_HEPMC3:BOOL=$USE_HEPMC3 -DEVTGEN_PYTHIA:BOOL=$USE_PYTHIA -DEVTGEN_PHOTOS:BOOL=$USE_PHOTOS -DEVTGEN_TAUOLA:BOOL=$USE_TAUOLA -DEVTGEN_SHERPA:BOOL=$USE_SHERPA -DEVTGEN_SHERPA3:BOOL=$USE_SHERPA3 -DEVTGEN_BUILD_TESTS:BOOL=$BUILD_TESTS -DEVTGEN_BUILD_VALIDATIONS:BOOL=$BUILD_TESTS -DEVTGEN_BUILD_DOXYGEN:BOOL=$BUILD_DOXYGEN $CI_PROJECT_DIR
	- cmake --build .
	- cmake --build . --target install

	check_formatting:
	extends:
	- .check_template
	script:
	- source applyFormatting.sh
	- git status --porcelain \| grep '^ M' \|\| exit 0
	- echo "The above files do not comply with formatting requirements"
	- exit 1

	check_test_json_files:
	extends:
	- .check_template
	script:
	- cd test
	- python checkJsonFiles.py

	check_member_varnames:
	extends:
	- .check_template
	script:
	- mkdir build && cd build
	- cmake -DEVTGEN_HEPMC3:BOOL=ON -DEVTGEN_PYTHIA:BOOL=ON -DEVTGEN_PHOTOS:BOOL=ON -DEVTGEN_TAUOLA:BOOL=ON -DEVTGEN_SHERPA:BOOL=ON -DEVTGEN_BUILD_TESTS:BOOL=OFF -DEVTGEN_BUILD_DOXYGEN:BOOL=ON $CI_PROJECT_DIR
	- cmake --build . --target doc
	- curl -o renameMemberVars.py "https://gitlab.cern.ch/evtgen/utilities/raw/master/RenameMemberVars/renameMemberVars.py?inline=false"
	- python renameMemberVars.py check $CI_PROJECT_DIR doxygen/html

	build_clang_opt:
	variables:
	LCG_COMPILER: "$CLANG_VERSION-opt"
	extends:
	- .production_os_image
	- .build_template

	build_gcc_opt:
	variables:
	LCG_COMPILER: "$GCC_VERSION-opt"
	extends:
	- .production_os_image
	- .build_template

	build_hepmc2_gcc_opt:
	variables:
	LCG_VERSION: "LCG_103_LHCB_7"
	LCG_COMPILER: "$GCC_VERSION-opt"
	USE_HEPMC3: "OFF"
	USE_SHERPA: "OFF"
	extends:
	- .legacy_os_image
	- .build_template

	build_gcc_dbg:
	variables:
	LCG_COMPILER: "$GCC_VERSION-dbg"
	BUILD_TYPE: "Debug"
	BUILD_TESTS: "ON"
	extends:
	- .production_os_image
	- .build_template
	after_script:
	- cd install/share/EvtGen/test
	- ./makeModelsJson $CI_PROJECT_DIR > Models.json
	- ./makeSrcDepsJson $CI_PROJECT_DIR $CI_PROJECT_DIR/build > SrcDeps.json
	artifacts:
	paths:
	- install
	expire_in: 1 day
	when: always

	build_legacy_clang_opt:
	variables:
	LCG_COMPILER: "$CLANG_VERSION-opt"
	extends:
	- .legacy_os_image
	- .build_template

	build_legacy_gcc_opt:
	variables:
	LCG_COMPILER: "$GCC_VERSION-opt"
	extends:
	- .legacy_os_image
	- .build_template

	.test_template:
	stage: test
	variables:
	LCG_COMPILER: "$GCC_VERSION-dbg"
	extends:
	- .production_os_image
	- .lcg_setup
	needs:
	- job: build_gcc_dbg
	artifacts: true
	script:
	- cd install/share/EvtGen/test
	- test ${#CHANGED_FILES} -eq 0 && echo "No files have changed, no need to run tests" && exit 0
	- ./runTests.py Models.json SrcDeps.json $CHANGED_FILES
	artifacts:
	paths:
	- install/share/EvtGen/test
	expire_in: 1 day
	when: always

	# TODO - need to copy in the appropriate files to act as references for this run
	# and save the resulting files for possible future use as references
	# - where to store the reference files? CI cache? EOS?
	# - need to add a test case for merge requests

	test_on_push_to_master:
	extends:
	- .test_template
	rules:
	- if: '$CI_PIPELINE_SOURCE == "push" && $CI_COMMIT_BRANCH == "master"'
	before_script:
	- !reference [.lcg_setup, before_script]
	- CHANGED_FILES=`git diff --name-only $CI_COMMIT_BEFORE_SHA $CI_COMMIT_SHA`

	test_on_push_to_branch:
	extends:
	- .test_template
	rules:
	- if: '$CI_PIPELINE_SOURCE == "push" && $CI_COMMIT_BRANCH != "master"'
	before_script:
	- !reference [.lcg_setup, before_script]
	- CHANGED_FILES=`git diff --name-only origin/master...$CI_COMMIT_SHA`

	test_on_schedule:
	extends:
	- .test_template
	rules:
	- if: '$CI_PIPELINE_SOURCE == "schedule" && $CI_COMMIT_BRANCH == "master"'
	before_script:
	- !reference [.lcg_setup, before_script]
	- CHANGED_FILES=`git tag --list 'R[[:digit:]][[:digit:]]-[[:digit:]][[:digit:]]-[[:digit:]][[:digit:]]' \| tail -n 1 \| xargs -I{} git diff --name-only {} $CI_COMMIT_SHA`
	diff --git a/CMakeLists.txt b/CMakeLists.txt
	index c25d1e9..2c3f66d 100644
	--- a/CMakeLists.txt
	+++ b/CMakeLists.txt
	@@ -1,166 +1,171 @@

	########################################################################
	# 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/>. #
	########################################################################

	# Top level CMakeLists.txt for EvtGen

	# Enforce an out-of-source build.
	# Should be the first action in the top level CMakeLists.txt
	if(${CMAKE_SOURCE_DIR} STREQUAL ${CMAKE_BINARY_DIR})
	message(STATUS "EvtGen requires an out-of-source build.")
	message(STATUS "Please remove the following files from ${CMAKE_BINARY_DIR}:")
	message(STATUS " CMakeCache.txt")
	message(STATUS " CMakeFiles")
	message(STATUS "Once these files are removed, create a separate directory")
	message(STATUS "and run CMake from there, pointing it to:")
	message(STATUS " ${CMAKE_SOURCE_DIR}")
	message(FATAL_ERROR "in-source build detected")
	endif()

	# Also require a minimum version of CMake
	cmake_minimum_required(VERSION 3.11.0)

	# Project setup
	if(${CMAKE_VERSION} VERSION_LESS 3.12.0)
	project(EvtGen
	VERSION 3.0.0
	DESCRIPTION "Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons."
	LANGUAGES C CXX
	)
	else()
	project(EvtGen
	VERSION 3.0.0
	DESCRIPTION "Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons."
	LANGUAGES C CXX
	HOMEPAGE_URL "https://evtgen.hepforge.org"
	)
	endif()

	# Prepend this project's custom module path(s) to CMAKE_MODULE_PATH
	set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake/Modules ${CMAKE_MODULE_PATH})
	# Include needed modules to perform any custom setup
	# Install paths
	include(GNUInstallDirs)
	# Compilation/linking flags and related settings
	include(EvtGenCompilerFlags)
	# EvtGen external dependencies
	option(EVTGEN_HEPMC3 "Enable/disable linking with HepMC3" ON)
	option(EVTGEN_PYTHIA "Enable/disable linking with Pythia8" OFF)
	option(EVTGEN_PHOTOS "Enable/disable linking with Photos++" OFF)
	option(EVTGEN_TAUOLA "Enable/disable linking with Tauola++" OFF)
	option(EVTGEN_SHERPA "Enable/disable linking with Sherpa" OFF)
	+option(EVTGEN_SHERPA3 "Prefer Sherpa 3 to Sherpa 2" ON)
	option(EVTGEN_SUPPRESS_EXTERNAL_WARNINGS "Enable/disable suppression of warnings from external dependencies" OFF)
	if (EVTGEN_HEPMC3)
	message(STATUS "EvtGen: Linking with HepMC3")
	else()
	message(STATUS "EvtGen: Linking with HepMC2")
	endif()
	message(STATUS "EvtGen: Optional linking with HepMC3 EVTGEN_HEPMC3 ${EVTGEN_HEPMC3}")
	message(STATUS "EvtGen: Optional linking with Pythia8 EVTGEN_PYTHIA ${EVTGEN_PYTHIA}")
	message(STATUS "EvtGen: Optional linking with Photos++ EVTGEN_PHOTOS ${EVTGEN_PHOTOS}")
	message(STATUS "EvtGen: Optional linking with Tauola++ EVTGEN_TAUOLA ${EVTGEN_TAUOLA}")
	-message(STATUS "EvtGen: Optional linking with Sherpa EVTGEN_SHERPA ${EVTGEN_SHERPA}")
	+if (EVTGEN_SHERPA3)
	+ message(STATUS "EvtGen: Optional linking with Sherpa3 EVTGEN_SHERPA ${EVTGEN_SHERPA}")
	+else()
	+ message(STATUS "EvtGen: Optional linking with Sherpa2 EVTGEN_SHERPA ${EVTGEN_SHERPA}")
	+endif()
	include(EvtGenExternalDependencies)

	# Setup clang-tidy checks
	option(EVTGEN_RUN_CLANG_TIDY "Enable/disable running clang-tidy static analysis checks during build" OFF)
	option(EVTGEN_CLANG_TIDY_FIX "Enable/disable automatic fixes in clang-tidy checks" OFF)
	set(EVTGEN_CLANG_TIDY_CHECKS "modernize-use-nullptr" CACHE STRING "Specify which clang-tidy checks to run (comma separated)")
	message(STATUS "EvtGen: Optional running of clang-tidy checks EVTGEN_RUN_CLANG_TIDY ${EVTGEN_RUN_CLANG_TIDY}")
	message(STATUS "EvtGen: Optional auto-fixes by clang-tidy checks EVTGEN_CLANG_TIDY_FIX ${EVTGEN_CLANG_TIDY_FIX}")
	message(STATUS "EvtGen: Optional list of clang-tidy checks to run EVTGEN_CLANG_TIDY_CHECKS ${EVTGEN_CLANG_TIDY_CHECKS}")
	include(EvtGenSetupClangTidyChecks)

	# Now build the library
	add_subdirectory(src)

	# Copy the particle property and decay tables
	configure_file(DECAY.DEC ${CMAKE_CURRENT_BINARY_DIR}/DECAY.DEC COPYONLY)
	configure_file(DECAY.XML ${CMAKE_CURRENT_BINARY_DIR}/DECAY.XML COPYONLY)
	configure_file(evt.pdl ${CMAKE_CURRENT_BINARY_DIR}/evt.pdl COPYONLY)
	# Copy the older decay tables
	configure_file(DECAY_2009.DEC ${CMAKE_CURRENT_BINARY_DIR}/DECAY_2009.DEC COPYONLY)
	configure_file(DECAY_2009.XML ${CMAKE_CURRENT_BINARY_DIR}/DECAY_2009.XML COPYONLY)
	configure_file(DECAY_2010.DEC ${CMAKE_CURRENT_BINARY_DIR}/DECAY_2010.DEC COPYONLY)
	configure_file(DECAY_2010.XML ${CMAKE_CURRENT_BINARY_DIR}/DECAY_2010.XML COPYONLY)

	# Build the executables in the test and validation directories
	option(EVTGEN_BUILD_DOC "Enable/disable building of userguide documentation in 'doc' directory" OFF)
	option(EVTGEN_BUILD_DOXYGEN "Enable/disable building of Doxygen documentation in 'doxygen' directory" OFF)
	option(EVTGEN_BUILD_TESTS "Enable/disable building of executables in 'test' directory" OFF)
	option(EVTGEN_BUILD_VALIDATIONS "Enable/disable building of executables in 'validation' directory" OFF)
	message(STATUS "EvtGen: Building of userguide documentation in 'doc' directory EVTGEN_BUILD_DOC ${EVTGEN_BUILD_DOC}")
	message(STATUS "EvtGen: Building of Doxygen documentation in 'doxygen' directory EVTGEN_BUILD_DOXYGEN ${EVTGEN_BUILD_DOXYGEN}")
	message(STATUS "EvtGen: Building of executables in 'test' directory EVTGEN_BUILD_TESTS ${EVTGEN_BUILD_TESTS}")
	message(STATUS "EvtGen: Building of executables in 'validation' directory EVTGEN_BUILD_VALIDATIONS ${EVTGEN_BUILD_VALIDATIONS}")
	if(EVTGEN_BUILD_DOC)
	add_subdirectory(doc)
	endif()
	if(EVTGEN_BUILD_DOXYGEN)
	add_subdirectory(doxygen)
	endif()
	if(EVTGEN_BUILD_TESTS)
	add_subdirectory(test)
	endif()
	if(EVTGEN_BUILD_VALIDATIONS)
	add_subdirectory(validation)
	endif()

	# Install the include directories
	install(DIRECTORY EvtGen DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
	install(DIRECTORY EvtGenBase DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
	install(DIRECTORY EvtGenExternal DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
	install(DIRECTORY EvtGenModels DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

	# Install the particle properties and decay tables
	install(FILES DECAY.DEC DECAY.XML evt.pdl DECAY_2010.XML DECAY_2010.DEC DECAY_2009.XML DECAY_2009.DEC DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME})

	# Install the COPYING and AUTHORS files
	install(FILES AUTHORS COPYING DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME})

	# Install the README and History files
	install(FILES README.md Pythia8_README.md Tauola_README.md History.md DESTINATION ${CMAKE_INSTALL_DOCDIR})

	# Generate CMake config files, which can be used by other projects
	include(CMakePackageConfigHelpers)

	set(INCLUDE_INSTALL_DIR ${CMAKE_INSTALL_INCLUDEDIR})
	set(LIB_INSTALL_DIR ${CMAKE_INSTALL_LIBDIR})
	set(DATA_INSTALL_DIR ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME})

	configure_package_config_file(cmake/Templates/EvtGenConfig.cmake.in
	${CMAKE_CURRENT_BINARY_DIR}/EvtGenConfig.cmake
	INSTALL_DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}/cmake
	PATH_VARS INCLUDE_INSTALL_DIR LIB_INSTALL_DIR DATA_INSTALL_DIR
	)

	write_basic_package_version_file(${CMAKE_CURRENT_BINARY_DIR}/EvtGenConfigVersion.cmake
	COMPATIBILITY AnyNewerVersion
	)

	install(FILES ${CMAKE_CURRENT_BINARY_DIR}/EvtGenConfig.cmake
	${CMAKE_CURRENT_BINARY_DIR}/EvtGenConfigVersion.cmake
	DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}/cmake
	)

	install(
	EXPORT "EvtGenTargets"
	NAMESPACE "EvtGen::"
	DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}/cmake
	)
	diff --git a/EvtGenExternal/EvtSherpaPhotons.hh b/EvtGenExternal/EvtSherpaPhotons.hh
	index 4c2db04..cf2dd0e 100644
	--- a/EvtGenExternal/EvtSherpaPhotons.hh
	+++ b/EvtGenExternal/EvtSherpaPhotons.hh
	@@ -1,87 +1,97 @@

	/***********************************************************************
	* 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_SHERPA
	#ifndef EVTSHERPAPHOTONS_HH
	#define EVTSHERPAPHOTONS_HH

	#include "EvtGenBase/EvtAbsRadCorr.hh"
	#include "EvtGenBase/EvtHepMCEvent.hh"
	#include "EvtGenBase/EvtId.hh"
	#include "EvtGenBase/EvtPDL.hh"
	#include "EvtGenBase/EvtParticle.hh"

	+#include "PHOTONS++/Main/Photons.H"
	#include "SHERPA/Main/Sherpa.H"

	#include <string>

	class EvtParticle;
	class EvtAbsExternalGen;

	/* Description: EvtGen's interface to PHOTOS for generation of
	* QED final state radiation.
	*/
	class EvtSherpaPhotons : public EvtAbsRadCorr {
	public:
	EvtSherpaPhotons( const bool useEvtGenRandom = true,
	const double infraredCutOff = 1.0e-7, const int mode = 2,
	const int useME = 0 );

	void initialise() override;

	void doRadCorr( EvtParticle* p ) override;

	+ // This function is needed to control the deletion of the static instances of Sherpa and PHOTONS++
	+ static void finalise();
	+
	private:
	// Provides vector of pointers to the configuration strings for Sherpa
	std::vector<char*> addParameters();

	// Updates the particle properties table of Sherpa
	void updateParticleLists();

	// Vector containing the configuration strings for Sherpa
	- // INIT_ONLY=6 intialises the Sherpa objects without launching simulation.
	- std::vector<std::string> m_configs{ "Sherpa", "INIT_ONLY=6" };
	+ std::vector<std::string> m_configs{ "Sherpa" };

	// Use EvtGen's random number generator
	bool m_useEvtGenRandom = true;

	// Default settings for PHOTONS++
	// Minimum photon Energy (infrared cut-off).
	double m_infraredCutOff = 1.0e-7;
	// Switches on (default) the hard emission corrections (mode 2),
	// otherwise follow the soft-photon only approach (mode 1).
	int m_mode = 2;
	// Switches off (default) or on the exact matrix-element corrections (works only for mode 2)
	int m_useME = 0;

	// Default photon properties
	const std::string m_photonType = "gamma";
	EvtId m_gammaId = EvtId( -1, -1 );
	long int m_gammaPDG = 22;

	- // The Sherpa instance.
	+ // Mutex to serialise all access to the shared Sherpa objects
	+ static std::mutex m_sherpa_mutex;
	+
	+ // The Sherpa instance
	static std::unique_ptr<SHERPA::Sherpa> m_sherpaGen;
	+
	+ // The PHOTONS++ instance
	+ static std::unique_ptr<PHOTONS::Photons> m_photonsGen;
	+
	+ // Flag whether initialisation has been performed
	static bool m_initialised;
	- static std::mutex m_sherpa_mutex;
	};

	#endif

	#endif
	diff --git a/History.md b/History.md
	index 1ac7611..2ba7dce 100644
	--- a/History.md
	+++ b/History.md
	@@ -1,1007 +1,1014 @@
	# History file for EvtGen

	From version 2.0.0,
	Tabc labels refer to [Maniphest tasks](https://phab.hepforge.org/maniphest/query/nkBRd9OhPCBN/), while
	Dxyz labels refer to [Differential code reviews](https://phab.hepforge.org/differential/query/YDY8EgjNGd.e/) on HepForge:

	https://phab.hepforge.org/Tabc

	https://phab.hepforge.org/Dxyz

	===
	## R03-00-00-betaX

	+3 Dec 2024 Fernando Abudinen
	+* D134: Add support for Sherpa 3
	+ - Modified EvtSherpaPhotons interface to support both Sherpa 2 and 3 releases
	+ - Introduced finalise function inside EvtSherpaPhotons class to control the order of static destruction
	+ - Updated CI to pick Sherpa 3, use Sherpa 2 for legacy builds
	+ - Updated setupEvtGen for Sherpa 3
	+
	2 Dec 2024 Thomas Latham
	* Remove or convert C headers to C++ equivalents

	18 Nov 2024 Thomas Latham
	* Update CI to use LCG\_105 on el9 with gcc13/clang16 by default
	- Also update CI to avoid spurious failure of test stage
	- Update clang-format config to clang16 and apply resulting changes

	25 Oct 2024 Thomas Latham
	* Fixes, updates and improvements to setupEvtGen script

	7 Oct 2024 Thomas Latham
	* D133: Fix bug in rotation to helicity basis for tensor particles

	4 Oct 2024 Heather Ratcliffe, Thomas Latham, Fernando Abudinen
	* D132: Introduce thread safety
	- Protect statics: make them const or thread\_local
	- Improve const-correctness of many functions
	- Make thread\_local the RNG, particle property table, external factory, decay tables, and sym tables
	- Add explicit mutexing to non-thread-safe externals (Sherpa, Tauola, Photos)
	- Propagate the EvtGen seed to the Fortran RNG in Tauola
	- Some fixes to EvtBtoXsEtap and EvtBtoXsgamma
	- Reorganise test code for multi-threading

	1 Oct 2024 Fernando Abudinen
	* D131: Tauola tests
	- Translated TAUOLA tests in validation to new json format
	- Covered all 22 supported TAUOLA models

	18 Sep 2024 Thomas Latham
	* D130: Fix recently introduced bug in EvtDalitzReso for parent with significant natural width
	- Add a test for eta\_c -> K+ K- pi0 under GENERIC\_DALITZ model
	- Modify code to restore correct behaviour
	- Make other general improvements to EvtDalitzReso class

	14 Aug 2024 Fernando Abudinen
	* D128: Improvements in EvtId and EvtPDL.
	- Modernisation of EvtId class improving Boolean functions.
	- Conversion of EvtPDL into thread-local singleton.
	- Propagation of changes in EvtPDL to EvtGen initialisation.

	9 Aug 2024 Fernando Abudinen
	* D129: Added tests for PYTHIA interface
	- Introduced new tests for decays with PYTHIA.
	- Converted preexisting tests to json format based on new testing framework.

	1 Aug 2024 Fernando Abudinen
	* D127: Improve FSR mutexing by using `std::lock_guard`.

	3 July 2024 Fernando Abudinen
	* D126: Fix FSR initialisation for thread safety.
	- Modified PHOTOS and Sherpa instances to become global static members
	and made sure that they are initialised only once.
	- Made sure that calls to these instances are mutexed.
	- Made sure that all interface class members are initialised for each thread.

	1 July 2024 Fernando Abudinen
	* D125: Removed `EventHandler` in `EvtSherpaPhotons` interface to reduce initialisation overhead.

	24 June 2024 John Back
	* D124: Add EvtD0ToKspipi DP model, courtesy of Camille Normand (LHCb).

	21 June 2024 Thomas Latham
	* D123: Removed `EvtPatches.hh`. Removed unused `clock_t` typedef and redundant `UNUSED` macro.

	20 June 2024 Thomas Latham
	* D122: Fix compilation warnings under gcc13 and clang16

	20 June 2024 Fernando Abudinen
	* D121: Introduced function `setSeed` to set seed for random number engines.

	19 June 2024 Thomas Latham
	* D120: Remove "reference" field from JSON test files.
	Instead simply prepend "Ref/" to output ROOT file name.

	22 May 2024 Fernando Abudinen
	* D118: Removed obsolete `ISR` attribute for photons. It was redundant as it was supposed to be
	the contrary of the `FSR` attribute. It was never used and could lead to confusion.
	- The `FSR` attribute is now set to 0 by default at initialisation for `EvtPhotonParticle`.

	21 May 2024 Fernando Abudinen
	* D115: Introduced Sherpa's PHOTONS++ generator
	as alternative for final-state radiation simulation.
	- Added Sherpa dependencies and linking in compilation files
	- Introduced new `EvtSherpaPhotons` interface
	- Introduced new `EvtSherpaRandom` interface to propagate EvtGen's random number into Sherpa
	- Implemented `fsrGenerator` as input argument in `testDecayModel` module

	8 May 2024 Michal Kreps, Thomas Latham, Andrii Verbytskyi
	* D58: Rename member variables to follow convention of `m_` prefix
	- In addition:
	* Introduce pre-build check in the CI that convention is being followed in all classes
	* Fix shadowing of variables in some classes
	* Fixes to Doxygen configuration
	* Enable shadowing and overloaded-virtual compiler warnings

	23 Apr 2024 Thomas Latham
	* D117: Drop redundant EVTGEN\_CPP11 preprocessor switch

	22 Apr 2024 Fernando Abudinen
	* D116: Replaced PHOTOS with FSR flags in examples and tests.

	17 Apr 2024 Fernando Abudinen
	* D114: Simplified pre-processor behaviour for PHOTOS.

	12 Apr 2024 Fernando Abudinen
	* D113: Introduced new FSR flag and deprecated PHOTOS flag.
	Implemented proper use of neverRadCorr.

	12 Apr 2024 Fernando Abudinen
	* D112: Removed EvtPhotosEngine and moved its functionality EvtPHOTOS class.
	Moved initialisation of fsrEngine to setRadCorrEngine.

	11 Apr 2024 Thomas Latham
	* D109: Remove broken or obsolete models
	- Models removed:
	* BHADRONIC
	* BQTOLLLLHYPERCP
	* BQTOLLLL
	* D\_MULTIBODY
	* KSTARSTARGAMMA
	* LAMBDAB2LAMBDAV and subsidiary models LAMBDA2PPIFORLAMBDAB2LAMBDAV and V2VPVMFORLAMBDAB2LAMBDAV
	* Lb2Lhh
	* PROPSLPOLE
	* VPHOTOV
	* VPHOTOVISR
	* VUB\_AC
	- Classes in EvtGenBase already unused or used only by the above models also removed
	* Evt3Rank3C
	* EvtAbsBinning
	* EvtAmpAmpPdf
	* EvtAmpIndex
	* EvtAmpSubIndex
	* EvtBreitWignerPdf
	* EvtDecayParm
	* EvtIntegPdf1D
	* EvtMBreitWigner
	* EvtMHelAmp
	* EvtMNode
	* EvtMParticle
	* EvtMLineShape
	* EvtMRes
	* EvtMTree
	* EvtMTrivialLS
	* EvtPointPred
	* EvtPto3PAmpSmpResolution
	* EvtSpinAmp

	9 Apr 2024 Fernando Abudinen
	* D111: Added tests for decays with FSR and implemented const correctness for getP4LabBeforeFSR() function.

	19 Feb 2024 Fernando Abudinen
	* D108: Bugfix protect polarisation vector for vector and tensor particles against wrong indices

	19 Feb 2024 Fernando Abudinen
	* D107: Implement const correctness for epsPhoton and epsParentPhoton

	3 Feb 2024 John Back
	* T230: Add EvtBcVPPHad model for Bc to Jpsi p pbar pi decays and generalised particle ordering
	for EvtBcVHad decay files, courtesy of Aleksei Luchinsky, Dmitrii Pereima & Vanya Belyaev (LHCb).
	Updated EvtPhiDalitz model to use helicity amplitudes and fixed the indices used in
	EvtVector3R::dot(), courtesy of Arnau Brossa Gonzalo & Antonio Romero Vidal (LHCb).

	31 Oct 2023 Fernando Abudinen
	* D102: Bugfix probmax issue for X38722-+\_PSI\_GAMMA model
	- Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.
	- Introduced weighting to compensate for phase-space changes due to resonance width.
	- Fixed bugs in rho couplings and loops over photon and vector states.

	16 Oct 2023 Fernando Abudinen
	* D99: Add tests for PHOTOS
	- EvtParticle::getAttribute function made const
	- Added variables for testing of FSR to testing framework
	- Added dedicated tests for FSR with Photos
	- Turned off Photos for all other tests

	16 Oct 2023 Thomas Latham
	* D98: Modernise EvtIdSet and other improvements
	- Modernise and greatly simplify EvtIdSet implementation
	- Fixes in EvtPropSLPole to avoid unnecessary dynamic allocations
	- Other minor fixes and tidy-ups
	* Credit to Heather Ratcliffe and Chris Brady for providing
	and/or inspiring several of these improvements

	22 Aug 2023 Fernando Abudinen
	* D97: Bugfix probmax issue and introduced pole compensation for VTOSLL model
	- Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.

	22 Aug 2023 Tom Latham
	* CMake updates
	- Update default C++ standard to 17
	- Suppress 'up-to-date' messages during install
	- Add 'EvtGen' prefix to file names of custom CMake modules
	- Remove unused CMake module
	* Update CI config
	- Update default LCG version to 103
	- Add builds for el9 OS using LCG 104
	- Allow switching on/off building against each external
	- Test stage: attempt to improve selection of commits to diff in different cases
	- Update CVMFS tags as per https://cern.service-now.com/service-portal?id=outage&n=OTG0079356
	* Fix script for generating src dependencies
	- Adapted to behaviour of newer CMake version in LCG 103
	* Apply clang-format
	- Few small changes after update of CI clang version to 12
	- Document version of clang-format to be used
	* Updates to install script
	- Update to latest Pythia8 and HepMC3 versions
	- Fix to ensure dependencies are picked up from the local install

	22 Aug 2023 Andrii Verbytskyi
	* Patch for finding Pythia8 xmldoc path

	21 Aug 2023 Tom Latham
	* D96: Work around change of interface for setRndmEnginePtr in Pythia8 310

	22 Jun 2023 Ludovico Massaccesi
	* T219: Fixed amplitudes for D_DALITZ model of D0 -> pi+ pi- pi0.

	1 March 2023 Fernando Abudinen
	* D92: Bugfix probmax issue for TENSOR daughter in EvtSSD_DirectCP model.
	Calculation of amplitude in EvtSSD_DirectCP model moved to dedicated calcAmp function.
	Got rid of a few static variables along the way.

	8 Feb 2023 Alexei Sibidanov and Tom Latham
	* D90: Apply clang format and enable checking of formatting in gitlab CI

	3 Feb 2023 John Back
	* D91: Check for non-zero momentum for EvtSLBaryonAmp parent spin density matrix.
	Print out integrals of JSON test histograms.

	16 Dec 2022 John Back
	* D89: Added probabilities for B_c -> V pi+ and V pi+ pi0 EvtBcVHad modes.

	16 Dec 2022 Alexei Sibidanov
	* D88: Applied clang-tidy modernize-use-nullptr

	13 Dec 2022 Fernando Abudinen and Tom Latham
	* Various improvements to testing framework

	16 Nov 2022 Tom Latham
	* T123: Provide documention of how to contribute bug reports, feature requests, new/modified code

	8 Sep 2022 Fernando Abudinen, John Back, Michal Kreps, Tom Latham, Alex Ward
	* T108: Implement JSON test framework for all decay models

	9 June 2022 Michal Kreps
	* D84: Improve efficiency of RareLbToLll decay model for final states with e+e- pair.

	===
	## R02-02-00

	12 May 2022 Michal Kreps
	* D83: Fix double counting of charmonia with K0 decays for B0.

	11 May 2022 Tom Latham
	* D80: Add CMake options for enabling clang-tidy static analysis checks during build

	14 Apr 2022 John Back
	* D82: EvtDecayProb: initialise data members and remove empty destructor

	14 Apr 2022 Michal Kreps
	* D81: Derive EvtFlatSqDalitz from EvtDecayIncoherent since we directly provide
	final kinematics

	2nd Mar 2022 John Back
	* D78: Add Bc -> J/psi K+ pi- pi+ pi- pi+, Bc -> J/psi K+ K- pi+ pi- pi+ &
	Bc -> J/psi 4pi+ 3pi- decay modes to the BC_VHAD model, courtesy of
	Aleksei Luchinsky, Anna Danilina, Dmitrii Pereima & Vanya Belyaev (LHCb)

	===
	## R02-01-01

	8th Sep 2021 Michal Kreps
	* Update location of web page for Pythia8 download in setup script.

	8th Sep 2021 Markus Prim, Lu Cao, Chaoyi Lyu and Michel De Cian (Michal Kreps)
	* D73: Add new model for semileptonic B decays with BCL and BGL form-factors

	8th June 2021 Michal Kreps
	* T110, D71: Fix B+ --> eta' l nu BF which was order of magnitude too high.
	Balance the decrease by increasing B+ --> D0 l nu, which is after change
	still bit smaller than PDG 2021.

	8th Jun 2021 Michal Kreps
	* D71: Fix B+ --> eta' l nu BF.

	20th Apr 2021 Tom Lathem
	* D68: Fix compilation with Pythia 8.304

	17th Mar 2021 Michal Kreps
	* D62: Improve PI0DALITZ model to dynamically work out maximum probability to
	make it usuable also for eta --> llgamma decays.
	Remove ETA2MUMUGAMMA since it is a pure one-to-one copy of PI0DALITZ.

	15th Jan 2021 Michal Kreps
	* D47: Model to generate 4-body phase-space decays in restricted part of the m12-m34 space

	12th Jan 2021 Michal Kreps
	* D48: Fix bug in calculation of the daughter momentum in decay model EvtBsMuMuKK

	7th Jan 2021 Michal Kreps
	* D43: Allow to pass particle properties table in form of stringstream to constructor
	of EvtGen for use case where these are created on fly.

	10th Dec 2020 Michal Kreps
	* D36: EvtFlatSqDalitz model to be more efficient and to avoid cut-off around the edges

	21st Aug 2020 John Back
	* T109: Add EvtEtaLLPiPi model for eta' -> l l pi pi decays (l = e or mu),
	- courtesy of Aleksei Luchinsky (LHCb).

	29th Jun 2020 Michal Kreps
	* T103: Add missing include in EvtGenBase/EvtMatrix.hh.

	15th May 2020 Michal Kreps
	* D28: Add EvtThreeBodyPhsp (rectangular Dalitz plot selection) and
	EvtPhspDecaytimeCut (minimum decay time requirement) models.
	* D27: Fix initialisation of constants in EvtBTo3hCP model.

	===
	## R02-00-00

	24th Apr 2020 Michal Kreps
	* Update particle properties file evt.pdl to 2019 version of RPP by PDG.

	23rd Apr 2020 Tom Latham
	* Apply copyright and licence notices to all relevant files.

	17th Apr 2020 Tom Latham
	* Add text of GNU GPLv3 in COPYING, add (preliminary) list of authors in
	AUTHORS, and clean-up all source files, prior to applying copyright and
	licence notices.

	9th Apr 2020 Tom Latham
	* Improve, and document use of, standalone installation script.
	* Apply clang-format formatting to all C++ source files.

	8th Apr 2020 Tom Latham
	* One small modernisation change to EvtPhotosEngine to match that already
	applied in EvtTauolaEngine.

	8th Apr 2020 John Back
	* Fixed NonReson amplitude and the 4-momentum boosts used for angles in
	EvtLambdacPHH,
	- courtesy of Elisabeth Niel (LHCb).

	7th Apr 2020 Gerhard Raven, Tom Latham, Michal Kreps and John Back
	* Incorporate C++ modernization changes from Gerhard Raven (LHCb).
	- Merged modernize branch into master.

	9th Mar 2020 John Back
	* Add EvtAbsExternalGen::getDecayProb() to allow external generators to
	return a probability that can be used in EvtDecayProb (default = 1).

	6th Mar 2020 Andrii Verbytskyi and Tom Latham
	* Implement HepMC3 support: EvtHepMCEvent, external engines & cmake files.

	15th Nov 2019 John Back
	* Added EvtPsi2JpsiPiPi model for psi2S -> J/psi pi+ pi- decays based on hep-ph/1507.07985,
	- courtesy of Aleksei Luchinsky (LHCb).

	21st August 2019 Michal Kreps
	* Added the EvtDToKpienu decay model for D -> K pi e nu decays from BESIII,
	- courtesy of Liaoyuan Dong.

	16th July 2019 John Back
	* Correct imaginary sign for EvtComplex /= (EvtComplex c) operator.

	3rd July 2019 John Back
	* Added the EvtLambdacPHH decay model for Lc -> p K pi decays with K*(890),
	Delta++(1232) and Lambda(1520) resonances, based on the Fermilab E791
	analysis hep-ex/9912003v1,
	- courtesy of Elisabeth Niel and Patrick Robbe (LHCb).
	* Modified EvtResonance2 to accept other barrier factor radii.

	3rd July 2019 Michal Kreps
	* Make sure minimum mass for resonances with non-zero widths is larger than
	1e-4 GeV in EvtRelBreitWignerBarrierFact.

	3rd May 2019 John Back
	* Corrected EvtSLDiBaryonAmp bugs/issues in the BToDiBaryonlnupQCD model:
	- parity, amplitude terms and B momentum reference frame variables.
	* Corrected treament of spinor indices in EvtRareLb2Lll,
	- courtesy of Tom Blake and Michal Kreps (LHCb).
	* Updated the EvtBcVHad model to also handle Bc -> psi Ks K decays,
	- courtesy of Aleksei Luchinsky (LHCb).
	* Add new decay model EvtBsMuMuKK (BS_MUMUKK) for Bs to J/psi (mu+mu-) K+K-,
	- courtesy of Veronika Chobanova, Jeremy Dalseno, Diego Martinez Santos and
	Marcos Romero Lamas (LHCb).
	* Fix infinite loop during initialisation of the EvtBTo3hCP model via
	EvtCBTo3piP00 and EvtCBTo3piMPP,
	- courtesy of Peter Richardson (Durham).

	15th March 2019 Tom Latham
	* Implement cmake build system, replacing the old config method.

	30th Jan 2019 John Back
	* Fix modernization compiler errors and warnings.

	29th Jan 2019 Michal Kreps
	* Allow reading decay files which are missing end-of-line before end-of-file.

	21st December 2018 John Back
	* Imported C++ modernization changes from Gerhard Raven (LHCb).

	7th December 2018 John Back
	* Added the EvtBLLNuL (BLLNUL) model that generates rare B -> ell ell nu ell
	decays, where ell = e or mu,
	- courtesy of Anna Danilina and Nikolai Nikitin (LHCb).
	* Removed the EvtB2MuMuMuNu (BUTOMMMN) model, since its now replaced
	by the more general BLLNuL one.

	5th November 2018 John Back
	* Added the BToDiBaryonlnupQCD model for generating B to p N* l nu decays,
	where N can be any (exited) charged baryon (spin 1/2 or 3/2),
	- courtesy of Mark Smith and Ryan Newcombe (LHCb), with added code optimisations.

	17th October 2018 John Back
	* Added various decay models from LHCb EvtGenExtras package:
	- EvtBcVHad ("BC_VHAD"),
	- Evtbs2llGammaMNT ("BSTOGLLMNT"),
	- Evtbs2llGammaISRFSR ("BSTOGLLISRFSR"),
	- EvtbTosllMS ("BTOSLLMS"),
	- EvtbTosllMSExt ("BTOSLLMSEXT"),
	- EvtLb2Baryonlnu ("Lb2Baryonlnu"),
	- EvtLb2plnuLCSR ("Lb2plnuLCSR"),
	- EvtLb2plnuLQCD ("Lb2plnuLQCD"),
	- EvtFlatSqDalitz ("FLATSQDALITZ"),
	- EvtPhspFlatLifetime ("PHSPFLATLIFETIME").

	5th October 2018 John Back
	* Updated setupEvtGen.sh to work with the new HepForge Phabricator site.

	13th March 2018 John Back
	* Updated EvtPythiaEngine to correctly handle updates of various particle
	properties so that Pythia uses the same information as EvtGen (evt.pdl)
	for the generic and alias PYTHIA decay model.

	12th March 2018 John Back
	* Updated EvtBcXMuNu models (X = Scalar, Vector, Tensor) to generate
	Bc to D0(star) mu nu decays, with associated form factors in EvtBCXFF,
	- courtesy of Aleksei Luchinsky (LHCb).
	* Also generalised the calculation
	of their maximum probabilities by reusing the CalcMaxProb method in
	EvtSemiLeptonicAmp, which now allows for different Q^2 binning
	(default remains at 25 bins).

	===
	## R01-07-00

	13th December 2017 John Back
	* New tag incorporating all changes below.
	* Recommended external packages are
	(as used in the setupEvtGen.sh script):
	- HepMC 2.06.09,
	- pythia 8.230,
	- Photos++ 3.61
	- Tauola++ 1.1.6c.

	12th December 2017 John Back
	* Changed Pythia validation example decay files to use Pythia8 codes.

	6th December 2017 John Back
	* Modified the examples to use DECAY.DEC (see 25th April 2016) instead of
	DECAY_2010.DEC. Changed EvtExternalGenList to assume Pythia8 codes are
	used in decay files by default, which is the case for DECAY.DEC. Also
	updated the setupEvtGen.sh script to work with Pythia 8.2x versions.

	29th November 2017 John Back
	* Modified EvtSVP, EvtVVP and EvtTVP models to handle both radiative and
	two-lepton decays,
	- courtesy of Aleksei Luchinsky (LHCb).

	14th July 2017 John Back
	* Only create external generator objects if they don't already exist in
	EvtExternalGenFactory.
	* Modified configure script to work with Pythia 8.2x

	5th July 2017 Michal Kreps
	* Register the VTOSLL model.

	14th June 2017 John Back
	* Add isNeutralKaon() boolean function and corrected comments in EvtDDalitz.

	8th May 2017 Michal Kreps
	* Fix bug in EvtbTosllVectorAmp to recognise Bs --> K*bar mu mu decay as
	b --> d ll transition.

	8th May 2017 Michal Kreps
	* Significantly simplify way how we decide on decay mode and daughters
	ordering in DDalitz model.
	- With new code by definition all orderings of
	daughters in the decay file will yield same output.

	4th May 2017 John Back
	* Further fixes to DDalitz particle ordering (including charge-conjugates):
	- Mode 5: D0 -> K- K0bar K+ and K+ K- K0bar
	- Mode 12: D0 -> pi0 pi- pi+ and pi+ pi0 pi-
	- Removed unneeded index ordering checks for mode 10 (D+ -> pi- pi+ pi+)
	and mode 11 (Ds+ -> pi- pi+ pi+)

	27th April 2017 John Back
	* Fixed DDalitz particle ordering for mode 7: D+ -> pi+ K- K+ and K+ pi+ K-
	and their charge-conjugates

	7th April 2017 John Back
	* Modified EvtGenExternal/EvtPythiaEngine to ensure that the EvtGen-based
	instances of Pythia8 (for generic and alias decays) use the same
	particle properties as defined by EvtGen,
	- courtesy Patrick Robbe (LHCb).

	5th April 2017 Michal Kreps
	* Fixed indexing in copy constructor of Evt3Rank3C, which would otherwise
	produce an infinite loop;
	- bug report from David Grellscheid.

	3rd November 2016 John Back
	* Modified EvtFlatQ2 model to work for all B -> X lepton lepton modes, as
	well as adding an extra phase space factor to correct for the dip at low
	q^2, which is enabled by using "FLATQ2 1" instead of just "FLATQ2" in the decay file,
	- courtesy of Marcin Chrzaszcz & Thomas Blake (LHCb).


	13th October 2016 John Back
	* Added the TauolaCurrentOption decfile keyword to select the hadronic
	current in Tauola; default is the BaBar-tuned current option (int = 1).
	* EvtParticles can store double attributes using the functions
	setAttributeDouble(name, double) and getAttributeDouble(name), which can
	be useful for storing and retrieving amplitude weights, for example.
	- The analogous EvtParticle integer attribute interface remains unchanged:
	setAttribute(name, int) and getAttribute(name).

	13th September 2016 John Back
	* Modified EvtTauolaEngine to use internal Tauola spin matrices for
	tau pair events by temporarily setting the PDG id of the mother as a
	boson, keeping the same 4-momentum.
	* The BaBar hadronic currents are now used by default.
	* Also added the ability to change some Tauola parameters
	using the "Define" keyword in decay files.
	* Added an example decay file
	illustrating the new features: validation/TauolaFiles/Btautau.dec

	9th September 2016 Michal Kreps
	* Reimplement code in EvtBTo3pi.F, EvtBTo3piMPP.F, EvtBTo3piP00.F and
	EvtBToKpipi.F in C++ in order to remove dependence on Fortran compiler.
	- With this, there is no internal Fortran code in EvtGen.

	===
	## R01-06-00

	1st June 2016 John Back
	* New tag incorporating all changes below.
	* Recommended external packages are
	- HepMC 2.06.09
	- pythia 8.186
	- Photos++ 3.61
	- Tauola++ 1.1.5

	28th April 2016 Michal Kreps
	* For Ds+ --> 2pi+ pi- there was double counting of branching fraction
	resulting in total branching fraction being 1.5 times larger than measured
	one.
	- Fix by revisiting submodes, which now fill total Ds --> 3pi.

	25th April 2016 Michal Kreps
	* Added DECAY.DEC/XML, which contain updated semileptonic charm and beauty
	branching fractions using the 2014 PDG, tuned to minimize disagreements
	between measurements and EvtGen for both inclusive and exclusive decays.
	* Updated the evt.pdl particle properties file to the PDG 2014 edition.
	* Implemented new LQCD form factors for Lb --> L mu mu from arXiv paper
	1602.01399 (EvtRareLbToLllFFlQCD); old LQCD form factors are removed.

	18th March 2016 John Back
	* Fixed incorrect spinor algebra used in S1 -> 1/2 S2, 1/2 -> S3 S4 decays
	in EvtDiracParticle::rotateToHelicityBasis() functions,
	- courtesy of Luis Miguel Garcia Martin and the IFIC Valencia LHCb group.

	19th Feburary 2016 John Back
	* Fixed bug in the definition of the initial spinor term Sinit in
	EvtRareLbToLll::HadronicAmpRS(),
	- from Tom Blake (LHCb).

	12th February 2016 John Back
	* From LHCb, added extensions to the EvtHQET2(FF) model for semitauonic
	decays from Brian Hamilton, which needs a patch to EvtSemiLeptonicAmp
	from Jack Wimberley to ensure that the q^2 range is physical when
	finding the maximum amplitude probability.

	2nd December 2015 John Back
	* From LHCb, added EvtKStopizmumu model for KS -> pi0 mu mu decays based on
	JHEP08(1998)004,
	- courtesy of Veronika Chobanova, Diego Martinez Santos and Jeremy Dalseno.
	* Added EvtConst::Fermi for Fermi coupling constant.

	===
	## R01-05-00

	21st October 2015 John Back
	* New tag incorporating all changes below.
	* Recommended external packages are
	- HepMC 2.06.09
	- pythia 8.186
	- Photos++ 3.61
	- Tauola++ 1.1.5

	* Added the EvtB2MuMuMuNu model for simulating the very rare four-leptonic
	decays B- -> mu+ mu- anti-nu_mu mu-,
	- courtesy Nikolai Nikitin.

	16th October 2015 John Back
	* Updated the configure script to automatically select the library names
	for PHOTOS++; version 3.56 and below uses Fortran, version 3.61 and above
	uses C++ only (default). Avoid using v3.60, since it does not work.
	This needs the PHOTOS libraries built before EvtGen is configured.
	Modified setupEvtGen.sh to use Photos++ v3.61.

	7th October 2015 John Back
	* Updated EvtGenExternal/EvtPhotosEngine to check that additional particles
	from the outgoing vertex are indeed (FSR) photons, since later versions of
	PHOTOS introduce pair emission, where particles may not always be photons.
	* Added the genRootDecayChain.cc validation program to create ROOT files
	containing information about the complete decay tree. Two example test
	decay files BKstarGamma.dec and BuDst0rhop.dec can be used with this; the
	first tests PHOTOS, the second looks at sequential decay chain storage.
	The plotBKstarGamma.C ROOT macro can be used for B -> K* gamma plots.

	2nd October 2015 John Back
	* Modified EvtSVPHelAmp and added a new EvtSVPHelCPMix model, implementing
	the complete mixing phenomenology of Bs to vector gamma decays,
	- courtesy of Clara Remon (LHCb).
	* EvtD0mixDalitz code: cleanup, inverted q/p for decays of D0bar (simplifies
	user decay files) and fixed y parameter bug,
	- courtesy of Jordi Tico (LHCb).
	* Changed the initialisation order of the infrared cut-off in EvtPhotosEngine.
	This actually has no effect, since the exponentiation function sets it to the
	same 1e-7 value, but it is now in the correct order if we need to update it.
	* Removed all remaining obsolete pragma (Win32) warnings from some classes.

	23rd September 2015 Michal Kreps
	* Reimplement the real Spence function in C++ and removed its fortran
	implementation.

	15th September 2015 Michal Kreps
	* Fixed accessed uninitialised memory in EvtPDL.cpp, line 213.
	* Modified the configure and setupEvtGen.sh scripts to work on Mac; needed
	Mac compilation patch files added to the new "platform" subdirectory.

	10th September 2015 John Back
	* Updated setupEvtGen.sh to use the recommended external packages:
	- HepMC 2.06.09, pythia 8.186, Photos++ 3.56 and Tauola++ 1.1.5.
	* Fixed form-factor calculations for the BTOSLLBALL model 6 used to generate
	b -> sll decays,
	- courtesy of Christoph Langenbruch and David Loh (LHCb).
	- Affects B->K*ll, B->rholl and B->omegall, particularly the electron modes.
	* In the validation directory, added runPhotosTest.sh for testing FSR in
	Upsilon(4S) -> e+ e- decays, and changed the plot comparison scripts to
	use the 2nd directory "oldRootFiles" (which could be a soft-link) for
	including ROOT histograms made from a previous version of EvtGen.

	27th August 2015 John Back
	* Added Mersenne-Twister random number generator (RNG) EvtMTRandomEngine.
	- It requires c++11 compiler features (>= gcc 4.7), which should
	automatically be enabled by the configure script.
	- Introduced the preprocessor environment variable EVTGEN_CPP11 for c++11
	features.
	- EvtMTRandomEngine is the default RNG for the validation and test examples
	if c++11 features are enabled.
	* Added a phase-space test validation/genPHSP.sh and PhaseSpacePlots.C to
	visually check the flatness of Dalitz plots in order to ensure that the
	RNG is not producing biased results that depend on particle ordering.
	* Added the models EvtbsToLLLLAmp and EvtbsToLLLLHyperCP for
	B0_q -> l+ l- l+ l- decays (SM and one supersymmetric scenario),
	- courtesy of Nikolai Nikitin and Konstantin Toms.
	- Documentation provided in doc/evt_BQTOLLLL_model.pdf and
	doc/evt_BQTOLLLLHYPERCP_model.pdf.

	* Changed the installation and set-up script name to be just setupEvtGen.sh;
	it uses the VERSION variable to specify the required tag. List of tags
	are available using either "svn ls -v http://svn.cern.ch/guest/evtgen/tags"
	or by going to http://svn.cern.ch/guest/evtgen/tags in a web browser.

	12th June 2015 John Back
	* Changed the width of chi_b1 in evt.pdl from 9.8928 GeV (!) to zero.

	1st May 2015 John Back
	* Added Bc -> scalar ell nu (EvtBcSMuNu) and Bc -> tensor ell nu
	(EvtBcTMuNu) decays,
	- courtesy of Jack Wimberley, LHCb.
	- Also included the chi_c1 mode for EvtBcVMuNu.

	===
	## R01-04-00

	2nd April 2015 John Back
	* Removed the EvtStdlibRandomEngine class since this can produce biases
	to kinematic distributions when one or more of the daughters is a
	resonance, such as B0 -> K pi psi
	- (thanks to Antonio Augusto Alves Jr who discovered this issue).
	- EvtSimpleRandomEngine is now the default random number generator in the
	validation and test examples.
	* Incorporated several additions and modifications from LHCb:
	a) From Michal Kreps, Tom Blake & Christoph Langenbruch,
	added rare Lb --> Lambda^(*) ell ell models described in
	arXiv:1108.6129, with various form factors from Gutsche et al.
	(arXiv:1301.3737) and lattice QCD (arXiv:1212.4827)
	b) From Andrew Crocombe, implemented Bs --> K* form factors
	from Ball-Zwicky and z-parametrization form factors from
	arXiv:1006.4945 for EvtbTosllBallFF
	c) Christoph Langenbruch fixed the Bs -> phi ll form factors in
	EvtbTosllBallFF; T3 showed a non-physical pole at very low
	q2 which significantly affected the electron mode
	d) From Michal Kreps, removed semicolons from wrong places to
	clear warnings when compiled with the "-pedantic" option.

	9th October 2014 John Back
	* Change svnweb.cern.ch to svn.cern.ch in the setup script.

	1st April 2014 John Back
	* In EvtReport, modified the logging output severity status flags
	to have the "EVTGEN_" prefix, e.g. INFO becomes EVTGEN_INFO.
	* The global report() function has been renamed to EvtGenReport().

	31st March 2014 John Back
	* Added the ability to store named attributes for EvtParticles in the
	form of a map<string, int>. The setAttribute(name, value) stores the
	required value, while getAttribute(name) retrieves the integer value.
	This is used in EvtPhotosEngine to specify the final-state radiation
	"FSR" attribute to 1 for any additional photons (EvtPhotonParticles)
	created by Photos++. It also stores the "ISR" attribute, but this
	is always set to zero, since only FSR photons are created.
	If the named attribute does not exist, then getAttribute() returns zero.

	29th January 2014 Daniel Craik
	* Removed mass assertion on GS shape in EvtDalitzReso to allow it to also
	be used for charged rho resonances.

	27th January 2014 John Back
	* Minor corrections to Vub models to remove further gcc 4.8 warnings.
	* Updated configure script to work for MacOS clang (from Genser team).

	===
	## R01-03-00

	9th January 2014 John Back
	* New tag version "1.3.0", incorporating all changes below.
	* Replaced auto-install script to work with this version as well as
	the latest versions of all external generator packages.
	* Updated README to mention the new CERN-based web pages for Photos++
	and Tauola++.

	8th January 2014 John Back
	* Fix gcc 4.6 and 4.8 compilation warnings,
	- courtesy of Patrick Robbe (LHCb);
	- main changes are removal of unused variables.
	* Changed the EvtPythiaEngine class and configure script to use new
	Pythia 8 header locations; Pythia 8.180 or above is now required.

	7th January 2014 John Back
	* Modified EvtBCVFF to correct the Kiselev form factors
	- from Jack Wimberley (LHCb).

	9th October 2013 Daniel Craik
	* Added Gounaris-Sakurai and Gaussian shapes to EvtGenericDalitz
	and set sensible defaults for LASS parameters.

	19th September 2013 John Back
	* Modified EvtGenExternal/EvtPythiaEngine to keep track of any new
	particles that are added to the default Pythia database to avoid
	duplicating particle/anti-particle entries that could override
	previously defined Pythia decay chains.

	18th September 2013 John Back
	* Added Mac OS flags for the configure script and src/Makefile.

	15th July 2013 Daniel Craik
	* Added flag to turn on scaling of LASS amplitude by M/q in EvtDalitzReso

	15th July 2013 Daniel Craik
	* EvtParserXML now accepts file names containing environment variables,
	exponential non-resonant shape in EvtDalitzReso now defined as exp(-alpha*m^2),
	LASS shape in EvtDalitzReso now takes a cutoff parameter

	4th July 2013 Daniel Craik
	* Added LASS, exponential non-resonant and linear non-resonant shapes to EvtGenericDalitz.

	3rd July 2013 Daniel Craik
	* Fixed auto-install script for R01-02-00.

	1st July 2013 Daniel Craik
	* Added auto-install script for R01-02-00.

	===
	## R01-02-00

	15th May 2013 John Back
	* New tag, version "1.2.0", incorporating all changes below.

	14th May 2013 Michal Kreps
	* Added Blatt-Weisskopf barrier factors up to L=5 in
	EvtGenBase/EvtBlattWeisskopf::compute().

	14th May 2013 John Back
	* Added additional entries (appended at the end) to the evt.pdl particle
	data file
	- courtesy of Romulus Godang and Belle II colleagues.

	14th March 2013 John Back
	* Added the method EvtParticle::getPDGId() to get the PDG integer for a
	particle directly (which just calls EvtPDL::getStdHep()).
	* Added a check in EvtPhotosEngine::doDecay to skip Photos if a given
	particle has too many daughters (>= 10) to avoid a problem with a
	hard coded upper limit in the Photos PHOENE subroutine.

	2nd February 2013 Daniel Craik
	* Updated EvtDalitzTable to estimate probMax when it is missing from a
	Dalitz model.

	1st February 2013 John Back
	* Added the ability to read in Pythia 6 commands in ascii decay files in
	EvtDecayTable::readDecayFile (this was already possible in xml files).
	* Modified the Photos++ engine default settings to be more suited to B
	decays (from LHCb defaults).

	31st January 2013 John Back
	* Added the ability to read in Pythia 8 commands in ascii decay files
	in EvtDecayTable::readDecayFile. They can be set using the syntax:
	"PythiaTypeParam module:variable=value", where Type = Generic, Alias or
	Both for specifying whether the parameter is for the generic or alias
	Pythia decay engines (or both). The 2nd argument must not contain any
	spaces. Fixed the list of commands strings being used in the
	EvtPythiaEngine class (i.e. Pythia parameters that can be set via decay
	files).

	31st January 2013 Daniel Craik
	* Added named parameters to various decay models.

	30th January 2013 John Back
	* Fixed some of the parameter arguments used in the EvtSVSCPiso model.

	24th January 2013 John Back
	* Set the Photos++ and Tauola++ models to use the EvtGen random number
	engine when useEvtGenRandom is set to true in the EvtExternalGenList
	constructor.

	23rd January 2013 John Back
	* Added EvtGenExternal/EvtPythiaRandom to allow the use of the EvtGen
	random number engine to also be used for the random engine for Pythia 8.
	* Added a boolean (useEvtGenRandom, default = true) within the
	EvtExternalGenList constructor to use this feature.

	18th December 2012 John Back
	* Corrected some wrong daughter ordering assignments for decay modes 5 and
	12 in EvtDDalitz. Updated validation/DalitzDecays.xml to also contain
	D decay mode 12, as well as various modes that may use K_S0 and K_L0.
	* Added validation/genDDalitzModes.sh and updated validation/compareDalitz.C to
	do a complete comparison of the D Dalitz modes with the xml versions.

	11th December 2012 Daniel Craik
	* Updated the Xml parser to support named model parameters.
	* Updated the generic Dalitz model to use named model parameters as an example.

	15th October 2012 John Back
	* Make EvtSimpleRandomEngine inherit from EvtRandomEngine to avoid
	crash in EvtGen.cpp when no random engine is defined
	- (from Bjoern Spruck).

	===
	## R01-01-00

	4th October 2012 John Back
	* New tag, version "1.1.0", incorporating all changes below.

	* Provide proper default constructors for EvtVector4R and
	EvtPhotonParticle. Modified the validation and test code to also
	compile/link in the case of no external generators being included.

	3rd October 2012 John Back
	* Corrected the t3 vector form factor values for the Ball-Zwicky 2005
	model (modelId = 6) in EvtbTosllBallFF::getVectorFF(), which
	were set to t3tilde instead.

	18th September 2012 John Back
	* Moved the external generator engines to a new sub-directory
	EvtGenExternal. Building the code now creates 2 libraries:
	libEvtGen.so (Base+Models) and libEvtGenExternal.so.
	- This allows anyone to ignore using the new external generators
	if required (by not creating/loading the 2nd library).
	* Added prefix option to the configure script/Makefile to allow the user
	to specify an installation directory for the include files, libraries,
	DECAY.DEC and evt.pdl files (for Genser).

	14th September 2012 Michal Kreps
	* Fixed the calculation of the angle between decay planes in the function
	EvtKine::EvtDecayAngleChi. Fixed typo in EvtLb2Lll decay model. Only
	some NP scenarious could be affected, SM one is definitely unaffected.

	13th September 2012 John Back
	* Added the use of the environment variables EVTGEN_PHOTOS, EVTGEN_PYTHIA
	and EVTGEN_TAUOLA to specify if the Photos, Pythia and/or Tauola engine
	classes are used or not. These variables are set by the configure script,
	depending if the library paths are specified for these generators.

	===
	## R01-00-01

	12th September 2012 John Back
	* New tag incorporating all changes below, since R01-00-00.

	11th September 2012 John Back
	* Modified the Photos and Tauola engine classes to use the new
	Photospp and Tauolapp namespaces that are present in the latest
	versions of Photos++(3.5) and Tauola++(1.0.7).
	* Updated the configure file to get the correct location of the Tauola++
	include files.
	* Added the D0->pi+pi-pi0 decay mode in EvtDDalitz
	- from Marco Gersabeck and Frederic Dreyer (LHCb).
	* Added new decay models/classes from Alexey Luchinsky (LHCb):
	EvtBcVMuNu, EvtTVP, EvtWnPi, EvtSVP, EvtXPsiGamma, EvtBcVNpi

	29th June 2012 John Back
	* Corrected mass(squared) variables filled in the Dalitz TTree in
	validation/genExampleRootFiles.

	15th May 2012 Daniel Craik
	* Updated EvtD0gammaDalitz to deal with D mesons from neutral B->DK
	* Added save function to validation/compareDalitz.C.

	11th May 2012 Daniel Craik
	* Replaced BaBar specific configuration for BlattWeisskopf birth factors.
	* Updated XML conversion script to handle new configuration.
	* Fixed some bugs in the XML conversion script related to particle
	modifications.

	9th May 2012 Daniel Craik
	* Added latex documentation for xml decay files.

	2nd May 2012 Daniel Craik
	* Added resDaughters attribute to the Dalitz resonance xml tag to
	simplify defining symmetric resonances. Updated validation xml files to
	use the new functionality.

	27th April 2012 Daniel Craik
	* Upgraded EvtGenericDalitz to use EvtDalitzReso for resonances.
	* Added validation to compare EvtGenericDalitz to all 11 EvtDDalitz modes.
	* Added a root macro to quickly compare two Dalitz decays for validation.

	24th April 2012 John Back
	* Solved two bugs in the EvtD0gammaDalitz model (from Jordi Tico, LHCb):
	configuration of the conjugated model, and using only the B charge
	to determine the model used, not the D flavour.

	17th April 2012 Daniel Craik
	* Updated the GenericDalitz validation code to use the same probMax
	values as DDalitz.
	* Added XML decay file parsing to EvtGen::readUDecay.
	- Dec files are still the default.

	30th March 2012 John Back
	* Update maximum probability values in EvtDDalitz::initProbMax()
	for all DDalitz modes.

	23rd March 2012 John Back
	* Added the EvtEta2MuMuGamma decay model from LHCb.

	21st March 2012 John Back
	* Added EvtD0gammaDalitz decay model from LHCb.

	20th March 2012 Daniel Craik
	* Added backwards compatibility for Pythia 6 commands in the XML configuration.
	* Updated decay file conversion tool to convert JetSetPar lines to pythia6Param
	tags.

	19th March 2012 Daniel Craik
	* Added infrastructure to pass commands to external generators.
	* XML config now takes Pythia8 configuration commands.

	16th March 2012 Daniel Craik
	* Added the ability to define particles from the PDL for Dalitz decay
	resonances instead of defining mass, width and spin seperately.
	* Renamed the lifetime attribute of Dalitz decay resonaces to width to avoid
	confusion.
	* Added further validation code for the generic Dalitz model.

	15th March 2012 Daniel Craik
	* Added validation code for xml decay files and the generic Dalitz model.

	===
	## R01-00-00

	6th March 2012 John Back
	* First official version for Genser (evtgen 1.0.0) that includes
	support for external generators: Pythia8, Photos++ and Tauola++.
	* This also includes a preliminary version of creating Dalitz plot
	decay models using EvtGenericDalitz.
	diff --git a/cmake/Modules/EvtGenExternalDependencies.cmake b/cmake/Modules/EvtGenExternalDependencies.cmake
	index 3840341..9851f41 100644
	--- a/cmake/Modules/EvtGenExternalDependencies.cmake
	+++ b/cmake/Modules/EvtGenExternalDependencies.cmake
	@@ -1,67 +1,66 @@

	########################################################################
	# 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/>. #
	########################################################################


	set(HEPMC2_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of HepMC 2 installation")
	set(HEPMC3_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of HepMC 3 installation")

	set(PYTHIA8_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of Pythia8 installation")
	set(Photos++_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of Photos++ installation")
	set(Tauola++_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of Tauola++ installation")
	set(PHOTOSPP_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of Photos++ installation - alternative spelling")
	set(TAUOLAPP_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of Tauola++ installation - alternative spelling")
	set(SHERPA_ROOT_DIR "${CMAKE_INSTALL_PREFIX}" CACHE PATH "Location of Sherpa installation")

	# The components we search for in the external generators depend on the version
	# of HepMC we're working with
	if(EVTGEN_HEPMC3)
	find_package(HepMC3 REQUIRED PATHS ${HEPMC3_ROOT_DIR})
	- if(EVTGEN_PYTHIA)
	- find_package(Pythia8 REQUIRED)
	- endif()
	if(EVTGEN_PHOTOS)
	# From version 3.64 Photos has HepMC3 support
	find_package(Photos++ REQUIRED COMPONENTS pp ppHepMC3)
	endif()
	if(EVTGEN_TAUOLA)
	# From version 1.1.8 Tauola has HepMC3 support
	find_package(Tauola++ REQUIRED COMPONENTS Fortran CxxInterface HepMC3)
	endif()
	- if(EVTGEN_SHERPA)
	- find_package(Sherpa 2.0.0 REQUIRED)
	- endif()
	else()
	find_package(HepMC2 REQUIRED)
	- if(EVTGEN_PYTHIA)
	- find_package(Pythia8 REQUIRED)
	- endif()
	if(EVTGEN_PHOTOS)
	# Photos has different library structures for versions before and after 3.58
	# so we need to search for either option: pp+ppHepMC or CxxInterface+Fortran
	find_package(Photos++ REQUIRED OPTIONAL_COMPONENTS pp ppHepMC CxxInterface Fortran)
	endif()
	if(EVTGEN_TAUOLA)
	# Older versions of Tauola don't have the HepMC component, the HepMC2 interface is in CxxInterface
	find_package(Tauola++ REQUIRED COMPONENTS Fortran CxxInterface OPTIONAL_COMPONENTS HepMC)
	endif()
	- if(EVTGEN_SHERPA)
	- find_package(Sherpa 2.0.0 REQUIRED)
	+endif()
	+
	+if(EVTGEN_PYTHIA)
	+ find_package(Pythia8 REQUIRED)
	+endif()
	+if(EVTGEN_SHERPA)
	+ if(EVTGEN_SHERPA3)
	+ find_package(SHERPA-MC 3.0.0 CONFIG REQUIRED)
	+ else()
	+ find_package(Sherpa 2.0.0 MODULE REQUIRED)
	endif()
	endif()
	diff --git a/setupEvtGen.sh b/setupEvtGen.sh
	index b19da7d..b8d0319 100644
	--- a/setupEvtGen.sh
	+++ b/setupEvtGen.sh
	@@ -1,359 +1,365 @@
	#!/bin/bash

	########################################################################
	# Copyright 1998-2024 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/>. #
	########################################################################

	# This script installs EvtGen with all external dependencies.
	# The variable INSTALL_PREFIX specifies the installation location.
	# The variable VERSION specifies the tag of EvtGen you want to use.
	# The list of available tags can be found by either going to the url
	# https://phab.hepforge.org/source/evtgen/tags/master
	# or issuing the command (without the need to clone the git repository)
	# git ls-remote --tags http://phab.hepforge.org/source/evtgen.git \| cut -d '/' -f3
	# The recommended versions of the external dependencies are given below.
	# Later versions should be OK as well, assuming their C++ interfaces do not change.
	# HepMC (either HepMC2 or HepMC3, the latter is recommended) is mandatory.
	# Whether to use a particular external dependency can be specified by setting
	# the corresponding USEXXX variable to "ON" or "OFF".
	# Note that some earlier EvtGen versions will not be compatible with all
	# external dependency versions given below, owing to C++ interface differences;
	# see the specific tagged version of the EvtGen/README file for guidance.
	# It is also not possible to compile Tauola++ on macOS at present, unless you
	# are building on a volume with a case sensitive file system, so this is
	# disabled by default.
	# Similarly, we have encountered problems building Sherpa on macOS, so that it
	# also disabled by default. NB Sherpa is anyway not used by EvtGen versions < 3.
	# To obtain this script use the "Download File" option on the right of the webpage:
	# https://phab.hepforge.org/source/evtgen/browse/master/setupEvtGen.sh?view=raw

	# Location in which to install
	INSTALL_PREFIX="/usr/local"

	# EvtGen version or tag number (or branch name)
	VERSION="R03-00-00-beta1"

	# EvtGen git URL
	# HepForge main repo
	EVTGEN_URL="https://phab.hepforge.org/source/evtgen.git"
	# CERN GitLab mirror
	#EVTGEN_URL="https://gitlab.cern.ch/evtgen/evtgen.git"

	# Options to pass to 'make', e.g. "-j 6"
	BUILDARGS=""

	# HepMC version numbers - change HEPMCMAJORVERSION to 2 in order to use HepMC2
	HEPMCMAJORVERSION="3"
	HEPMC2VER="2.06.11"
	HEPMC3VER="3.3.0"
	HEPMC2PKG="HepMC-$HEPMC2VER"
	HEPMC3PKG="HepMC3-$HEPMC3VER"
	HEPMC2TAR="hepmc$HEPMC2VER.tgz"
	HEPMC3TAR="$HEPMC3PKG.tar.gz"
	HEPMCBASEURL="https://hepmc.web.cern.ch/hepmc/releases"
	if [ "$HEPMCMAJORVERSION" -lt "3" ]
	then
	+ USEHEPMC3="OFF"
	HEPMCURL=$HEPMCBASEURL/$HEPMC2TAR
	+ HEPMCTAR=$HEPMC2TAR
	else
	+ USEHEPMC3="ON"
	HEPMCURL=$HEPMCBASEURL/$HEPMC3TAR
	+ HEPMCTAR=$HEPMC3TAR
	fi

	# Pythia version number with no decimal points, e.g. 8310 corresponds to version 8.310
	# This follows the naming convention of Pythia install tar files
	USEPYTHIA="ON"
	PYTHIAVER="8312"
	PYTHIAPKG="pythia$PYTHIAVER"
	PYTHIATAR="$PYTHIAPKG.tgz"
	PYTHIABASEURL="https://pythia.org/download"
	PYTHIAURL="$PYTHIABASEURL/pythia${PYTHIAVER:0:2}/$PYTHIATAR"

	# Photos++ version number
	USEPHOTOS="ON"
	PHOTOSVER="3.64"
	PHOTOSPKG="PHOTOS"
	PHOTOSDIR="PHOTOS.$PHOTOSVER"
	PHOTOSTAR="$PHOTOSDIR.tar.gz"
	PHOTOSBASEURL="https://photospp.web.cern.ch/resources"
	PHOTOSURL="$PHOTOSBASEURL/$PHOTOSDIR/$PHOTOSTAR"

	# Tauola++ version number
	USETAUOLA="ON"
	TAUOLAVER="1.1.8"
	TAUOLAPKG="TAUOLA"
	TAUOLADIR="TAUOLA.$TAUOLAVER"
	TAUOLATAR="$TAUOLADIR.tar.gz"
	TAUOLABASEURL="https://tauolapp.web.cern.ch/resources"
	TAUOLAURL="$TAUOLABASEURL/$TAUOLADIR/$TAUOLATAR"

	# Sherpa version number
	USESHERPA="ON"
	-SHERPAVER="2.2.16"
	+SHERPAMAJORVERSION="3"
	+if [ "$SHERPAMAJORVERSION" -lt "3" ]
	+then
	+ USESHERPA3="OFF"
	+ SHERPAVER="2.2.16"
	+else
	+ USESHERPA3="ON"
	+ SHERPAVER="3.0.0"
	+fi
	SHERPAPKG="sherpa-v$SHERPAVER"
	SHERPATAR="$SHERPAPKG.tar.gz"
	SHERPABASEURL="https://gitlab.com/sherpa-team/sherpa/-/archive"
	SHERPAURL="$SHERPABASEURL/v$SHERPAVER/$SHERPATAR"

	# Determine OS
	osArch=`uname`

	# macOS settings
	if [ "$osArch" == "Darwin" ]
	then
	export LANG=en_US.UTF-8
	export LC_ALL=en_US.UTF-8
	# need to disable Tauola
	USETAUOLA="OFF"
	USESHERPA="OFF"
	fi

	#This is for systems with cmake and cmake3
	if command -v cmake3; then
	CMAKE=cmake3
	else
	CMAKE=cmake
	fi

	echo Will install EvtGen version $VERSION and its dependencies in $INSTALL_PREFIX

	# Make sure that the dependencies are picked up from our local install
	delpath()
	{
	eval "$1=\$(echo \$$1 \| sed -e s%^$2\$%% -e s%^$2\:%% -e s%:$2\:%:%g -e s%:$2\\\$%%)"
	}
	addpath_front()
	{
	delpath $*
	eval "$1=$2:\$$1"
	}
	addpath_front CMAKE_PREFIX_PATH $INSTALL_PREFIX
	export CMAKE_PREFIX_PATH

	BUILD_BASE=`mktemp -d` \|\| exit 1

	echo Temporary build area is $BUILD_BASE

	cd $BUILD_BASE

	mkdir -p tarfiles
	mkdir -p sources
	mkdir -p builds

	echo Downloading EvtGen source from GIT, using URL: $EVTGEN_URL

	cd sources
	git clone $EVTGEN_URL evtgen
	cd evtgen
	git checkout $VERSION

	echo Downloading sources of external dependencies

	cd $BUILD_BASE/tarfiles

	echo Downloading HepMC source from: $HEPMCURL
	curl -O $HEPMCURL
	if [ $USEPYTHIA == "ON" ]
	then
	echo Downloading Pythia8 source from: $PYTHIAURL
	curl -O $PYTHIAURL
	fi
	if [ $USEPHOTOS == "ON" ]
	then
	echo Downloading Photos source from: $PHOTOSURL
	curl -O $PHOTOSURL
	fi
	if [ $USETAUOLA == "ON" ]
	then
	echo Downloading Tauola source from: $TAUOLAURL
	curl -O $TAUOLAURL
	fi
	if [ $USESHERPA == "ON" ]
	then
	echo Downloading Sherpa source from: $SHERPAURL
	curl -O $SHERPAURL
	fi

	cd $BUILD_BASE/sources

	echo Extracting external dependencies
	-if [ "$HEPMCMAJORVERSION" -lt "3" ]
	-then
	- tar -xzf $BUILD_BASE/tarfiles/$HEPMC2TAR
	-else
	- tar -xzf $BUILD_BASE/tarfiles/$HEPMC3TAR
	-fi
	+tar -xzf $BUILD_BASE/tarfiles/$HEPMCTAR
	if [ $USEPYTHIA == "ON" ]
	then
	tar -xzf $BUILD_BASE/tarfiles/$PYTHIATAR
	fi
	if [ $USEPHOTOS == "ON" ]
	then
	tar -xzf $BUILD_BASE/tarfiles/$PHOTOSTAR
	fi
	if [ $USETAUOLA == "ON" ]
	then
	tar -xzf $BUILD_BASE/tarfiles/$TAUOLATAR
	fi
	if [ $USESHERPA == "ON" ]
	then
	tar -xzf $BUILD_BASE/tarfiles/$SHERPATAR
	fi

	if [ "$osArch" == "Darwin" ]
	then
	# Patch PHOTOS and TAUOLA on Darwin (macOS)
	if [ $USEPHOTOS == "ON" ]
	then
	sed -i '' 's/soname/install_name/g' PHOTOS/Makefile
	patch -p0 < $BUILD_BASE/sources/evtgen/platform/photos_Darwin.patch
	fi
	if [ $USETAUOLA == "ON" ]
	then
	sed -i '' 's/soname/install_name/g' TAUOLA/Makefile
	patch -p0 < $BUILD_BASE/sources/evtgen/platform/tauola_Darwin.patch
	fi

	# Uncomment the lines below to force usage of Apple clang
	# export CC=clang
	# export CXX=clang++
	# sed -i '' 's/\-lstdc++/-lc++/g' PHOTOS/platform/make.inc.in
	# sed -i '' 's/\-lstdc++/-lc++/g' TAUOLA/platform/make.inc.in
	fi

	cd $BUILD_BASE

	if [ "$HEPMCMAJORVERSION" -lt "3" ]
	then
	# HepMC2 install

	echo Installing HepMC from $BUILD_BASE/sources/$HEPMC2PKG
	mkdir -p $BUILD_BASE/builds/HepMC2
	cd $BUILD_BASE/builds/HepMC2
	$CMAKE -DCMAKE_INSTALL_PREFIX:PATH=$INSTALL_PREFIX $BUILD_BASE/sources/$HEPMC2PKG -Dmomentum:STRING=GEV -Dlength:STRING=MM
	make $BUILDARGS
	make install

	if [ $USEPYTHIA == "ON" ]
	then
	echo Installing Pythia8 from $BUILD_BASE/sources/$PYTHIAPKG
	cd $BUILD_BASE/sources/$PYTHIAPKG
	./configure --enable-shared --prefix=$INSTALL_PREFIX
	make $BUILDARGS
	make install
	fi

	if [ $USEPHOTOS == "ON" ]
	then
	echo Installing PHOTOS from $BUILD_BASE/sources/$PHOTOSPKG
	cd $BUILD_BASE/sources/$PHOTOSPKG
	./configure --with-hepmc3= --with-hepmc=$INSTALL_PREFIX --prefix=$INSTALL_PREFIX
	make $BUILDARGS
	make install
	fi

	if [ $USETAUOLA == "ON" ]
	then
	echo Installing TAUOLA from $BUILD_BASE/sources/$TAUOLAPKG
	cd $BUILD_BASE/sources/$TAUOLAPKG
	./configure --without-hepmc3 --with-hepmc=$INSTALL_PREFIX --prefix=$INSTALL_PREFIX
	make $BUILDARGS
	make install
	fi

	else
	# HepMC3 install

	echo Installing HepMC3 from $BUILD_BASE/sources/$HEPMC3PKG
	mkdir -p $BUILD_BASE/builds/HepMC3
	cd $BUILD_BASE/builds/HepMC3
	$CMAKE -DHEPMC3_ENABLE_ROOTIO:BOOL=OFF -DHEPMC3_ENABLE_PYTHON:BOOL=OFF -DCMAKE_INSTALL_PREFIX:PATH=$INSTALL_PREFIX $BUILD_BASE/sources/$HEPMC3PKG
	make $BUILDARGS
	make install

	if [ $USEPYTHIA == "ON" ]
	then
	echo Installing Pythia8 from $BUILD_BASE/souces/$PYTHIAPKG
	cd $BUILD_BASE/sources/$PYTHIAPKG
	./configure --enable-shared --prefix=$INSTALL_PREFIX
	make $BUILDARGS
	make install
	fi

	if [ $USEPHOTOS == "ON" ]
	then
	echo Installing PHOTOS from $BUILD_BASE/sources/$PHOTOSPKG
	cd $BUILD_BASE/sources/$PHOTOSPKG
	./configure --without-hepmc --with-hepmc3=$INSTALL_PREFIX --prefix=$INSTALL_PREFIX
	make $BUILDARGS
	make install
	fi

	if [ $USETAUOLA == "ON" ]
	then
	echo Installing TAUOLA from $BUILD_BASE/sources/$TAUOLAPKG
	cd $BUILD_BASE/sources/$TAUOLAPKG
	./configure --without-hepmc --with-hepmc3=$INSTALL_PREFIX --prefix=$INSTALL_PREFIX
	make $BUILDARGS
	make install
	fi
	fi

	if [ $USESHERPA == "ON" ]
	then
	- echo Installing Sherpa from $BUILD_BASE/sources/$SHERPAPKG
	- cd $BUILD_BASE/sources/$SHERPAPKG
	- autoreconf -i
	- ./configure --with-sqlite3=install --prefix=$INSTALL_PREFIX
	- make $BUILDARGS
	- make install
	+ echo Installing Sherpa $SHERPAMAJORVERSION from $BUILD_BASE/sources/$SHERPAPKG
	+ if [ "$SHERPAMAJORVERSION" -lt "3" ]
	+ then
	+ cd $BUILD_BASE/sources/$SHERPAPKG
	+ autoreconf -i
	+ ./configure --with-sqlite3=install --prefix=$INSTALL_PREFIX
	+ make $BUILDARGS
	+ make install
	+ else
	+ mkdir -p $BUILD_BASE/builds/Sherpa
	+ cd $BUILD_BASE/builds/Sherpa
	+ $CMAKE -DCMAKE_INSTALL_PREFIX:PATH=$INSTALL_PREFIX $BUILD_BASE/sources/$SHERPAPKG
	+ $CMAKE --build .
	+ $CMAKE --install .
	+ fi
	fi

	echo Installing EvtGen from $BUILD_BASE/sources/evtgen
	mkdir -p $BUILD_BASE/builds/evtgen
	cd $BUILD_BASE/builds/evtgen
	-if [ "$HEPMCMAJORVERSION" -lt "3" ]
	-then
	- $CMAKE -DCMAKE_INSTALL_PREFIX:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_HEPMC3:BOOL=OFF -DHEPMC2_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_PYTHIA:BOOL=$USEPYTHIA -DPYTHIA8_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_PHOTOS:BOOL=$USEPHOTOS -DPHOTOSPP_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_TAUOLA:BOOL=$USETAUOLA -DTAUOLAPP_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_SHERPA:BOOL=$USESHERPA -DSHERPA_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- $BUILD_BASE/sources/evtgen
	-else
	- $CMAKE -DCMAKE_INSTALL_PREFIX:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_HEPMC3:BOOL=ON -DHEPMC3_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_PYTHIA:BOOL=$USEPYTHIA -DPYTHIA8_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_PHOTOS:BOOL=$USEPHOTOS -DPHOTOSPP_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_TAUOLA:BOOL=$USETAUOLA -DTAUOLAPP_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- -DEVTGEN_SHERPA:BOOL=$USESHERPA -DSHERPA_ROOT_DIR:PATH=$INSTALL_PREFIX \
	- $BUILD_BASE/sources/evtgen
	-fi
	+$CMAKE -DCMAKE_INSTALL_PREFIX:PATH=$INSTALL_PREFIX \
	+ -DEVTGEN_HEPMC3:BOOL=$USEHEPMC3 -DHEPMC${HEPMCMAJORVERSION}_ROOT_DIR:PATH=$INSTALL_PREFIX \
	+ -DEVTGEN_PYTHIA:BOOL=$USEPYTHIA -DPYTHIA8_ROOT_DIR:PATH=$INSTALL_PREFIX \
	+ -DEVTGEN_PHOTOS:BOOL=$USEPHOTOS -DPHOTOSPP_ROOT_DIR:PATH=$INSTALL_PREFIX \
	+ -DEVTGEN_TAUOLA:BOOL=$USETAUOLA -DTAUOLAPP_ROOT_DIR:PATH=$INSTALL_PREFIX \
	+ -DEVTGEN_SHERPA:BOOL=$USESHERPA -DSHERPA_ROOT_DIR:PATH=$INSTALL_PREFIX \
	+ -DEVTGEN_SHERPA3:BOOL=$USESHERPA3 \
	+ $BUILD_BASE/sources/evtgen
	make $BUILDARGS
	make install

	echo Setup done.
	echo To complete, set the Pythia8 data path:
	if [ "$PYTHIAVER" -lt "8200" ]
	then
	echo PYTHIA8DATA=$INSTALL_PREFIX/xmldoc
	else
	echo PYTHIA8DATA=$INSTALL_PREFIX/share/Pythia8/xmldoc
	fi

	echo If installation fully successful you can remove the temporary build area $BUILD_BASE
	cd $BUILD_BASE
	diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
	index fb95832..1f140e6 100644
	--- a/src/CMakeLists.txt
	+++ b/src/CMakeLists.txt
	@@ -1,292 +1,307 @@

	########################################################################
	# 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/>. #
	########################################################################


	# Use glob to find the sources for the main and external libraries
	file(GLOB EVTGEN_SOURCES
	${CMAKE_CURRENT_SOURCE_DIR}/*.cpp
	${CMAKE_CURRENT_SOURCE_DIR}/EvtGenBase/*.cpp
	${CMAKE_CURRENT_SOURCE_DIR}/EvtGenModels/*.cpp
	)

	file(GLOB EVTGEN_EXTERNAL_SOURCES
	${CMAKE_CURRENT_SOURCE_DIR}/EvtGenExternal/*.cpp
	)

	# Add the main EvtGen library...
	# First build all the object files (with PIC enabled) so they can be used to
	# build both the shared and static libs
	add_library(objlib OBJECT ${EVTGEN_SOURCES})
	set_target_properties(objlib PROPERTIES POSITION_INDEPENDENT_CODE 1)
	target_include_directories(objlib PRIVATE ${CMAKE_SOURCE_DIR})
	if(EVTGEN_HEPMC3)
	target_compile_definitions(objlib PRIVATE EVTGEN_HEPMC3)
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib SYSTEM PRIVATE ${HEPMC3_INCLUDE_DIR})
	else()
	target_include_directories(objlib PRIVATE ${HEPMC3_INCLUDE_DIR})
	endif()
	else()
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib SYSTEM PRIVATE ${HEPMC2_INCLUDE_DIR})
	else()
	target_include_directories(objlib PRIVATE ${HEPMC2_INCLUDE_DIR})
	endif()
	endif()

	# Now make the shared library from the object files
	add_library(EvtGen SHARED $<TARGET_OBJECTS:objlib>)
	set_target_properties(EvtGen PROPERTIES OUTPUT_NAME EvtGen)
	set_target_properties(EvtGen PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR} )
	set_target_properties(EvtGen PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR})
	target_include_directories(EvtGen PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}> $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
	if(EVTGEN_HEPMC3)
	target_compile_definitions(EvtGen PUBLIC EVTGEN_HEPMC3)
	target_include_directories(EvtGen PUBLIC ${HEPMC3_INCLUDE_DIR})
	target_link_libraries(EvtGen PUBLIC ${HEPMC3_LIB} ${HEPMC3_SEARCH_LIB})
	else()
	target_include_directories(EvtGen PUBLIC ${HEPMC2_INCLUDE_DIR})
	target_link_libraries(EvtGen PUBLIC ${HEPMC2_LIBRARIES})
	endif()

	# Now make the static library from the object files
	add_library(EvtGenStatic STATIC $<TARGET_OBJECTS:objlib>)
	set_target_properties(EvtGenStatic PROPERTIES OUTPUT_NAME EvtGen)
	set_target_properties(EvtGenStatic PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR})
	target_include_directories(EvtGenStatic PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}> $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
	if(EVTGEN_HEPMC3)
	target_compile_definitions(EvtGenStatic PUBLIC EVTGEN_HEPMC3)
	target_include_directories(EvtGenStatic PUBLIC ${HEPMC3_INCLUDE_DIR})
	target_link_libraries(EvtGenStatic PUBLIC ${HEPMC3_LIB} ${HEPMC3_SEARCH_LIB})
	else()
	target_include_directories(EvtGenStatic PUBLIC ${HEPMC2_INCLUDE_DIR})
	target_link_libraries(EvtGenStatic PUBLIC ${HEPMC2_LIBRARIES})
	endif()



	# Add the EvtGenExternal library (if required)...
	if( EVTGEN_PYTHIA OR EVTGEN_PHOTOS OR EVTGEN_TAUOLA OR EVTGEN_SHERPA)
	# Again, first build all the object files (with PIC enabled) so they can be
	# used to build both the shared and static libs
	add_library(objlib_ext OBJECT ${EVTGEN_EXTERNAL_SOURCES})
	set_target_properties(objlib_ext PROPERTIES POSITION_INDEPENDENT_CODE 1)
	target_include_directories(objlib_ext PRIVATE ${CMAKE_SOURCE_DIR})
	if(EVTGEN_HEPMC3)
	target_compile_definitions(objlib_ext PRIVATE EVTGEN_HEPMC3)
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib_ext SYSTEM PRIVATE ${HEPMC3_INCLUDE_DIR})
	else()
	target_include_directories(objlib_ext PRIVATE ${HEPMC3_INCLUDE_DIR})
	endif()
	else()
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib_ext SYSTEM PRIVATE ${HEPMC2_INCLUDE_DIR})
	else()
	target_include_directories(objlib_ext PRIVATE ${HEPMC2_INCLUDE_DIR})
	endif()
	endif()
	if(EVTGEN_PYTHIA)
	target_compile_definitions(objlib_ext PRIVATE EVTGEN_PYTHIA)
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib_ext SYSTEM PRIVATE ${PYTHIA8_INCLUDE_DIRS})
	else()
	target_include_directories(objlib_ext PRIVATE ${PYTHIA8_INCLUDE_DIRS})
	endif()
	endif()
	if(EVTGEN_PHOTOS)
	target_compile_definitions(objlib_ext PRIVATE EVTGEN_PHOTOS)
	if(Photos++_pp_FOUND)
	target_compile_definitions(objlib_ext PRIVATE EVTGEN_PHOTOS_NEWLIBS)
	endif()
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib_ext SYSTEM PRIVATE ${Photos++_INCLUDE_DIRS})
	else()
	target_include_directories(objlib_ext PRIVATE ${Photos++_INCLUDE_DIRS})
	endif()
	endif()
	if(EVTGEN_TAUOLA)
	target_compile_definitions(objlib_ext PRIVATE EVTGEN_TAUOLA)
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib_ext SYSTEM PRIVATE ${Tauola++_INCLUDE_DIRS})
	else()
	target_include_directories(objlib_ext PRIVATE ${Tauola++_INCLUDE_DIRS})
	endif()
	endif()
	if(EVTGEN_SHERPA)
	target_compile_definitions(objlib_ext PRIVATE EVTGEN_SHERPA)
	- if( ${SHERPA_VERSION} VERSION_GREATER 2.2.15 )
	- target_compile_definitions(objlib_ext PRIVATE EVTGEN_SHERPA_USEMACRO)
	+ if(EVTGEN_SHERPA3)
	+ target_compile_definitions(objlib_ext PRIVATE EVTGEN_SHERPA_USEMACRO EVTGEN_SHERPA3)
	+ get_target_property(SHERPA_INCLUDE_DIRS SHERPA::SherpaMain INTERFACE_INCLUDE_DIRECTORIES)
	+ else()
	+ if( ${SHERPA_VERSION} VERSION_GREATER 2.2.15 )
	+ target_compile_definitions(objlib_ext PRIVATE EVTGEN_SHERPA_USEMACRO)
	+ endif()
	endif()
	if (EVTGEN_SUPPRESS_EXTERNAL_WARNINGS)
	target_include_directories(objlib_ext SYSTEM PRIVATE ${SHERPA_INCLUDE_DIRS})
	else()
	target_include_directories(objlib_ext PRIVATE ${SHERPA_INCLUDE_DIRS})
	endif()
	endif()

	# Now make the shared library from the object files
	add_library(EvtGenExternal SHARED $<TARGET_OBJECTS:objlib_ext>)
	set_target_properties(EvtGenExternal PROPERTIES OUTPUT_NAME EvtGenExternal)
	set_target_properties(EvtGenExternal PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR} )
	set_target_properties(EvtGenExternal PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR})
	target_include_directories(EvtGenExternal PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}> $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
	if(EVTGEN_PYTHIA)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_PYTHIA)
	target_include_directories(EvtGenExternal PUBLIC ${PYTHIA8_INCLUDE_DIRS})
	target_link_libraries(EvtGenExternal PUBLIC ${PYTHIA8_LIBRARIES})
	endif()
	if(EVTGEN_SHERPA)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_SHERPA)
	- if( ${SHERPA_VERSION} VERSION_GREATER 2.2.15 )
	- target_compile_definitions(EvtGenExternal PRIVATE EVTGEN_SHERPA_USEMACRO)
	+ if(EVTGEN_SHERPA3)
	+ target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_SHERPA_USEMACRO EVTGEN_SHERPA3)
	+ target_link_libraries(EvtGenExternal PUBLIC SHERPA::All)
	+ else()
	+ if( ${SHERPA_VERSION} VERSION_GREATER 2.2.15 )
	+ target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_SHERPA_USEMACRO)
	+ endif()
	+ target_include_directories(EvtGenExternal PUBLIC ${SHERPA_INCLUDE_DIRS})
	+ target_link_libraries(EvtGenExternal PUBLIC ${SHERPA_LIBRARIES})
	endif()
	- target_include_directories(EvtGenExternal PUBLIC ${SHERPA_INCLUDE_DIRS})
	- target_link_libraries(EvtGenExternal PUBLIC ${SHERPA_LIBRARIES})
	endif()
	if(EVTGEN_HEPMC3)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_HEPMC3)
	target_include_directories(EvtGenExternal PUBLIC ${HEPMC3_INCLUDE_DIR})
	target_link_libraries(EvtGenExternal PUBLIC ${HEPMC3_LIB} ${HEPMC3_SEARCH_LIB})
	if(EVTGEN_PHOTOS)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_PHOTOS)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_PHOTOS_NEWLIBS)
	target_include_directories(EvtGenExternal PUBLIC ${Photos++_INCLUDE_DIRS})
	# From version 3.64 Photos has HepMC3 support
	target_link_libraries(EvtGenExternal PUBLIC ${Photos++_pp_LIBRARY} ${Photos++_ppHepMC3_LIBRARY})
	endif()
	if(EVTGEN_TAUOLA)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_TAUOLA)
	target_include_directories(EvtGenExternal PUBLIC ${Tauola++_INCLUDE_DIRS})
	# From version 1.1.8 Tauola has HepMC3 support
	target_link_libraries(EvtGenExternal PUBLIC ${Tauola++_CxxInterface_LIBRARY} ${Tauola++_Fortran_LIBRARY} ${Tauola++_HepMC3_LIBRARY})
	endif()
	else()
	target_include_directories(EvtGenExternal PUBLIC ${HEPMC2_INCLUDE_DIR})
	target_link_libraries(EvtGenExternal PUBLIC ${HEPMC2_LIBRARIES})
	if(EVTGEN_PHOTOS)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_PHOTOS)
	target_include_directories(EvtGenExternal PUBLIC ${Photos++_INCLUDE_DIRS})
	# Photos has different library structures for versions before and after 3.58
	# so we need to check which one we have
	if(Photos++_pp_FOUND AND Photos++_ppHepMC_FOUND)
	message(STATUS "EvtGen: PHOTOS has pp and ppHepMC components")
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_PHOTOS_NEWLIBS)
	target_link_libraries(EvtGenExternal PUBLIC ${Photos++_pp_LIBRARY} ${Photos++_ppHepMC_LIBRARY})
	elseif(Photos++_CxxInterface_FOUND AND Photos++_Fortran_FOUND)
	message(STATUS "EvtGen: PHOTOS has CxxInterface and Fortran components")
	target_link_libraries(EvtGenExternal PUBLIC ${Photos++_CxxInterface_LIBRARY} ${Photos++_Fortran_LIBRARY})
	else()
	message(ERROR "EvtGen: PHOTOS lacks pp+ppHepMC or CxxInterface+Fortran components")
	endif()
	endif()
	if(EVTGEN_TAUOLA)
	target_compile_definitions(EvtGenExternal PUBLIC EVTGEN_TAUOLA)
	target_include_directories(EvtGenExternal PUBLIC ${Tauola++_INCLUDE_DIRS})
	target_link_libraries(EvtGenExternal PUBLIC ${Tauola++_CxxInterface_LIBRARY} ${Tauola++_Fortran_LIBRARY})
	# Older versions of Tauola don't have the HepMC component, the HepMC2 interface is in CxxInterface
	if(Tauola++_HepMC_FOUND)
	target_link_libraries(EvtGenExternal PUBLIC ${Tauola++_HepMC_LIBRARY})
	endif()
	endif()
	endif()

	# Now make the static library from the object files
	add_library(EvtGenExternalStatic STATIC $<TARGET_OBJECTS:objlib_ext>)
	set_target_properties(EvtGenExternalStatic PROPERTIES OUTPUT_NAME EvtGenExternal)
	set_target_properties(EvtGenExternalStatic PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR})
	target_include_directories(EvtGenExternalStatic PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}> $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
	if(EVTGEN_PYTHIA)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_PYTHIA)
	target_include_directories(EvtGenExternalStatic PUBLIC ${PYTHIA8_INCLUDE_DIRS})
	target_link_libraries(EvtGenExternalStatic PUBLIC ${PYTHIA8_LIBRARIES})
	endif()
	if(EVTGEN_SHERPA)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_SHERPA)
	- if( ${SHERPA_VERSION} VERSION_GREATER 2.2.15 )
	- target_compile_definitions(EvtGenExternalStatic PRIVATE EVTGEN_SHERPA_USEMACRO)
	+ if(EVTGEN_SHERPA3)
	+ target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_SHERPA_USEMACRO EVTGEN_SHERPA3)
	+ target_link_libraries(EvtGenExternalStatic PUBLIC SHERPA::All)
	+ else()
	+ if( ${SHERPA_VERSION} VERSION_GREATER 2.2.15 )
	+ target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_SHERPA_USEMACRO)
	+ endif()
	+ target_include_directories(EvtGenExternalStatic PUBLIC ${SHERPA_INCLUDE_DIRS})
	+ target_link_libraries(EvtGenExternalStatic PUBLIC ${SHERPA_LIBRARIES})
	endif()
	- target_include_directories(EvtGenExternalStatic PUBLIC ${SHERPA_INCLUDE_DIRS})
	- target_link_libraries(EvtGenExternalStatic PUBLIC ${SHERPA_LIBRARIES})
	endif()
	if(EVTGEN_HEPMC3)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_HEPMC3)
	target_include_directories(EvtGenExternalStatic PUBLIC ${HEPMC3_INCLUDE_DIR})
	target_link_libraries(EvtGenExternalStatic PUBLIC ${HEPMC3_LIB} ${HEPMC3_SEARCH_LIB})
	if(EVTGEN_PHOTOS)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_PHOTOS)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_PHOTOS_NEWLIBS)
	target_include_directories(EvtGenExternalStatic PUBLIC ${Photos++_INCLUDE_DIRS})
	# From version 3.64 Photos has HepMC3 support
	target_link_libraries(EvtGenExternalStatic PUBLIC ${Photos++_pp_LIBRARY} ${Photos++_ppHepMC3_LIBRARY})
	endif()
	if(EVTGEN_TAUOLA)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_TAUOLA)
	target_include_directories(EvtGenExternalStatic PUBLIC ${Tauola++_INCLUDE_DIRS})
	# From version 1.1.8 Tauola has HepMC3 support
	target_link_libraries(EvtGenExternalStatic PUBLIC ${Tauola++_CxxInterface_LIBRARY} ${Tauola++_Fortran_LIBRARY} ${Tauola++_HepMC3_LIBRARY})
	endif()
	else()
	target_include_directories(EvtGenExternalStatic PUBLIC ${HEPMC2_INCLUDE_DIR})
	target_link_libraries(EvtGenExternalStatic PUBLIC ${HEPMC2_LIBRARIES})
	if(EVTGEN_PHOTOS)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_PHOTOS)
	target_include_directories(EvtGenExternalStatic PUBLIC ${Photos++_INCLUDE_DIRS})
	# Photos has different library structures for versions before and after 3.58
	# so we need to check which one we have
	if(Photos++_pp_FOUND AND Photos++_ppHepMC_FOUND)
	#message(STATUS "EvtGen: PHOTOS has pp and ppHepMC components")
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_PHOTOS_NEWLIBS)
	target_link_libraries(EvtGenExternalStatic PUBLIC ${Photos++_pp_LIBRARY} ${Photos++_ppHepMC_LIBRARY})
	elseif(Photos++_CxxInterface_FOUND AND Photos++_Fortran_FOUND)
	#message(STATUS "EvtGen: PHOTOS has CxxInterface and Fortran components")
	target_link_libraries(EvtGenExternalStatic PUBLIC ${Photos++_CxxInterface_LIBRARY} ${Photos++_Fortran_LIBRARY})
	else()
	message(ERROR "EvtGen: PHOTOS lacks pp+ppHepMC or CxxInterface+Fortran components")
	endif()
	endif()
	if(EVTGEN_TAUOLA)
	target_compile_definitions(EvtGenExternalStatic PUBLIC EVTGEN_TAUOLA)
	target_include_directories(EvtGenExternalStatic PUBLIC ${Tauola++_INCLUDE_DIRS})
	target_link_libraries(EvtGenExternalStatic PUBLIC ${Tauola++_CxxInterface_LIBRARY} ${Tauola++_Fortran_LIBRARY})
	# Older versions of Tauola don't have the HepMC component, the HepMC2 interface is in CxxInterface
	if(Tauola++_HepMC_FOUND)
	target_link_libraries(EvtGenExternalStatic PUBLIC ${Tauola++_HepMC_LIBRARY})
	endif()
	endif()
	endif()
	endif()

	# Install the libraries
	install(
	TARGETS EvtGen EvtGenStatic
	EXPORT "EvtGenTargets"
	LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
	ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
	)

	if( EVTGEN_PYTHIA OR EVTGEN_PHOTOS OR EVTGEN_TAUOLA OR EVTGEN_SHERPA)
	install(
	TARGETS EvtGenExternal EvtGenExternalStatic
	EXPORT "EvtGenTargets"
	LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
	ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
	)
	endif()
	diff --git a/src/EvtGenExternal/EvtSherpaPhotons.cpp b/src/EvtGenExternal/EvtSherpaPhotons.cpp
	index dcecee1..e34f0f9 100644
	--- a/src/EvtGenExternal/EvtSherpaPhotons.cpp
	+++ b/src/EvtGenExternal/EvtSherpaPhotons.cpp
	@@ -1,303 +1,369 @@

	/***********************************************************************
	* 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_SHERPA

	#include "EvtGenExternal/EvtSherpaPhotons.hh"

	#include "EvtGenBase/EvtPhotonParticle.hh"
	#include "EvtGenBase/EvtRandom.hh"
	#include "EvtGenBase/EvtReport.hh"
	#include "EvtGenBase/EvtVector4R.hh"

	#include "ATOOLS/Org/MyStrStream.H"
	+#include "ATOOLS/Phys/Flavour.H"
	+#ifdef EVTGEN_SHERPA3
	+#include "ATOOLS/Phys/KF_Table.H"
	+#else
	+#include "ATOOLS/Org/Data_Reader.H"
	+#endif
	+#include "ATOOLS/Phys/Blob.H"
	#include "ATOOLS/Phys/Particle.H"
	-#include "SHERPA/Initialization/Initialization_Handler.H"
	-#include "SHERPA/SoftPhysics/Soft_Photon_Handler.H"

	#include <cstring>
	#include <iostream>
	#include <vector>

	using std::endl;

	+std::mutex EvtSherpaPhotons::m_sherpa_mutex;
	std::unique_ptr<SHERPA::Sherpa> EvtSherpaPhotons::m_sherpaGen;
	+std::unique_ptr<PHOTONS::Photons> EvtSherpaPhotons::m_photonsGen;
	bool EvtSherpaPhotons::m_initialised = false;
	-std::mutex EvtSherpaPhotons::m_sherpa_mutex;

	EvtSherpaPhotons::EvtSherpaPhotons( const bool useEvtGenRandom,
	const double infraredCutOff, const int mode,
	const int useME ) :
	m_useEvtGenRandom{ useEvtGenRandom },
	m_infraredCutOff{ infraredCutOff },
	m_mode{ mode },
	m_useME{ useME }
	{
	}

	void EvtSherpaPhotons::initialise()
	{
	// Sherpa initialisation is not thread safe, so we mutex it
	const std::lock_guard<std::mutex> lock( m_sherpa_mutex );

	m_gammaId = EvtPDL::getId( m_photonType );
	m_gammaPDG = EvtPDL::getStdHep( m_gammaId );

	if ( m_initialised ) {
	return;
	}

	+#ifdef EVTGEN_SHERPA3
	+ const std::string equal_string = ":";
	+ const std::string yfs_prefix = "YFS:\n";
	+ const std::string mode_string = m_mode == 2
	+ ? "Full"
	+ : ( m_mode == 1 ? "Soft" : "None" );
	+
	+ // Beam parameters are necessary in Sherpa 3 for initialisation
	+ // Otherwise m_s (CMS energy) is nan in YFS_Base and creates SegFault
	+ // This has no impact on PHOTONS++. USING LHC values below.
	+ m_configs.push_back( "BEAMS:2212" );
	+ m_configs.push_back( "BEAM_ENERGIES:7000" );
	+
	+#else
	+ const std::string equal_string = "=";
	+ const std::string yfs_prefix = "YFS_";
	+ const std::string mode_string = ATOOLS::ToString( m_mode );
	+#endif
	+
	+ // INIT_ONLY=6 intialises the Sherpa objects without launching simulation.
	+ m_configs.push_back( "INIT_ONLY" + equal_string + "6" );
	+
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< "Setting up Sherpa PHOTONS++ generator for FSR." << endl;

	- /* Specify whether we are going to use EvtGen's random number generator
	+ /* Specify whether we are going to use EvtGen's random number generator
	* (via EvtSherpaRandom interface) for Sherpa's PHOTONS++ simulation. */
	if ( m_useEvtGenRandom ) {
	- m_configs.push_back( "SHERPA_LDADD=EvtGenExternal" );
	- m_configs.push_back( "EXTERNAL_RNG=EvtSherpaRandom" );
	+ m_configs.push_back( "SHERPA_LDADD" + equal_string + "EvtGenExternal" );
	+ m_configs.push_back( "EXTERNAL_RNG" + equal_string + "EvtSherpaRandom" );
	}

	/*
	* Internal PHOTONS++ settings.
	- * Documentation at
	+ * Documentation at
	* https://sherpa.hepforge.org/doc/Sherpa.html#QED-Corrections
	*/

	/*
	- * YFS_IR_CUTOFF sets the infrared cut-off which serves as a
	+ * YFS_IR_CUTOFF sets the infrared cut-off which serves as a
	* minimum photon energy in this frame for explicit photon generation.
	* The default is YFS_IR_CUTOFF = 1E-3 (GeV) but we set it to 1E-7 to equalize
	* with the cutoff used for PHOTOS.
	*/
	- m_configs.push_back( "YFS_IR_CUTOFF=" + ATOOLS::ToString( m_infraredCutOff ) );
	+ m_configs.push_back( yfs_prefix + "IR_CUTOFF" + equal_string +
	+ ATOOLS::ToString( m_infraredCutOff ) );

	/*
	- * The keyword YFS_MODE = [0,1,2] determines the mode of operation.
	- * YFS_MODE = 0 switches Photons off.
	- * YFS_MODE = 1 sets the mode to "soft only", meaning soft emissions will be
	- * treated correctly to all orders but no hard emission corrections will be included.
	- * YFS_MODE = 2 these hard emission corrections will also be included up to
	- * first order in alpha_QED. This is the default setting in Sherpa.
	+ * The keyword YFS_MODE = [0,1,2] determines the mode of operation.
	+ * YFS_MODE = 0 switches Photons off.
	+ * YFS_MODE = 1 sets the mode to "soft only", meaning soft emissions will be
	+ * treated correctly to all orders but no hard emission corrections will be included.
	+ * YFS_MODE = 2 these hard emission corrections will also be included up to
	+ * first order in alpha_QED. This is the default setting in Sherpa.
	*/
	-
	- m_configs.push_back( "YFS_MODE=" + ATOOLS::ToString( m_mode ) );
	+ m_configs.push_back( yfs_prefix + "MODE" + equal_string + mode_string );

	/*
	- * The switch YFS_USE_ME = [0,1] tells Photons how to correct hard emissions to
	- * first order in alpha_QED. If YFS_USE_ME = 0, then Photons will use collinearly
	- * approximated real emission matrix elements. Virtual emission matrix elements of
	- * order alpha_QED are ignored. If, however, YFS_USE_ME=1, then exact real and/or
	- * virtual emission matrix elements are used wherever possible.
	- * These are presently available for V->FF, V->SS, S->FF, S->SS, S->Slnu,
	- * S->Vlnu type decays, Z->FF decays and leptonic tau and W decays. For all other
	- * decay types general collinearly approximated matrix elements are used.
	- * In both approaches all hadrons are treated as point-like objects.
	+ * The switch YFS_USE_ME = [0,1] tells Photons how to correct hard emissions to
	+ * first order in alpha_QED. If YFS_USE_ME = 0, then Photons will use collinearly
	+ * approximated real emission matrix elements. Virtual emission matrix elements of
	+ * order alpha_QED are ignored. If, however, YFS_USE_ME=1, then exact real and/or
	+ * virtual emission matrix elements are used wherever possible.
	+ * These are presently available for V->FF, V->SS, S->FF, S->SS, S->Slnu,
	+ * S->Vlnu type decays, Z->FF decays and leptonic tau and W decays. For all other
	+ * decay types general collinearly approximated matrix elements are used.
	+ * In both approaches all hadrons are treated as point-like objects.
	* The default setting is YFS_USE_ME = 1. This switch is only effective if YFS_MODE = 2.
	*/
	- m_configs.push_back( "YFS_USE_ME=" + ATOOLS::ToString( m_useME ) );
	-
	- // TODO: Virtual photon splitting into l+l- will be part of ME corrections from Sherpa 3.0.0 on,
	- // Once released, this should be switched off by default or
	- // taken into account while retrieving the event.
	+ m_configs.push_back( yfs_prefix + "USE_ME" + equal_string +
	+ ATOOLS::ToString( m_useME ) );
	+
	+#ifdef EVTGEN_SHERPA3
	+ // Virtual photon splitting is part of ME corrections from Sherpa 3.0.0 on,
	+ // We switch this off below.
	+ // Support for photon splitting should be implemented for all FSR plug-ins,
	+ // taking into account produced pairs while retrieving the event.
	+ m_configs.push_back( yfs_prefix + "PHOTON_SPLITTER_MODE" + equal_string +
	+ ATOOLS::ToString( 0 ) );
	+#endif

	// Set up run-time arguments for Sherpa.
	std::vector<char*> argv = this->addParameters();

	// Create instance and initialise Sherpa.
	+#ifdef EVTGEN_SHERPA3
	+ m_sherpaGen = std::make_unique<SHERPA::Sherpa>( argv.size(), &argv[0] );
	+ m_sherpaGen->InitializeTheRun();
	+ // Create the PHOTONS++ instance
	+ m_photonsGen = std::make_unique<PHOTONS::Photons>();
	+#else
	m_sherpaGen = std::make_unique<SHERPA::Sherpa>();
	m_sherpaGen->InitializeTheRun( argv.size(), &argv[0] );
	+ // Create the PHOTONS++ instance
	+ // Needs to create and pass a Data_reader object as in the Soft_Photon_Handler
	+ std::unique_ptr<ATOOLS::Data_Reader> dataread =
	+ std::make_unique<ATOOLS::Data_Reader>( " ", ";", "!", "=" );
	+ m_photonsGen = std::make_unique<PHOTONS::Photons>( dataread.get() );
	+#endif

	this->updateParticleLists();

	m_initialised = true;
	}

	std::vector<char*> EvtSherpaPhotons::addParameters()
	{
	std::vector<char*> argv;
	argv.reserve( m_configs.size() );

	for ( auto& config : m_configs ) {
	if ( config != "Sherpa" ) {
	EvtGenReport( EVTGEN_INFO, "EvtGen" )
	<< " EvtSherpaPhotons::initialise: Setting parameter '"
	<< config << "'" << endl;
	}
	argv.push_back( config.data() );
	}

	return argv;
	}

	void EvtSherpaPhotons::updateParticleLists()
	{
	- /* Use the EvtGen decay and pdl tables (decay/user.dec and evt.pdl)
	- * to update Sherpa's KF_Table which contains all defined particles
	- * and their properties.
	+ /* Use the EvtGen decay and pdl tables (decay/user.dec and evt.pdl)
	+ * to update Sherpa's KF_Table which contains all defined particles
	+ * and their properties.
	* Loop over all entries in the EvtPDL particle data table.
	*/
	const int nPDL = EvtPDL::entries();

	for ( int iPDL = 0; iPDL < nPDL; iPDL++ ) {
	const EvtId particleId = EvtPDL::getEntry( iPDL );

	// Sherpa and EvtGen should use the same PDG codes for particles
	const int PDGCode = EvtPDL::getStdHep( particleId );
	const long unsigned int particleCode = std::abs( PDGCode );

	// If particle is found in Sherpa's KF_Table, then continue.
	// Updating the pole mass is not necessary because we provide the 4-momenta
	// and the final mass for final-state radiation generation.
	ATOOLS::KFCode_ParticleInfo_Map::const_iterator kf_it(
	ATOOLS::s_kftable.find( particleCode ) );
	if ( kf_it != ATOOLS::s_kftable.end() ) {
	continue;
	}

	// If the PDG is negative, use the charge conjugated ID because
	// the KF_Table is defined only for positive PDG numbers.
	const EvtId particleIdConj = EvtPDL::chargeConj( particleId );

	// Get mass, width, charge and spin
	const double mass = EvtPDL::getMeanMass( particleId );
	const double width = EvtPDL::getWidth( particleId );
	const int charge = PDGCode < 0 ? EvtPDL::chg3( particleIdConj )
	: EvtPDL::chg3( particleId );
	const int spin = EvtSpinType::getSpin2(
	EvtPDL::getSpinType( particleId ) );

	// Use as reference the name of the particle (not anti-particle).
	const std::string idName = PDGCode < 0 ? particleIdConj.getName()
	: particleId.getName();

	// Create new entry in the KF_Table. The 0 and 1 in the constructor below
	// stand for m_on = 0 and m_stable = 1. These properties are irrelevant
	// for final-state radiation and set to default values.
	+#ifdef EVTGEN_SHERPA3
	+ const double radius{ 0.0 };
	+ ATOOLS::s_kftable[particleCode] =
	+ new ATOOLS::Particle_Info( particleCode, mass, radius, width,
	+ charge, spin, 0, 1, idName, idName );
	+#else
	ATOOLS::s_kftable[particleCode] = new ATOOLS::Particle_Info(
	particleCode, mass, width, charge, spin, 0, 1, idName, idName );
	+#endif
	}
	}

	void EvtSherpaPhotons::doRadCorr( EvtParticle* theParticle )
	{
	if ( !theParticle ) {
	return;
	}

	/* Skip running FSR if the particle has no daughters, since we can't add FSR.
	* Also skip FSR if the particle has too many daughters (>= 10) as done for PHOTOS.
	*/
	const int nDaughters( theParticle->getNDaug() );
	if ( nDaughters == 0 \|\| nDaughters >= 10 ) {
	return;
	}

	// Create a blob.
	std::unique_ptr<ATOOLS::Blob> blob = std::make_unique<ATOOLS::Blob>();

	- // Tell Sherpa that the blob needs FSR (that is extra QED)
	- blob->SetStatus( ATOOLS::blob_status::needs_extraQED );
	-
	// Create mother particle and add it to blob
	ATOOLS::Flavour mother_flav( EvtPDL::getStdHep( theParticle->getId() ) );
	mother_flav.SetStable( false );
	const double motherM0 = theParticle->mass();
	const ATOOLS::Vec4D mother_mom( motherM0, 0., 0., 0. );

	// Add mother to blob with m_number=-1. The m_number will start from 0 for daughters.
	- ATOOLS::Particle* mother_part = new ATOOLS::Particle( -1, mother_flav,
	- mother_mom );
	+ std::unique_ptr<ATOOLS::Particle> mother_part =
	+ std::make_unique<ATOOLS::Particle>( -1, mother_flav, mother_mom );
	mother_part->SetFinalMass( motherM0 );
	mother_part->SetStatus( ATOOLS::part_status::decayed );
	// Set m_info="I" (initial) for mother. The m_info will be "F" (final) for daughters.
	mother_part->SetInfo( 'I' );

	- blob->AddToInParticles( mother_part );
	+ blob->AddToInParticles( mother_part.get() );
	+
	+ mother_part.release();

	// Add daughters to the blob
	for ( int iDaug = 0; iDaug < nDaughters; iDaug++ ) {
	const EvtParticle* theDaughter = theParticle->getDaug( iDaug );
	ATOOLS::Flavour daughter_flav( EvtPDL::getStdHep( theDaughter->getId() ) );
	daughter_flav.SetStable( true );

	const double daugE = theDaughter->getP4().get( 0 );
	const double daugPx = theDaughter->getP4().get( 1 );
	const double daugPy = theDaughter->getP4().get( 2 );
	const double daugPz = theDaughter->getP4().get( 3 );
	const double daugM0 = theDaughter->mass();

	const ATOOLS::Vec4D daughter_mom( daugE, daugPx, daugPy, daugPz );

	- ATOOLS::Particle* daughter_part =
	- new ATOOLS::Particle( iDaug, daughter_flav, daughter_mom );
	+ std::unique_ptr<ATOOLS::Particle> daughter_part =
	+ std::make_unique<ATOOLS::Particle>( iDaug, daughter_flav,
	+ daughter_mom );
	daughter_part->SetFinalMass( daugM0 );

	daughter_part->SetStatus( ATOOLS::part_status::active );
	daughter_part->SetInfo( 'F' );

	- blob->AddToOutParticles( daughter_part );
	+ blob->AddToOutParticles( daughter_part.get() );
	+
	+ daughter_part.release();
	}

	{
	// AddRadiation function is not yet thread_safe, so we mutex it
	const std::lock_guard<std::mutex> lock( m_sherpa_mutex );

	- const SHERPA::Initialization_Handler* inithandler =
	- m_sherpaGen->GetInitHandler();
	- SHERPA::Soft_Photon_Handler* softphotonhandler =
	- inithandler->GetSoftPhotonHandler();
	-
	// Simulate the radiation
	- softphotonhandler->AddRadiation( blob.get() );
	+ m_photonsGen->AddRadiation( blob.get() );
	}

	// Get number of final-state particles
	const int nFinal( blob->NOutP() );

	// Retrieve the event from Sherpa blob if FSR photons were added
	if ( nFinal > nDaughters ) {
	for ( int iLoop = 0; iLoop < nFinal; iLoop++ ) {
	const ATOOLS::Particle* outParticle = blob->OutParticle( iLoop );
	const ATOOLS::Vec4D daughter_mom = outParticle->Momentum();
	const long int pdgId = outParticle->Flav();
	const EvtVector4R newP4( daughter_mom[0], daughter_mom[1],
	daughter_mom[2], daughter_mom[3] );
	const char daugInfo = outParticle->Info();

	if ( iLoop < nDaughters ) {
	// Update momenta of initial particles
	EvtParticle* daugParticle = theParticle->getDaug( iLoop );
	if ( daugParticle ) {
	daugParticle->setP4WithFSR( newP4 );
	}
	} else if ( pdgId == m_gammaPDG && daugInfo == 'S' ) {
	// Add FSR photons. PHOTONS flags them with info 'S'.
	// Create a new photon particle and add it to the list of daughters
	EvtPhotonParticle* gamma = new EvtPhotonParticle();
	gamma->init( m_gammaId, newP4 );
	gamma->setFSRP4toZero();
	gamma->setAttribute( "FSR", 1 );
	gamma->addDaug( theParticle );
	}
	}
	}
	}

	+void EvtSherpaPhotons::finalise()
	+{
	+ // Mutex this since Sherpa instances are static
	+ const std::lock_guard<std::mutex> lock( m_sherpa_mutex );
	+
	+ if ( !m_initialised ) {
	+ return;
	+ }
	+
	+ // We explicitly destroy the static instances to make sure that this occurs in the right order.
	+ // This is needed due to the static destruction order fiasco.
	+ m_photonsGen.reset();
	+ m_sherpaGen.reset();
	+ m_initialised = false;
	+}
	+
	#endif
	diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
	index 04446b5..c1751b7 100644
	--- a/test/testDecayModel.cc
	+++ b/test/testDecayModel.cc
	@@ -1,1832 +1,1844 @@

	/***********************************************************************
	* 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/>. *
	***********************************************************************/

	#include "testDecayModel.hh"

	#include "EvtGen/EvtGen.hh"

	#include "EvtGenBase/EvtAbsRadCorr.hh"
	#include "EvtGenBase/EvtConst.hh"
	#include "EvtGenBase/EvtDecayBase.hh"
	#include "EvtGenBase/EvtId.hh"
	#include "EvtGenBase/EvtKine.hh"
	#include "EvtGenBase/EvtMTRandomEngine.hh"
	#include "EvtGenBase/EvtPDL.hh"
	#include "EvtGenBase/EvtParticle.hh"
	#include "EvtGenBase/EvtParticleFactory.hh"
	#include "EvtGenBase/EvtRandom.hh"
	#include "EvtGenBase/EvtVector4R.hh"

	#ifdef EVTGEN_EXTERNAL
	#include "EvtGenExternal/EvtExternalGenList.hh"
	+#include "EvtGenExternal/EvtSherpaPhotons.hh"
	#endif

	#include "TROOT.h"

	#include "tbb/blocked_range.h"
	#include "tbb/global_control.h"
	#include "tbb/parallel_reduce.h"

	#include <chrono>
	#include <fstream>
	#include <future>
	#include <iostream>
	#include <list>
	#include <memory>

	using nlohmann::json;

	std::once_flag TestDecayModel::m_createDecFile_threadlock;

	TestHistos::TestHistos( const std::string& parentName, const json& config )
	{
	// Histogram information
	const std::size_t nHistos{ config.size() };

	m_1DhistVect.reserve( nHistos );
	m_2DhistVect.reserve( nHistos );

	for ( const auto& hInfo : config ) {
	const auto varTitle{ hInfo.at( "title" ).get<std::string>() };

	const auto varName{ hInfo.at( "variable" ).get<std::string>() };
	// Integer values that define what particles need to be used
	// for invariant mass combinations or helicity angles etc
	const auto d1{ hInfo.at( "d1" ).get<int>() };
	const auto d2{ hInfo.at( "d2" ).get<int>() };

	const auto nBins{ hInfo.at( "nbins" ).get<int>() };
	const auto xmin{ hInfo.at( "xmin" ).get<double>() };
	const auto xmax{ hInfo.at( "xmax" ).get<double>() };

	std::string histName( varName.c_str() );
	if ( d1 != 0 ) {
	histName += "_";
	histName += std::to_string( d1 );
	}
	if ( d2 != 0 ) {
	histName += "_";
	histName += std::to_string( d2 );
	}

	if ( !hInfo.contains( "variableY" ) ) {
	auto hist = std::make_unique<TH1D>( histName.c_str(),
	varTitle.c_str(), nBins, xmin,
	xmax );
	m_1DhistVect.emplace_back(
	std::make_pair( HistInfo{ varName, d1, d2 }, std::move( hist ) ) );
	} else {
	const auto varNameY{ hInfo.at( "variableY" ).get<std::string>() };
	const auto d1Y{ hInfo.at( "d1Y" ).get<int>() };
	const auto d2Y{ hInfo.at( "d2Y" ).get<int>() };

	const auto nBinsY{ hInfo.at( "nbinsY" ).get<int>() };
	const auto ymin{ hInfo.at( "ymin" ).get<double>() };
	const auto ymax{ hInfo.at( "ymax" ).get<double>() };

	histName += "_";
	histName += varNameY;
	if ( d1Y != 0 ) {
	histName += "_";
	histName += std::to_string( d1Y );
	}
	if ( d2Y != 0 ) {
	histName += "_";
	histName += std::to_string( d2Y );
	}
	auto hist = std::make_unique<TH2D>( histName.c_str(),
	varTitle.c_str(), nBins, xmin,
	xmax, nBinsY, ymin, ymax );
	m_2DhistVect.emplace_back(
	std::make_pair( HistInfo{ varName, d1, d2, varNameY, d1Y, d2Y },
	std::move( hist ) ) );
	}
	}

	// Add a mixed/unmixed histogram
	// Useful for the case where the parent is a neutral K, D or B
	const std::array<std::string, 14> parentsThatMix{
	"B_s0", "anti-B_s0", "B_s0L", "B_s0H", "B0", "anti-B0", "B0L",
	"B0H", "D0", "anti-D0", "K0", "anti-K0", "K_S0", "K_L0" };
	if ( std::find( parentsThatMix.begin(), parentsThatMix.end(), parentName ) !=
	parentsThatMix.end() ) {
	const std::string varTitle{ parentName + " mixed" };
	m_mixedHist = std::make_unique<TH1D>( "mixed", varTitle.c_str(), 2, 0.0,
	2.0 );
	// TODO maybe set bin labels?
	}
	}

	TestHistos::TestHistos( const TestHistos& rhs )
	{
	m_1DhistVect.reserve( rhs.m_1DhistVect.size() );
	for ( auto& [info, hist] : rhs.m_1DhistVect ) {
	auto newHist = std::unique_ptr<TH1>{ static_cast<TH1*>( hist->Clone() ) };
	m_1DhistVect.push_back( std::make_pair( info, std::move( newHist ) ) );
	}

	m_2DhistVect.reserve( rhs.m_2DhistVect.size() );
	for ( auto& [info, hist] : rhs.m_2DhistVect ) {
	auto newHist = std::unique_ptr<TH2>{ static_cast<TH2*>( hist->Clone() ) };
	m_2DhistVect.push_back( std::make_pair( info, std::move( newHist ) ) );
	}

	if ( rhs.m_mixedHist ) {
	m_mixedHist.reset( static_cast<TH1*>( rhs.m_mixedHist->Clone() ) );
	}
	}

	TestHistos::TestHistos( TestHistos&& rhs ) noexcept
	{
	this->swap( rhs );
	}

	TestHistos& TestHistos::operator=( const TestHistos& rhs )
	{
	TestHistos tmp{ rhs };
	this->swap( tmp );
	return *this;
	}

	TestHistos& TestHistos::operator=( TestHistos&& rhs ) noexcept
	{
	this->swap( rhs );
	return *this;
	}

	void TestHistos::swap( TestHistos& rhs ) noexcept
	{
	m_1DhistVect.swap( rhs.m_1DhistVect );
	m_2DhistVect.swap( rhs.m_2DhistVect );
	std::swap( m_mixedHist, rhs.m_mixedHist );
	}

	void TestHistos::add( const TestHistos& rhs )
	{
	// handle the special case where we have been default constructed and the rhs has not
	if ( m_1DhistVect.empty() && m_2DhistVect.empty() && !m_mixedHist ) {
	( *this ) = rhs;
	return;
	}

	// TODO - should really check that the sets of histograms are the same between left and right

	const std::size_t n1DHists{ rhs.m_1DhistVect.size() };
	for ( std::size_t i{ 0 }; i < n1DHists; ++i ) {
	m_1DhistVect[i].second->Add( rhs.m_1DhistVect[i].second.get() );
	}

	const std::size_t n2DHists{ rhs.m_2DhistVect.size() };
	for ( std::size_t i{ 0 }; i < n2DHists; ++i ) {
	m_2DhistVect[i].second->Add( rhs.m_2DhistVect[i].second.get() );
	}

	if ( m_mixedHist && rhs.m_mixedHist ) {
	m_mixedHist->Add( rhs.m_mixedHist.get() );
	}
	}

	void TestHistos::normalise()
	{
	for ( auto& [_, hist] : m_1DhistVect ) {
	const double area{ hist->Integral() };
	if ( area > 0.0 ) {
	hist->Scale( 1.0 / area );
	}
	}
	for ( auto& [_, hist] : m_2DhistVect ) {
	const double area{ hist->Integral() };
	if ( area > 0.0 ) {
	hist->Scale( 1.0 / area );
	}
	}
	if ( m_mixedHist ) {
	const double area{ m_mixedHist->Integral() };
	if ( area > 0.0 ) {
	m_mixedHist->Scale( 1.0 / area );
	}
	}
	}

	void TestHistos::save( TFile* outputFile )
	{
	outputFile->cd();

	for ( auto& [info, hist] : m_1DhistVect ) {
	hist->SetDirectory( outputFile );
	hist->Write();
	hist.release();
	}

	for ( auto& [info, hist] : m_2DhistVect ) {
	hist->SetDirectory( outputFile );
	hist->Write();
	hist.release();
	}

	if ( m_mixedHist ) {
	m_mixedHist->SetDirectory( outputFile );
	m_mixedHist->Write();
	m_mixedHist.release();
	}
	}

	TestDecayModel::TestDecayModel( const json& config ) :
	m_config{ checkMandatoryFields( config )
	? readConfig( config )
	: throw std::runtime_error{
	"ERROR : json does not contain all required fields" } }
	{
	}

	bool TestDecayModel::checkMandatoryFields( const json& config )
	{
	const std::array<std::string, 7> mandatoryFields{ "parent", "daughters",
	"models", "parameters",
	"outfile", "events",
	"histograms" };
	const std::array<std::string, 7> mandatoryHistoFields{
	"title", "variable", "d1", "d2", "nbins", "xmin", "xmax" };
	const std::array<std::string, 6> extra2DHistoFields{ "variableY", "d1Y",
	"d2Y", "nbinsY",
	"ymin", "ymax" };

	bool allMandatoryFields{ true };

	for ( const auto& field : mandatoryFields ) {
	if ( !config.contains( field ) ) {
	std::cerr << "ERROR : json does not contain required field: " << field
	<< std::endl;
	allMandatoryFields = false;
	continue;
	}
	if ( field == "histograms" ) {
	const json& jHistos{ config.at( "histograms" ) };
	for ( const auto& hInfo : jHistos ) {
	for ( const auto& hField : mandatoryHistoFields ) {
	if ( !hInfo.contains( hField ) ) {
	std::cerr
	<< "ERROR : json does not contain required field for histogram definition: "
	<< hField << std::endl;
	allMandatoryFields = false;
	}
	}
	if ( hInfo.contains( extra2DHistoFields[0] ) ) {
	for ( const auto& hField : extra2DHistoFields ) {
	if ( !hInfo.contains( hField ) ) {
	std::cerr
	<< "ERROR : json does not contain required field for 2D histogram definition: "
	<< hField << std::endl;
	allMandatoryFields = false;
	}
	}
	}
	}
	}
	}

	return allMandatoryFields;
	}

	TestConfig TestDecayModel::readConfig( const json& config )
	{
	TestConfig cfg;

	// Get all the mandatory fields first
	cfg.parentName = config.at( "parent" ).get<std::string>();
	cfg.daughterNames = config.at( "daughters" ).get<std::vector<std::string>>();
	cfg.modelNames = config.at( "models" ).get<std::vector<std::string>>();
	cfg.modelParameters =
	config.at( "parameters" ).get<std::vector<std::vector<std::string>>>();
	cfg.nEvents = config.at( "events" ).get<std::size_t>();

	// Histogram information
	cfg.testHistograms = TestHistos{ cfg.parentName, config.at( "histograms" ) };

	// Then check for optional fields, setting default values if not present

	if ( config.contains( "grand_daughters" ) &&
	config.at( "grand_daughters" ).is_array() ) {
	cfg.grandDaughterNames = config.at( "grand_daughters" )
	.get<std::vector<std::vector<std::string>>>();
	}

	if ( config.contains( "extras" ) && config.at( "extras" ).is_array() ) {
	cfg.extraCommands = config.at( "extras" ).get<std::vector<std::string>>();
	}

	// Set the number of threads to use, 1 by default
	cfg.nThreads = 1;
	if ( config.contains( "threads" ) ) {
	cfg.nThreads = config.at( "threads" ).get<std::size_t>();
	}

	// Set the FSR generator. Use PHOTOS by default.
	cfg.fsrGenerator = FSRGenerator::PHOTOS;
	if ( config.contains( "fsr_generator" ) ) {
	cfg.fsrGenerator = config.at( "fsr_generator" ).get<FSRGenerator>();
	}

	// Set the type of threading to use, stdlib by default
	cfg.threadModel = ThreadModel::StdLib;
	if ( config.contains( "thread_model" ) ) {
	cfg.threadModel = config.at( "thread_model" ).get<ThreadModel>();
	}

	// Set the RNG seed base (defaults to zero), to which the event number is added
	cfg.rngSeed = 0;
	if ( config.contains( "rng_seed" ) ) {
	cfg.rngSeed = config.at( "rng_seed" ).get<std::size_t>();
	}

	// Set reference and output file names
	// Insert fsrGenerator name if FSR simulation is not deactivated
	const bool noFSR = std::find( cfg.extraCommands.begin(),
	cfg.extraCommands.end(),
	"noFSR" ) != cfg.extraCommands.end();

	const std::string fileNameEnd =
	noFSR ? ".root" : "_" + to_string( cfg.fsrGenerator ) + ".root";

	const auto outFileStrSize = config.at( "outfile" ).get<std::string>().size() -
	5;

	cfg.outFileName =
	config.at( "outfile" ).get<std::string>().substr( 0, outFileStrSize ) +
	fileNameEnd;
	cfg.refFileName = "Ref/" + cfg.outFileName;
	cfg.decFileName = cfg.outFileName.substr( 0, cfg.outFileName.size() - 5 ) +
	".dec";

	cfg.debugFlag = ( config.contains( "debug_flag" ) &&
	config.at( "debug_flag" ).is_boolean() )
	? config.at( "debug_flag" ).get<bool>()
	: false;

	if ( config.contains( "do_conjugate_decay" ) &&
	config.at( "do_conjugate_decay" ).is_array() ) {
	cfg.doConjDecay =
	config.at( "do_conjugate_decay" ).get<std::vector<bool>>();
	}
	if ( cfg.doConjDecay.size() != cfg.modelNames.size() ) {
	cfg.doConjDecay.resize( cfg.modelNames.size(), false );
	}

	return cfg;
	}

	void TestDecayModel::run()
	{
	TestHistos theHistos;

	const auto start{ std::chrono::steady_clock::now() };

	if ( m_config.nThreads > 1 ) {
	// Run multi-threaded using the specified thread model
	switch ( m_config.threadModel ) {
	case ThreadModel::StdLib:
	theHistos = runStdThreads();
	break;
	case ThreadModel::TBB:
	theHistos = runTBBThreads();
	break;
	}
	} else {
	// Run in the main thread
	theHistos = runDecayBody( 0, m_config.nEvents );
	}

	+#ifdef EVTGEN_SHERPA
	+ // The Sherpa and PHOTONS++ instances are static global instances inside EvtGen
	+ // as they are not thread safe (they are mutexed). From Sherpa 3.0.0 on,
	+ // the internal Sherpa settings are also a static global instance.
	+ // Because of the static destruction order fiasco,
	+ // we need to control the order in which the static instances are destroyed.
	+ // This is done by hand inside the function below. If the function below is not called,
	+ // a SegFault occurs at the end of the program execution when using Sherpa for FSR.
	+ EvtSherpaPhotons::finalise();
	+#endif
	+
	const auto end{ std::chrono::steady_clock::now() };
	const std::chrono::duration<double, std::milli> elapsed_ms{ ( end - start ) };
	const std::chrono::duration<double, std::milli> elapsed_ms_per_event{
	elapsed_ms / m_config.nEvents };

	std::cout << "Took " << elapsed_ms.count() << " ms to generate "
	<< m_config.nEvents << " events using " << m_config.nThreads
	<< " " << to_string( m_config.threadModel ) << " threads ("
	<< elapsed_ms_per_event.count() << " ms per event)" << std::endl;

	// Normalise histograms.
	theHistos.normalise();

	// Compare with reference histograms
	compareHistos( theHistos, m_config.refFileName );

	// Create the root output file and write the histograms to it
	// Only save the mixed/unmixed histogram for neutral K, D, B
	std::unique_ptr<TFile> outFile{
	TFile::Open( m_config.outFileName.c_str(), "recreate" ) };
	theHistos.save( outFile.get() );
	outFile->Close();
	std::cout << "Created output file: " << m_config.outFileName.c_str()
	<< std::endl;
	}

	TestHistos TestDecayModel::runStdThreads() const
	{
	// Determine the number of threads and the number of events per thread
	const std::size_t nThreads{ std::min( m_config.nThreads, m_config.nEvents ) };
	const std::size_t nEventsPerThread{ ( m_config.nEvents % nThreads )
	? m_config.nEvents / nThreads + 1
	: m_config.nEvents / nThreads };

	// Create the store for the results from each thread
	std::vector<std::future<TestHistos>> allHistos;
	allHistos.reserve( nThreads );

	// Launch the threads
	std::size_t firstEvent{ 0 };
	std::size_t nEvents{ nEventsPerThread };

	for ( std::size_t iThread{ 0 }; iThread < nThreads; ++iThread ) {
	// The last thread may need to generate slightly fewer events to get the required total
	if ( ( firstEvent + nEvents ) >= m_config.nEvents ) {
	nEvents = m_config.nEvents - firstEvent;
	}

	std::cout << "Thread " << iThread << " will generate " << nEvents
	<< " events" << std::endl;

	allHistos.emplace_back(
	std::async( std::launch::async, [this, firstEvent, nEvents]()
	{ return runDecayBody( firstEvent, nEvents ); } ) );

	firstEvent += nEvents;
	}

	// Now wait for each thread to finish
	bool complete{ false };
	do {
	// Set the flag to completed and set it back if we find incomplete threads
	complete = true;
	for ( auto& future : allHistos ) {
	auto status = future.wait_for( std::chrono::seconds( 10 ) );
	if ( status != std::future_status::ready ) {
	complete = false;
	}
	}
	} while ( !complete );

	// Accumulate the histograms from all threads
	TestHistos theHistos{ allHistos[0].get() };
	for ( std::size_t iThread{ 1 }; iThread < nThreads; ++iThread ) {
	theHistos.add( allHistos[iThread].get() );
	}

	return theHistos;
	}

	TestHistos TestDecayModel::runTBBThreads() const
	{
	tbb::global_control gc{ tbb::global_control::parameter::max_allowed_parallelism,
	m_config.nThreads };

	TestHistos null_hist;

	return tbb::parallel_reduce(
	tbb::blocked_range<std::size_t>( 0, m_config.nEvents ), null_hist,
	[this]( const tbb::blocked_range<std::size_t>& range,
	const TestHistos& init ) -> TestHistos
	{
	std::cout << "Thread "
	<< tbb::this_task_arena::current_thread_index()
	<< " will generate " << range.size()
	<< " events: " << range.begin() << " - " << range.end()
	<< std::endl;
	TestHistos tmp{ init };
	tmp.add( runDecayBody( range.begin(), range.size() ) );
	return tmp;
	},
	[]( const TestHistos& lhs, const TestHistos& rhs ) -> TestHistos
	{
	TestHistos tmp{ lhs };
	tmp.add( rhs );
	return tmp;
	} );
	}

	TestHistos TestDecayModel::runDecayBody( const std::size_t firstEvent,
	const std::size_t nEvents ) const
	{
	// Initialise the EvtGen object and hence the EvtPDL tables
	// The EvtGen object is used by generateEvents, while the
	// latter are also used within createDecFile

	// Define the random number generator
	static thread_local auto randomEngine{ std::make_unique<EvtMTRandomEngine>() };

	// TODO - need to streamline the extra models stuff
	static thread_local EvtAbsRadCorr* radCorrEngine{ nullptr };
	std::list<EvtDecayBase*> extraModels;

	static thread_local bool initialised{ false };

	#ifdef EVTGEN_EXTERNAL
	if ( !initialised ) {
	bool convertPythiaCodes( false );
	bool useEvtGenRandom( true );
	bool seedTauolaFortran( true );
	EvtExternalGenList genList( convertPythiaCodes, "", "gamma",
	useEvtGenRandom, seedTauolaFortran );
	switch ( m_config.fsrGenerator ) {
	case FSRGenerator::PHOTOS:
	radCorrEngine = genList.getPhotosModel();
	break;
	case FSRGenerator::SherpaPhotons1:
	radCorrEngine = genList.getSherpaPhotonsModel( 1e-7, 1, 0 );
	break;
	case FSRGenerator::SherpaPhotons20:
	radCorrEngine = genList.getSherpaPhotonsModel( 1e-7, 2, 0 );
	break;
	case FSRGenerator::SherpaPhotons21:
	radCorrEngine = genList.getSherpaPhotonsModel( 1e-7, 2, 1 );
	break;
	}
	extraModels = genList.getListOfModels();
	}
	#endif

	static thread_local EvtGen theGen{ "../DECAY.DEC", "../evt.pdl",
	randomEngine.get(), radCorrEngine,
	&extraModels };

	// Creates a decay file based on json file input
	// We don't want this to be called by every thread!
	std::call_once( m_createDecFile_threadlock, [this]() { createDecFile(); } );

	// Read the decay file
	if ( !initialised ) {
	theGen.readUDecay( m_config.decFileName.c_str() );
	initialised = true;
	}

	// Define the histograms to be saved
	TestHistos theHistos{ m_config.testHistograms };

	// Generate events and fill histograms
	generateEvents( theGen, theHistos, firstEvent, nEvents );

	return theHistos;
	}

	void TestDecayModel::createDecFile() const
	{
	// Create (or overwrite) the decay file
	std::ofstream decFile{ m_config.decFileName };

	// Create daughter aliases if needed
	std::vector<std::string> aliasPrefix;
	for ( std::size_t daughter_index{ 0 };
	daughter_index < m_config.daughterNames.size(); daughter_index++ ) {
	if ( !m_config.grandDaughterNames.empty() &&
	!m_config.grandDaughterNames[daughter_index].empty() ) {
	decFile << "Alias My" << m_config.daughterNames[daughter_index] << " "
	<< m_config.daughterNames[daughter_index] << std::endl;
	if ( m_config.doConjDecay[daughter_index + 1] ) {
	const EvtId daugID{
	EvtPDL::getId( m_config.daughterNames[daughter_index] ) };
	const EvtId daugConjID{ EvtPDL::chargeConj( daugID ) };
	const std::string conjName{ daugConjID.getName() };
	std::string conjName_Alias{ daugConjID.getName() };
	if ( std::find( std::begin( m_config.daughterNames ),
	std::end( m_config.daughterNames ),
	daugConjID.getName() ) !=
	std::end( m_config.daughterNames ) ) {
	conjName_Alias = conjName_Alias + "_" + daughter_index;
	}
	decFile << "Alias My" << conjName_Alias << " " << conjName
	<< std::endl;
	decFile << "ChargeConj My"
	<< m_config.daughterNames[daughter_index] << " My"
	<< conjName_Alias << std::endl;
	} else if ( m_config.doConjDecay[0] ) {
	decFile << "ChargeConj My"
	<< m_config.daughterNames[daughter_index] << " My"
	<< m_config.daughterNames[daughter_index] << std::endl;
	}
	aliasPrefix.push_back( "My" );
	} else {
	aliasPrefix.push_back( "" );
	}
	}

	for ( const auto& iExtra : m_config.extraCommands ) {
	decFile << iExtra << std::endl;
	}

	// Parent decay
	decFile << "Decay " << m_config.parentName << std::endl;
	decFile << "1.0";

	for ( std::size_t daughter_index{ 0 };
	daughter_index < m_config.daughterNames.size(); daughter_index++ ) {
	decFile << " " << aliasPrefix[daughter_index]
	<< m_config.daughterNames[daughter_index];
	}

	decFile << " " << m_config.modelNames[0];

	for ( const auto& par : m_config.modelParameters[0] ) {
	decFile << " " << par;
	}

	decFile << ";" << std::endl;
	decFile << "Enddecay" << std::endl;
	if ( m_config.doConjDecay[0] ) {
	EvtId parID{ EvtPDL::getId( m_config.parentName ) };
	EvtId parConjID{ EvtPDL::chargeConj( parID ) };
	decFile << "CDecay " << parConjID.getName() << std::endl;
	}

	// Daughter decays into granddaughters
	for ( std::size_t daughter_index{ 0 };
	daughter_index < m_config.grandDaughterNames.size(); daughter_index++ ) {
	if ( m_config.grandDaughterNames[daughter_index].empty() )
	continue;
	decFile << "Decay " << aliasPrefix[daughter_index]
	<< m_config.daughterNames[daughter_index] << std::endl;
	decFile << "1.0";
	for ( std::size_t grandDaughter_index{ 0 };
	grandDaughter_index <
	m_config.grandDaughterNames[daughter_index].size();
	grandDaughter_index++ ) {
	decFile << " "
	<< m_config.grandDaughterNames[daughter_index][grandDaughter_index];
	}
	decFile << " " << m_config.modelNames[daughter_index + 1];
	for ( const auto& par : m_config.modelParameters[daughter_index + 1] ) {
	decFile << " " << par;
	}
	decFile << ";" << std::endl;
	decFile << "Enddecay" << std::endl;
	if ( m_config.doConjDecay[daughter_index + 1] ) {
	EvtId daugID{
	EvtPDL::getId( m_config.daughterNames[daughter_index] ) };
	EvtId daugConjID{ EvtPDL::chargeConj( daugID ) };
	std::string conjName_Alias{ daugConjID.getName() };
	if ( std::find( std::begin( m_config.daughterNames ),
	std::end( m_config.daughterNames ),
	daugConjID.getName() ) !=
	std::end( m_config.daughterNames ) ) {
	conjName_Alias = conjName_Alias + "_" + daughter_index;
	}
	decFile << "CDecay " << aliasPrefix[daughter_index]
	<< conjName_Alias << std::endl;
	}
	}

	decFile << "End" << std::endl;

	decFile.close();
	}

	void TestDecayModel::generateEvents( EvtGen& theGen, TestHistos& theHistos,
	const std::size_t firstEvent,
	const std::size_t nEvents ) const
	{
	// Generate the decays
	const EvtId parId{ EvtPDL::getId( m_config.parentName.c_str() ) };
	const EvtId conjId{ m_config.doConjDecay[0] ? EvtPDL::chargeConj( parId )
	: parId };

	for ( std::size_t i{ firstEvent }; i < ( firstEvent + nEvents ); i++ ) {
	if ( i % 1000 == 0 ) {
	std::cout << "Event " << firstEvent + nEvents - i << std::endl;
	}

	// seed the RNG based on the event number
	EvtRandom::setSeed( m_config.rngSeed + i );

	// Initial 4-momentum and particle
	EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 );
	EvtParticle* parent{ nullptr };
	if ( EvtRandom::Flat() < 0.5 ) {
	parent = EvtParticleFactory::particleFactory( parId, pInit );
	} else {
	parent = EvtParticleFactory::particleFactory( conjId, pInit );
	}

	theGen.generateDecay( parent );

	// Check for mixing (and fill histogram)
	EvtParticle* prodParent{ nullptr };
	TH1* mixedHist{ theHistos.getMixedHist() };
	if ( mixedHist ) {
	if ( parent->getNDaug() == 1 ) {
	prodParent = parent;
	parent = prodParent->getDaug( 0 );

	mixedHist->Fill( 1 );
	} else {
	mixedHist->Fill( 0 );
	}
	}

	// To debug

	if ( m_config.debugFlag ) {
	std::cout << "Parent PDG code: " << parent->getPDGId()
	<< " has daughters " << parent->getNDaug() << std::endl;
	for ( std::size_t iDaughter{ 0 }; iDaughter < parent->getNDaug();
	iDaughter++ ) {
	std::cout << "Parent PDG code of daughter " << iDaughter
	<< " : " << parent->getDaug( iDaughter )->getPDGId()
	<< " has daughters "
	<< parent->getDaug( iDaughter )->getNDaug()
	<< std::endl;

	for ( std::size_t iGrandDaughter{ 0 };
	iGrandDaughter < parent->getDaug( iDaughter )->getNDaug();
	iGrandDaughter++ ) {
	std::cout << "Parent PDG code of grand daughter "
	<< iGrandDaughter << " : "
	<< parent->getDaug( iDaughter )
	->getDaug( iGrandDaughter )
	->getPDGId()
	<< " has daughters "
	<< parent->getDaug( iDaughter )
	->getDaug( iGrandDaughter )
	->getNDaug()
	<< std::endl;
	}
	}
	}

	int leadingChargedDaughter = -1;
	const std::string fsrStr = "_FSRPhotons";
	const std::string ewStr = "_EnergyWeight";

	const std::string perDaughter = "_perDaughter";

	// Store information
	for ( auto& [info, hist] : theHistos.get1DHistos() ) {
	if ( !hist ) {
	continue;
	}

	const std::string varName = info.getName();
	std::string reducedVarName = varName;

	const std::string::size_type findFSRstr = varName.find( fsrStr );

	const std::string::size_type findPerDaugStr = varName.find(
	perDaughter );

	// If the variable name does not have the substrings '_perDaughter' and '_FSRPhotons', then add an entry per event and continue
	if ( findFSRstr == std::string::npos &&
	findPerDaugStr == std::string::npos ) {
	const double value{
	getValue( parent, varName, info.getd1(), info.getd2() ) };

	hist->Fill( value );
	continue;
	}
	// If the variable name has the substring '_perDaughter' and does not have the substring '_FSRPhotons', then add an entry per daughter and continue
	else if ( findFSRstr == std::string::npos ) {
	reducedVarName.erase( findPerDaugStr, perDaughter.length() );

	// Figure out whether the variable indicates to be saved for a requested particle type
	const std::string::size_type findIsStr = reducedVarName.find(
	"_is" );
	const std::string requestedType = findIsStr == std::string::npos
	? ""
	: reducedVarName.substr(
	findIsStr + 3 );

	if ( !requestedType.empty() ) {
	reducedVarName.erase( findIsStr,
	reducedVarName.length() - findIsStr );
	}

	for ( int iDaug{ 0 }; iDaug < (int)parent->getNDaug(); iDaug++ ) {
	const double value{
	getValue( parent, reducedVarName, iDaug + 1, 0 ) };

	const int partGroup = getPartGroup(
	parent->getDaug( iDaug )->getPDGId() );

	if ( requestedType.empty() \|\|
	isPartType( partGroup, requestedType ) ) {
	hist->Fill( value );
	}
	}
	continue;
	}

	// If otherwise the variable contains the substring '_FSRPhotons', then add an entry per photon
	if ( leadingChargedDaughter == -1 ) {
	leadingChargedDaughter = findChargedDaughterWithMaxE( parent );
	}

	reducedVarName.erase( findFSRstr, fsrStr.length() );

	const std::string::size_type findEwStr = reducedVarName.find( ewStr );

	// If variable name has substring '_EnergyWeight' then add an entry weighted by the photon energy
	bool energyWeight = false;
	if ( findEwStr != std::string::npos ) {
	reducedVarName.erase( findEwStr, ewStr.length() );
	energyWeight = true;
	}

	for ( int iDaug{ 0 }; iDaug < (int)parent->getNDaug(); iDaug++ ) {
	const EvtParticle* iDaughter = parent->getDaug( iDaug );

	if ( iDaughter->getAttribute( "FSR" ) == 1 ) {
	const double value{ getValue( parent, reducedVarName,
	iDaug + 1,
	leadingChargedDaughter + 1 ) };
	if ( energyWeight ) {
	const double photonEnergy = iDaughter->getP4().get( 0 );
	hist->Fill( value, photonEnergy );
	} else {
	hist->Fill( value );
	}
	}
	}
	}

	for ( auto& [info, hist] : theHistos.get2DHistos() ) {
	if ( !hist ) {
	continue;
	}

	const double valueX{ getValue( parent,
	info.getName( HistInfo::Axis::X ),
	info.getd1( HistInfo::Axis::X ),
	info.getd2( HistInfo::Axis::X ) ) };
	const double valueY{ getValue( parent,
	info.getName( HistInfo::Axis::Y ),
	info.getd1( HistInfo::Axis::Y ),
	info.getd2( HistInfo::Axis::Y ) ) };
	hist->Fill( valueX, valueY );
	}

	if ( m_config.debugFlag ) {
	if ( prodParent ) {
	prodParent->printTree();
	} else {
	parent->printTree();
	}
	}

	// Cleanup
	if ( prodParent ) {
	prodParent->deleteTree();
	} else {
	parent->deleteTree();
	}
	}
	}

	int TestDecayModel::findChargedDaughterWithMaxE( const EvtParticle* parent ) const
	{
	/* This function returns the index of the charged daughter with the highest energy
	* following the sign convention below. */

	double max_E = 0;
	double max_index = 0;

	const int parentCh3 = EvtPDL::chg3( parent->getId() );

	for ( int iDaug{ 0 }; iDaug < (int)parent->getNDaug(); iDaug++ ) {
	const EvtParticle* iDaughter = parent->getDaug( iDaug );
	const int dauCh3 = EvtPDL::chg3( iDaughter->getId() );
	const double daugE = iDaughter->getP4LabBeforeFSR().get( 0 );
	// Sign convention: take negative daughter if mother is neutral,
	// otherwise the one with the same sign as the mother.
	if ( ( ( parentCh3 == 0 && dauCh3 < 0 ) \|\| ( parentCh3 * dauCh3 > 0 ) ) &&
	daugE > max_E ) {
	max_E = daugE;
	max_index = iDaug;
	}
	}
	return max_index;
	}

	double TestDecayModel::getValue( const EvtParticle* parent,
	const std::string& varName, const int d1,
	const int d2 ) const
	{
	double value{ std::numeric_limits<double>::quiet_NaN() };
	if ( !parent ) {
	return value;
	}

	const int NDaugMax( parent->getNDaug() );
	// If variable name contains "_daugX", we are interested in daughters of daughter X
	// Else we are interested in daughters
	const EvtParticle* selectedParent{ parent };
	std::string selectedVarName{ varName };
	if ( varName.find( "_daug" ) != std::string::npos ) {
	// Get daughter index from last character in string
	const int iDaughter{ varName.back() - '0' };
	selectedVarName = varName.substr( 0, varName.size() - 6 );
	if ( iDaughter > 0 && iDaughter <= NDaugMax ) {
	selectedParent = parent->getDaug( iDaughter - 1 );
	} else {
	return value;
	}
	}

	const int sel_NDaugMax( selectedParent->getNDaug() );
	const EvtParticle* par1{ d1 > 0 && d1 <= sel_NDaugMax
	? selectedParent->getDaug( d1 - 1 )
	: nullptr };
	const EvtParticle* par2{ d2 > 0 && d2 <= sel_NDaugMax
	? selectedParent->getDaug( d2 - 1 )
	: nullptr };
	const EvtParticle* par3{ sel_NDaugMax > 0
	? selectedParent->getDaug( sel_NDaugMax - 1 )
	: nullptr };

	// 4-momenta in parent rest frame
	const EvtVector4R p1{ par1 != nullptr ? par1->getP4() : EvtVector4R() };
	const EvtVector4R p2{ par2 != nullptr ? par2->getP4() : EvtVector4R() };
	const EvtVector4R p3{ par3 != nullptr ? par3->getP4() : EvtVector4R() };

	// 4-momenta in lab frame (1st parent in decay tree)
	const EvtVector4R p1_lab{ par1 != nullptr ? par1->getP4Lab() : EvtVector4R() };
	const EvtVector4R p2_lab{ par2 != nullptr ? par2->getP4Lab() : EvtVector4R() };
	const EvtVector4R p3_lab{ par3 != nullptr ? par3->getP4Lab() : EvtVector4R() };

	if ( !selectedVarName.compare( "id" ) ) {
	// StdHep ID of one of the daughters (controlled by d1) or the parent
	if ( par1 ) {
	value = par1->getPDGId();
	} else {
	value = selectedParent->getPDGId();
	}
	} else if ( !selectedVarName.compare( "particleType" ) ) {
	if ( par1 ) {
	value = getPartGroup( par1->getPDGId() );
	} else {
	value = getPartGroup( selectedParent->getPDGId() );
	}
	} else if ( !selectedVarName.compare( "nDaug" ) ) {
	// Number of daughters
	value = sel_NDaugMax;

	} else if ( !selectedVarName.compare( "parMass" ) ) {
	// Parent invariant mass
	value = selectedParent->mass();

	} else if ( !selectedVarName.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();
	}

	} else if ( !selectedVarName.compare( "massSq" ) ) {
	// Invariant mass
	if ( d2 != 0 ) {
	// Invariant 4-mass combination of particles d1 and d2
	value = ( p1 + p2 ).mass2();
	} else {
	// Invariant mass of particle d1 only
	value = p1.mass2();
	}
	} else if ( !selectedVarName.compare( "mPrime" ) ) {
	// pick up first unused particle rather than last one
	if ( sel_NDaugMax != 3 ) {
	return -1;
	}
	int unused{ 0 };
	for ( int ii{ 1 }; ii <= sel_NDaugMax; ++ii ) {
	if ( ii != d1 && ii != d2 ) {
	unused = ii;
	break;
	}
	}
	if ( unused == 0 ) {
	unused = sel_NDaugMax;
	}
	const auto parL{ selectedParent->getDaug( unused - 1 ) };
	// const auto& pL = parL->getP4();

	const double mB{ selectedParent->mass() };
	const double m1{ par1->mass() };
	const double m2{ par2->mass() };
	const double m3{ parL->mass() };
	const double m12{ ( p1 + p2 ).mass() };
	const double m12norm{
	2 * ( ( m12 - ( m1 + m2 ) ) / ( mB - ( m1 + m2 + m3 ) ) ) - 1 };
	value = acos( m12norm ) / EvtConst::pi;

	} else if ( !selectedVarName.compare( "thetaPrime" ) ) {
	// pick up first unused particle rather than last one
	if ( sel_NDaugMax != 3 ) {
	return -1;
	}
	int unused{ 0 };
	for ( int ii{ 1 }; ii <= sel_NDaugMax; ++ii ) {
	if ( ii != d1 && ii != d2 ) {
	unused = ii;
	break;
	}
	}
	if ( unused == 0 ) {
	unused = sel_NDaugMax;
	}
	const auto parL{ selectedParent->getDaug( unused - 1 ) };
	const auto& pL{ parL->getP4() };

	const double mB{ selectedParent->mass() };
	const double m1{ p1.mass() };
	const double m2{ p2.mass() };
	const double m3{ pL.mass() };
	double mBSq{ mB * mB };
	double m1Sq{ m1 * m1 };
	double m2Sq{ m2 * m2 };
	double m3Sq{ m3 * m3 };
	const double m12{ ( p1 + p2 ).mass() };
	const double m13{ ( p1 + p3 ).mass() };
	const double m12Sq{ m12 * m12 };
	const double m13Sq{ m13 * m13 };
	double en1{ ( m12Sq - m2Sq + m1Sq ) / ( 2.0 * m12 ) };
	double en3{ ( mBSq - m12Sq - m3Sq ) / ( 2.0 * m12 ) };
	double p1_12{ std::sqrt( en1 * en1 - m1Sq ) };
	double p3_12{ std::sqrt( en3 * en3 - m3Sq ) };
	double cosTheta{ ( -m13Sq + m1Sq + m3Sq + 2. * en1 * en3 ) /
	( 2. * p1_12 * p3_12 ) };
	value = acos( cosTheta ) / EvtConst::pi;

	} else if ( !selectedVarName.compare( "pSumSq" ) ) {
	// Invariant momentum sum squared
	value = ( p1 + p2 ).mass2();

	} else if ( !selectedVarName.compare( "pDiffSq" ) ) {
	// Invariant momentum difference squared
	value = ( p1 - p2 ).mass2();

	} else if ( !selectedVarName.compare( "mass3" ) ) {
	// Invariant mass of 3 daughters
	value = ( p1 + p2 + p3 ).mass();

	} else if ( !selectedVarName.compare( "mass3_specified" ) ) {
	// Invariant mass of 3 daughters, d1 is first daughter
	// second daughter is d1 + 1, d2 is last daughter
	const EvtParticle* daug2{ selectedParent->getDaug( d1 ) };
	const EvtParticle* daug3{ selectedParent->getDaug( d2 - 1 ) };

	const EvtVector4R p2_specified{ daug2 != nullptr ? daug2->getP4Lab()
	: EvtVector4R() };
	const EvtVector4R p3_specified{ daug3 != nullptr ? daug3->getP4Lab()
	: EvtVector4R() };

	value = ( p1 + p2_specified + p3_specified ).mass();

	} else if ( !selectedVarName.compare( "cosTheta3" ) ) {
	// Cosine of the polar angle of the momentum of d1 + d2 + d3
	const EvtVector4R p123{ p1 + p2 + p3 };
	if ( p123.d3mag() > 0.0 ) {
	value = p123.get( 3 ) / p123.d3mag();
	}

	} else if ( !selectedVarName.compare( "pLab" ) ) {
	// Momentum of particle d1 in lab frame
	value = p1_lab.d3mag();

	} else if ( !selectedVarName.compare( "p" ) ) {
	// Momentum of particle d1 in parent rest frame
	value = p1.d3mag();

	} else if ( !selectedVarName.compare( "pLabSq" ) ) {
	// Momentum squared of particle d1 (in lab frame)
	const double p1_lab_x{ p1_lab.get( 1 ) };
	const double p1_lab_y{ p1_lab.get( 2 ) };
	const double p1_lab_z{ p1_lab.get( 3 ) };
	value = p1_lab_x * p1_lab_x + p1_lab_y * p1_lab_y + p1_lab_z * p1_lab_z;

	} else if ( !selectedVarName.compare( "pSq" ) ) {
	// Momentum squared of particle d1 (in lab frame)
	const double p1_x{ p1.get( 1 ) };
	const double p1_y{ p1.get( 2 ) };
	const double p1_z{ p1.get( 3 ) };
	value = p1_x * p1_x + p1_y * p1_y + p1_z * p1_z;

	} else if ( !selectedVarName.compare( "pxLab" ) ) {
	// x momentum of particle d1 in lab frame
	value = p1_lab.get( 1 );

	} else if ( !selectedVarName.compare( "px" ) ) {
	// x momentum of particle d1 in parent frame
	value = p1.get( 1 );

	} else if ( !selectedVarName.compare( "pyLab" ) ) {
	// y momentum of particle d1 in lab frame
	value = p1_lab.get( 2 );

	} else if ( !selectedVarName.compare( "py" ) ) {
	// y momentum of particle d1 in parent frame
	value = p1.get( 2 );

	} else if ( !selectedVarName.compare( "pzLab" ) ) {
	// z momentum of particle d1 in lab frame
	value = p1_lab.get( 3 );

	} else if ( !selectedVarName.compare( "pz" ) ) {
	// z momentum of particle d1 in parent frame
	value = p1.get( 3 );

	} else if ( !selectedVarName.compare( "cosHel" ) \|\|
	!selectedVarName.compare( "absCosHel" ) \|\|
	!selectedVarName.compare( "cosHelParent" ) ) {
	// Cosine of helicity angle
	EvtVector4R p12;
	EvtVector4R p1Res;

	if ( !selectedVarName.compare( "cosHelParent" ) ) {
	p12 = selectedParent->getP4Lab();
	const EvtVector4R boost{ p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ),
	-p12.get( 3 ) };
	// Momentum of particle d1 in resonance frame, p1Res
	p1Res = boostTo( p1_lab, boost );
	} else if ( d2 != 0 ) {
	// Resonance center-of-mass system (d1 and d2)
	p12 = p1_lab + p2_lab;
	// 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
	p1Res = boostTo( p1_lab, boost );
	} else {
	// The resonance is d1
	p12 = p1;

	// Find its first daughter
	const EvtParticle* gpar{ par1 != nullptr ? par1->getDaug( 0 )
	: nullptr };

	p1Res = gpar != nullptr ? gpar->getP4() : EvtVector4R();
	}

	// Cosine of angle between p1Res and momentum of resonance in parent frame
	const double p1ResMag{ p1Res.d3mag() };
	const double p12Mag{ p12.d3mag() };
	if ( p1ResMag > 0.0 && p12Mag > 0.0 ) {
	value = -p1Res.dot( p12 ) / ( p1ResMag * p12Mag );
	}

	if ( !selectedVarName.compare( "absCosHel" ) ) {
	value = std::abs( value );
	}

	} else if ( !selectedVarName.compare( "cosHelTau" ) ) {
	// Works only for B -> X nu_1 tau -> pi nu_2.
	// Cosine of helicity angle between pi momentum and opposite W momentum -(nu_1 + tau)
	// in the tau rest frame. Index d1 must match with tau.
	// p3 (momentum of last daughter) is taken as the neutrino momentum

	// W momentum
	const EvtVector4R p_W{ p1_lab + p3_lab };

	// Index d2 must match the index of the pion (daughter of tau)
	const EvtParticle* pion{
	selectedParent->getDaug( d1 - 1 )->getDaug( d2 - 1 ) };
	const EvtVector4R p_pion{ pion != nullptr ? pion->getP4Lab()
	: EvtVector4R() };

	// Boost vector to tau frame
	const EvtVector4R boost{ p1_lab.get( 0 ), -p1_lab.get( 1 ),
	-p1_lab.get( 2 ), -p1_lab.get( 3 ) };

	// Boost both momenta to tau frame
	const EvtVector4R p_W_boosted{ boostTo( p_W, boost ) };
	const EvtVector4R p_pion_boosted{ boostTo( p_pion, boost ) };
	// Cosine of angle between opposite W momentum and pion momentum in tau frame
	const double p_W_boostedMag{ p_W_boosted.d3mag() };
	const double p_pion_boostedMag{ p_pion_boosted.d3mag() };

	if ( p_W_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
	value = -p_W_boosted.dot( p_pion_boosted ) /
	( p_W_boostedMag * p_pion_boostedMag );
	}

	} else if ( !selectedVarName.compare( "cosHelDiTau" ) \|\|
	!selectedVarName.compare( "cosHelDiTau_over05" ) \|\|
	!selectedVarName.compare( "cosHelDiTau_under05" ) ) {
	// Works for B -> tau (pi nu) tau -> (pi nu), B -> phi (-> KK) l l decays,
	// B -> 4 pions, or similar decays..
	// Cosine of helicity angle between pi momentum and opposite B momentum in tau rest frame.
	// Index d1 must match with tau

	if ( sel_NDaugMax < 2 \|\| sel_NDaugMax > 4 ) {
	return value;
	}

	// B momentum
	const EvtVector4R p_B{ selectedParent->getP4Lab() };

	// Index d2 must match the index of the pion (daughter of tau)
	const EvtParticle* pion_1{
	sel_NDaugMax <= 3
	? selectedParent->getDaug( d1 - 1 )->getDaug( d2 - 1 )
	: selectedParent->getDaug( d1 - 1 ) };
	const EvtVector4R p_pion_1{ pion_1 != nullptr ? pion_1->getP4Lab()
	: EvtVector4R() };

	const EvtVector4R p_first_daughter{
	sel_NDaugMax <= 3 ? selectedParent->getDaug( d1 - 1 )->getP4Lab()
	: selectedParent->getDaug( d1 - 1 )->getP4Lab() +
	selectedParent->getDaug( d1 )->getP4Lab() };

	// Boost vector to tau frame
	const EvtVector4R boost_1{ p_first_daughter.get( 0 ),
	-p_first_daughter.get( 1 ),
	-p_first_daughter.get( 2 ),
	-p_first_daughter.get( 3 ) };

	// Boost both momenta to tau frame
	const EvtVector4R p_B_boosted_1{ boostTo( p_B, boost_1 ) };
	const EvtVector4R p_pion_boosted_1{ boostTo( p_pion_1, boost_1 ) };
	// Cosine of angle between opposite W momentum and pion momentum in tau frame
	const double p_B_boosted_1_Mag{ p_B_boosted_1.d3mag() };
	const double p_pion_boosted_1_Mag{ p_pion_boosted_1.d3mag() };

	double hel1{ std::numeric_limits<double>::quiet_NaN() };
	double hel{ std::numeric_limits<double>::quiet_NaN() };

	if ( p_B_boosted_1_Mag > 0.0 && p_pion_boosted_1_Mag > 0.0 ) {
	hel1 = -p_B_boosted_1.dot( p_pion_boosted_1 ) /
	( p_B_boosted_1_Mag * p_pion_boosted_1_Mag );
	}

	if ( ( !selectedVarName.compare( "cosHelDiTau_over05" ) ) \|\|
	( !selectedVarName.compare( "cosHelDiTau_under05" ) ) ) {
	// Works for B -> tau (pi nu) tau -> (pi nu) or similar decays.
	// Cosine of helicity angle between pi momentum and opposite B momentum in tau rest frame
	// Index d1 must match with tau; cosHelicity above +0.5 or below -0.5

	// Index d2 must match the index of the pion (daughter of tau)
	const EvtParticle* pion{
	sel_NDaugMax <= 3 ? selectedParent->getDaug( 0 )->getDaug( d2 - 1 )
	: selectedParent->getDaug( 0 ) };
	const EvtVector4R p_pion{ pion != nullptr ? pion->getP4Lab()
	: EvtVector4R() };
	// Boost vector to tau frame
	const EvtVector4R p_second_daughter{
	sel_NDaugMax == 2
	? selectedParent->getDaug( 0 )->getP4Lab()
	: selectedParent->getDaug( 0 )->getP4Lab() +
	selectedParent->getDaug( 1 )->getP4Lab() };

	const EvtVector4R boost{ p_second_daughter.get( 0 ),
	-p_second_daughter.get( 1 ),
	-p_second_daughter.get( 2 ),
	-p_second_daughter.get( 3 ) };

	// Boost both momenta to tau frame
	const EvtVector4R p_B_boosted{ boostTo( p_B, boost ) };
	const EvtVector4R p_pion_boosted{ boostTo( p_pion, boost ) };
	// Cosine of angle between opposite W momentum and pion momentum in tau frame
	const double p_B_boostedMag{ p_B_boosted.d3mag() };
	const double p_pion_boostedMag{ p_pion_boosted.d3mag() };

	if ( p_B_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
	hel = -p_B_boosted.dot( p_pion_boosted ) /
	( p_B_boostedMag * p_pion_boostedMag );
	}
	}

	if ( !selectedVarName.compare( "cosHelDiTau" ) \|\|
	( hel > 0.5 && !selectedVarName.compare( "cosHelDiTau_over05" ) ) \|\|
	( hel < -0.5 && !selectedVarName.compare( "cosHelDiTau_under05" ) ) ) {
	value = hel1;
	}

	} else if ( !selectedVarName.compare( "cosAcoplanarityAngle" ) \|\|
	!selectedVarName.compare( "acoplanarityAngle" ) ) {
	// Acoplanarity angle or cosine for B -> tau (pi nu) tau -> (pi nu),
	// B -> phi (-> KK) l l decays, B -> 4 pions, or similar decays

	value = getCosAcoplanarityAngle( selectedParent, sel_NDaugMax, d1, d2 );
	if ( !selectedVarName.compare( "acoplanarityAngle" ) ) {
	value = std::acos( value ) * 180.0 / EvtConst::pi; // degrees
	}

	} else if ( !selectedVarName.compare( "cosTheta" ) ) {
	// Cosine of polar angle of first daughter in lab frame
	const double p1_lab_mag{ p1_lab.d3mag() };
	if ( p1_lab_mag > 0.0 ) {
	value = p1_lab.get( 3 ) / p1_lab_mag;
	}

	} else if ( !selectedVarName.compare( "phi" ) ) {
	// Azimuthal angle of first daughter in lab frame (degrees)
	const double p1_lab_mag{ p1_lab.d3mag() };
	if ( p1_lab_mag > 0.0 ) {
	value = atan2( p1_lab.get( 1 ), p1_lab.get( 2 ) ) * 180.0 /
	EvtConst::pi;
	}

	} else if ( !selectedVarName.compare( "openingAngle" ) ) {
	// Polar angle between first and second daughters in parent's frame

	const double cost{ p1.dot( p2 ) / ( p1.d3mag() * p2.d3mag() ) };

	value = acos( cost ) * 180.0 / EvtConst::pi;

	} else if ( !selectedVarName.compare( "decayangle" ) ) {
	// Polar angle between first and second daughters in lab frame

	const EvtVector4R p{ selectedParent->getP4() };
	const EvtVector4R q{ p1 + p2 };
	const EvtVector4R d{ p1 };

	const double cost{ EvtDecayAngle( p, q, d ) };

	value = acos( cost ) * 180.0 / EvtConst::pi;

	} else if ( !selectedVarName.compare( "decayangle3" ) ) {
	// Polar angle between combined first and second daughters
	// with the third daughter in lab frame

	const EvtVector4R p{ selectedParent->getP4() };
	const EvtVector4R q{ p1 + p2 + p3 };
	const EvtVector4R d{ p1 + p2 };

	const double cost{ EvtDecayAngle( p, q, d ) };

	value = acos( cost ) * 180.0 / EvtConst::pi;

	} else if ( !selectedVarName.compare( "decayangle_BTO4PI" ) ) {
	// Polar angle between first and second daughters in lab frame.
	// Used in PIPIPI (BTO4PI_CP) model

	const EvtVector4R p{ p1 + p2 + p3 };
	const EvtVector4R q{ p1 + p2 };
	const EvtVector4R d{ p1 };

	const double cost{ EvtDecayAngle( p, q, d ) };

	value = acos( cost ) * 180.0 / EvtConst::pi;

	} else if ( !selectedVarName.compare( "chi" ) ) {
	// Chi angle (degrees) using all 4 daughters and parent 4 mom

	const EvtParticle* daug1{ selectedParent->getDaug( d1 - 1 ) };
	const EvtParticle* daug2{ selectedParent->getDaug( d1 ) };
	const EvtParticle* daug3{ selectedParent->getDaug( d1 + 1 ) };
	const EvtParticle* daug4{ selectedParent->getDaug( d1 + 2 ) };

	const EvtVector4R p_parent{ selectedParent->getP4() };
	const EvtVector4R p_daug1{ daug1 != nullptr ? daug1->getP4()
	: EvtVector4R() };
	const EvtVector4R p_daug2{ daug2 != nullptr ? daug2->getP4()
	: EvtVector4R() };
	const EvtVector4R p_daug3{ daug3 != nullptr ? daug3->getP4()
	: EvtVector4R() };
	const EvtVector4R p_daug4{ daug4 != nullptr ? daug4->getP4()
	: EvtVector4R() };

	const double chi{
	EvtDecayAngleChi( p_parent, p_daug1, p_daug2, p_daug3, p_daug4 ) };
	value = chi * 180.0 / EvtConst::pi;

	} else if ( !selectedVarName.compare( "cosThetaResNorm" ) ) {
	// P -> R p4 -> (p1 p2 p3) p4, where the resonance R decays to p1 p2 p3.
	// Theta is the angle between the normal of the plane containing p1, p2 & p3
	// and the bachelor particle p4 in the rest frame of resonance R.
	// The normal vector is given by the cross product p3 x p1

	if ( sel_NDaugMax > 1 ) {
	const EvtParticle* res{ selectedParent->getDaug( 0 ) };
	const EvtParticle* bac{ selectedParent->getDaug( 1 ) };

	// Check resonance has 3 daughters
	if ( res != nullptr && res->getNDaug() == 3 ) {
	const EvtParticle* daug1 = res->getDaug( 0 );
	const EvtParticle* daug3 = res->getDaug( 2 );

	// 4-momenta in base parent P lab frame
	const EvtVector4R p4_Res{ res->getP4Lab() };
	const EvtVector4R p4_p4{ bac != nullptr ? bac->getP4Lab()
	: EvtVector4R() };
	const EvtVector4R p4_p1{ daug1 != nullptr ? daug1->getP4Lab()
	: EvtVector4R() };
	const EvtVector4R p4_p3{ daug3 != nullptr ? daug3->getP4Lab()
	: EvtVector4R() };

	// Boost 4-vector for resonance frame
	const EvtVector4R boost{ p4_Res.get( 0 ), -p4_Res.get( 1 ),
	-p4_Res.get( 2 ), -p4_Res.get( 3 ) };

	// Momentum of p1 and p3 in resonance frame
	const EvtVector4R p1Res{ boostTo( p4_p1, boost ) };
	const EvtVector4R p3Res{ boostTo( p4_p3, boost ) };

	// Plane normal vector (just uses 3-momentum components)
	const EvtVector4R norm{ p3Res.cross( p1Res ) };

	// Momentum of p4 in resonance frame
	const EvtVector4R p4Res{ boostTo( p4_p4, boost ) };

	// Cosine of the angle between the normal and p4 in the resonance frame
	const double normMag{ norm.d3mag() };
	const double p4ResMag{ p4Res.d3mag() };
	if ( normMag > 0.0 && p4ResMag > 0.0 ) {
	value = norm.dot( p4Res ) / ( normMag * p4ResMag );
	}
	}
	}

	} else if ( !selectedVarName.compare( "cosBetaRes" ) ) {
	// For resonance P (parent) -> p1 p2 p3, beta is the
	// angle between p1 & p3 in the (p1 + p2) rest frame
	if ( sel_NDaugMax > 2 ) {
	const EvtParticle* daug1 = selectedParent->getDaug( 0 );
	const EvtParticle* daug2 = selectedParent->getDaug( 1 );
	const EvtParticle* daug3 = selectedParent->getDaug( 2 );

	// 4-momenta in base parent frame
	const EvtVector4R p4_p1{ daug1 != nullptr ? daug1->getP4Lab()
	: EvtVector4R() };
	const EvtVector4R p4_p2{ daug2 != nullptr ? daug2->getP4Lab()
	: EvtVector4R() };
	const EvtVector4R p4_p3{ daug3 != nullptr ? daug3->getP4Lab()
	: EvtVector4R() };

	// p1 + p2
	const EvtVector4R p12{ p4_p1 + p4_p2 };

	// Boost 4-vector for p12 frame
	const EvtVector4R boost{ p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ),
	-p12.get( 3 ) };

	// Momentum of p1 & p3 in p12 frame
	const EvtVector4R p1_12{ boostTo( p4_p1, boost ) };
	const EvtVector4R p3_12{ boostTo( p4_p3, boost ) };

	// Cosine of angle between p1 & p3 in p12 frame
	const double p1_12Mag{ p1_12.d3mag() };
	const double p3_12Mag{ p3_12.d3mag() };
	if ( p1_12Mag > 0.0 && p3_12Mag > 0.0 ) {
	value = p1_12.dot( p3_12 ) / ( p1_12Mag * p3_12Mag );
	}
	}

	} else if ( !selectedVarName.compare( "E" ) ) {
	// Energy of first daughter in lab frame
	value = p1_lab.get( 0 );

	} else if ( !selectedVarName.compare( "E_over_Eparent" ) ) {
	// Energy of first daughter w.r.t parent energy (lab frame)
	value = p1_lab.get( 0 ) / selectedParent->getP4Lab().get( 0 );

	} else if ( !selectedVarName.compare( "E_over_Eparent_over05" ) ) {
	// First daughter E_over_Eparent (lab frame) if d2 granddaughter E ratio > 0.5
	const double E_over_Eparent_2{
	parent->getDaug( 0 )->getDaug( d2 - 1 )->getP4Lab().get( 0 ) /
	parent->getDaug( 0 )->getP4Lab().get( 0 ) };

	if ( E_over_Eparent_2 > 0.5 ) {
	value = p1_lab.get( 0 ) / selectedParent->getP4Lab().get( 0 );
	}

	} else if ( !selectedVarName.compare( "totE_over_Mparent" ) ) {
	// Energy of given daughters w.r.t parent mass (lab frame)
	value = ( p1_lab.get( 0 ) + p2_lab.get( 0 ) ) / selectedParent->mass();

	} else if ( !selectedVarName.compare( "prob" ) ) {
	// Decay probability
	const double* dProb{ selectedParent->decayProb() };
	if ( dProb ) {
	value = *dProb;
	}

	} else if ( !selectedVarName.compare( "lifetime" ) ) {
	// Lifetime of particle d1 (ps)
	if ( d1 ) {
	value = par1 != nullptr ? par1->getLifetime() * 1e12 / EvtConst::c
	: 0.0;
	} else {
	value = parent != nullptr
	? parent->getLifetime() * 1e12 / EvtConst::c
	: 0.0;
	}

	} else if ( !selectedVarName.compare( "deltaT" ) ) {
	// Lifetime difference between particles d1 and d2
	const double t1{
	par1 != nullptr ? par1->getLifetime() * 1e12 / EvtConst::c : 0.0 };
	const double t2{
	par2 != nullptr ? par2->getLifetime() * 1e12 / EvtConst::c : 0.0 };
	value = t1 - t2;

	} else if ( !selectedVarName.compare( "decTime" ) ) {
	// Decay flight time of particle d1 in picoseconds.
	// Decay vertex = position of (1st) decay particle of d1
	const EvtParticle* gpar{ par1 != nullptr ? par1->getDaug( 0 ) : nullptr };
	const EvtVector4R vtxPos{ gpar != nullptr ? gpar->get4Pos()
	: EvtVector4R() };
	const double p{ p1_lab.d3mag() };
	value = p > 0.0
	? 1e12 * vtxPos.d3mag() * p1_lab.mass() / ( p * EvtConst::c )
	: 0.0;
	} else if ( !selectedVarName.compare( "nFSRPhotons" ) ) {
	// Loop over all daughters and get number of FSR photons

	double nFSRPhotons{ 0 };

	for ( std::size_t iDaughter{ 0 };
	iDaughter < selectedParent->getNDaug(); iDaughter++ ) {
	const EvtParticle* iDaug = selectedParent->getDaug( iDaughter );

	if ( iDaug->getAttribute( "FSR" ) == 1 )
	nFSRPhotons += 1.0;
	}

	value = nFSRPhotons;

	} else if ( !selectedVarName.compare( "totalFSREnergy" ) ) {
	// Loop over all daughters and get number of FSR photons

	double totalFSREnergy{ 0 };

	for ( std::size_t iDaughter{ 0 };
	iDaughter < selectedParent->getNDaug(); iDaughter++ ) {
	const EvtParticle* iDaug = selectedParent->getDaug( iDaughter );

	if ( iDaug->getAttribute( "FSR" ) == 1 )
	totalFSREnergy += iDaug->getP4Lab().get( 0 );
	}

	value = totalFSREnergy;
	} else {
	std::cerr << "Warning: Did not recognise variable name "
	<< selectedVarName << std::endl;
	}

	return value;
	}

	void TestDecayModel::compareHistos( const TestHistos& theHistos,
	const std::string& refFileName ) const
	{
	// Compare histograms with the same name, calculating the chi-squared
	std::unique_ptr<TFile> refFile{ TFile::Open( refFileName.c_str(), "read" ) };

	if ( !refFile ) {
	std::cerr << "Could not open reference file " << refFileName << std::endl;
	return;
	}

	// TODO - should we plot the (signed) chisq histogram? and save it as pdf/png?

	// TODO - add comparison of mixedHist

	for ( auto& [_, hist] : theHistos.get1DHistos() ) {
	const std::string histName{ hist->GetName() };
	// Get equivalent reference histogram
	const TH1* refHist{
	dynamic_cast<TH1*>( refFile->Get( histName.c_str() ) ) };

	if ( refHist ) {
	double chiSq{ 0.0 };
	int nDof{ 0 };
	int iGood{ 0 };
	const double pValue{
	refHist->Chi2TestX( hist.get(), chiSq, nDof, iGood, "WW" ) };
	const double integral{ refHist->Integral() };
	std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq
	<< "/" << nDof << ", pValue = " << pValue
	<< ", integral = " << integral << std::endl;

	} else {
	std::cerr << "Could not find reference histogram " << histName
	<< std::endl;
	}
	}

	for ( auto& [_, hist] : theHistos.get2DHistos() ) {
	const std::string histName{ hist->GetName() };
	// Get equivalent reference histogram
	const TH2* refHist{
	dynamic_cast<TH2*>( refFile->Get( histName.c_str() ) ) };

	if ( refHist ) {
	double chiSq{ 0.0 };
	int nDof{ 0 };
	int iGood{ 0 };
	const double pValue{
	refHist->Chi2TestX( hist.get(), chiSq, nDof, iGood, "WW" ) };
	const double integral{ refHist->Integral() };
	std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq
	<< "/" << nDof << ", pValue = " << pValue
	<< ", integral = " << integral << std::endl;

	} else {
	std::cerr << "Could not find reference histogram " << histName
	<< std::endl;
	}
	}

	refFile->Close();
	}

	int TestDecayModel::getPartGroup( const int PDGId ) const
	{
	int group( -1 );

	const int absPDGId = std::abs( PDGId );

	if ( absPDGId >= 11 && absPDGId <= 16 ) {
	group = 0; // leptons
	} else if ( absPDGId == 22 ) {
	group = 1; // photon
	} else if ( absPDGId == 211 ) {
	group = 2; // pi+-
	} else if ( absPDGId == 111 ) {
	group = 3; // pi0
	} else if ( absPDGId == 321 ) {
	group = 4; // K+-
	} else if ( absPDGId == 311 \|\| absPDGId == 130 \|\| absPDGId == 310 ) {
	group = 5; // K0
	} else if ( absPDGId == 411 ) {
	group = 6; // D+-
	} else if ( absPDGId == 421 ) {
	group = 7; // D0
	} else if ( absPDGId == 2212 \|\| absPDGId == 2112 \|\| absPDGId == 2224 \|\|
	absPDGId == 2214 \|\| absPDGId == 2114 \|\| absPDGId == 1114 ) {
	group = 8; // light baryons
	} else if ( absPDGId >= 3112 && absPDGId <= 3334 ) {
	group = 9; // strange baryons
	} else if ( absPDGId != 0 ) {
	group = 10; // other particles
	}

	return group;
	}

	bool TestDecayModel::isPartType( const int group,
	const std::string& particleType ) const
	{
	if ( ( !particleType.compare( "LeptonOrPhoton" ) &&
	( group == 0 \|\| group == 1 ) ) \|\|
	( !particleType.compare( "Pion" ) && ( group == 2 \|\| group == 3 ) ) \|\|
	( !particleType.compare( "Kaon" ) && ( group == 4 \|\| group == 5 ) ) \|\|
	( !particleType.compare( "Dmeson" ) && ( group == 6 \|\| group == 7 ) ) \|\|
	( !particleType.compare( "Baryon" ) && ( group == 8 \|\| group == 9 ) ) \|\|
	( !particleType.compare( "Other" ) && group == 10 ) ) {
	return true;
	} else {
	return false;
	}
	}

	double TestDecayModel::getCosAcoplanarityAngle( const EvtParticle* selectedParent,
	const int sel_NDaugMax,
	const int d1, const int d2 ) const
	{
	// Given a two-body decay, the acoplanarity angle is defined as the angle between the
	// two decay planes (normal vectors) in the reference frame of the mother.
	// Each normal vector is the cross product of the momentum of one daughter (in the frame of the mother)
	// and the momentum of one of the granddaughters (in the reference frame of the daughter).
	// In case of 3 daughters, we build an intermediate daughter (2+3) from the last 2 daughters
	// (which are then treated as grand daughters), and in case of 4 daughters, we build
	// 2 intermediate daughters (1+2) and (3+4)

	double cosAco{ 0.0 };

	if ( sel_NDaugMax < 2 \|\| sel_NDaugMax > 4 ) {
	return cosAco;
	}

	const EvtParticle* daughter1{ selectedParent->getDaug( 0 ) };
	const EvtParticle* daughter2{ selectedParent->getDaug( 1 ) };
	const EvtParticle* daughter3{ selectedParent->getDaug( 2 ) };
	const EvtParticle* daughter4{ selectedParent->getDaug( 3 ) };

	const EvtParticle* grandDaughter1{ daughter1->getDaug( d1 - 1 ) };
	const EvtParticle* grandDaughter2{ daughter2->getDaug( d2 - 1 ) };

	const EvtVector4R parent4Vector{ selectedParent->getP4Lab() };

	const EvtVector4R daughter4Vector1{
	sel_NDaugMax <= 3 ? daughter1->getP4Lab()
	: daughter1->getP4Lab() + daughter2->getP4Lab() };

	const EvtVector4R daughter4Vector2{
	sel_NDaugMax == 2 ? daughter2->getP4Lab()
	: sel_NDaugMax == 3 ? daughter2->getP4Lab() + daughter3->getP4Lab()
	: daughter3->getP4Lab() + daughter4->getP4Lab() };

	const EvtVector4R grandDaughter4Vector1{
	sel_NDaugMax <= 3 ? grandDaughter1->getP4Lab() : daughter1->getP4Lab() };

	const EvtVector4R grandDaughter4Vector2{
	sel_NDaugMax == 2 ? grandDaughter2->getP4Lab()
	: sel_NDaugMax == 3 ? daughter2->getP4Lab()
	: daughter3->getP4Lab() };

	const EvtVector4R parentBoost{ parent4Vector.get( 0 ),
	-parent4Vector.get( 1 ),
	-parent4Vector.get( 2 ),
	-parent4Vector.get( 3 ) };

	const EvtVector4R daughter1Boost{ daughter4Vector1.get( 0 ),
	-daughter4Vector1.get( 1 ),
	-daughter4Vector1.get( 2 ),
	-daughter4Vector1.get( 3 ) };

	const EvtVector4R daughter2Boost{ daughter4Vector2.get( 0 ),
	-daughter4Vector2.get( 1 ),
	-daughter4Vector2.get( 2 ),
	-daughter4Vector2.get( 3 ) };

	// Boosting daughters to reference frame of the mother
	const EvtVector4R daughter4Vector1_boosted{
	boostTo( daughter4Vector1, parentBoost ) };
	const EvtVector4R daughter4Vector2_boosted{
	boostTo( daughter4Vector2, parentBoost ) };

	// Boosting each granddaughter to reference frame of its mother
	const EvtVector4R grandDaughter4Vector1_boosted{
	boostTo( grandDaughter4Vector1, daughter1Boost ) };
	const EvtVector4R grandDaughter4Vector2_boosted{
	boostTo( grandDaughter4Vector2, daughter2Boost ) };

	// We calculate the normal vectors of the decay two planes
	const EvtVector4R normalVector1{
	daughter4Vector1_boosted.cross( grandDaughter4Vector1_boosted ) };
	const EvtVector4R normalVector2{
	daughter4Vector2_boosted.cross( grandDaughter4Vector2_boosted ) };

	const double normalVector1Mag{ normalVector1.d3mag() };
	const double normalVector2Mag{ normalVector2.d3mag() };

	if ( normalVector1Mag > 0.0 && normalVector2Mag > 0.0 ) {
	cosAco = normalVector1.dot( normalVector2 ) /
	( normalVector1Mag * normalVector2Mag );
	}

	return cosAco;
	}

	int main( int argc, char* argv[] )
	{
	if ( argc < 2 \|\| argc > 3 ) {
	std::cerr << "Expecting at least one argument: test configuration json file"
	<< "\nAdditional argument supported for: general configuration json file"
	<< std::endl;
	return 1;
	}

	// Tweak ROOT behaviour
	ROOT::EnableThreadSafety();
	TH1::AddDirectory( kFALSE );

	// Load input file in json format
	const std::string testConfigFileName{ argv[1] };
	const std::string generalConfigFileName{
	( argc > 2 ) ? argv[2] : "jsonFiles/config/default.json" };

	json generalConfig;
	json testConfig;

	{
	std::ifstream inputStr{ generalConfigFileName };
	inputStr >> generalConfig;
	}

	{
	std::ifstream inputStr{ testConfigFileName };
	inputStr >> testConfig;
	}

	if ( testConfig.is_array() ) {
	for ( auto& cc : testConfig ) {
	cc.merge_patch( generalConfig );
	TestDecayModel test{ cc };
	test.run();
	}
	} else {
	testConfig.merge_patch( generalConfig );
	TestDecayModel test{ testConfig };
	test.run();
	}

	return 0;
	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:28 PM (22 h, 37 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4437442
Default Alt Text: (192 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions