D112.1759146623.diff
No OneTemporary
Actions

Size

35 KB

Referenced Files

None

Subscribers

None

D112.1759146623.diff
View Options

	diff --git a/EvtGenBase/EvtAbsRadCorr.hh b/EvtGenBase/EvtAbsRadCorr.hh
	--- a/EvtGenBase/EvtAbsRadCorr.hh
	+++ b/EvtGenBase/EvtAbsRadCorr.hh
	@@ -29,7 +29,8 @@
	public:
	EvtAbsRadCorr(){};
	virtual ~EvtAbsRadCorr(){};
	- virtual void doRadCorr( EvtParticle* p ) = 0;
	+ virtual void initialise() = 0;
	+ virtual void doRadCorr( EvtParticle* p ) const = 0;

	private:
	};
	diff --git a/EvtGenExternal/EvtExternalGenFactory.hh b/EvtGenExternal/EvtExternalGenFactory.hh
	--- a/EvtGenExternal/EvtExternalGenFactory.hh
	+++ b/EvtGenExternal/EvtExternalGenFactory.hh
	@@ -33,7 +33,6 @@
	enum genId
	{
	PythiaGenId = 0,
	- PhotosGenId,
	TauolaGenId
	};

	@@ -45,8 +44,6 @@

	void definePythiaGenerator( std::string xmlDir, bool convertPhysCodes,
	bool useEvtGenRandom = true );
	- void definePhotosGenerator( std::string photonType = "gamma",
	- bool useEvtGenRandom = true );
	void defineTauolaGenerator( bool useEvtGenRandom = true );

	//methods to add configuration commands to the pythia generators
	diff --git a/EvtGenExternal/EvtExternalGenList.hh b/EvtGenExternal/EvtExternalGenList.hh
	--- a/EvtGenExternal/EvtExternalGenList.hh
	+++ b/EvtGenExternal/EvtExternalGenList.hh
	@@ -40,10 +40,14 @@

	std::list<EvtDecayBase*> getListOfModels();

	- EvtAbsRadCorr* getPhotosModel();
	+ EvtAbsRadCorr* getPhotosModel( const double infraredCutOff = 1.0e-7,
	+ const double maxWtInterference = 64.0 );

	protected:
	private:
	+ std::string m_photonType;
	+
	+ bool m_useEvtGenRandom;
	};

	#endif
	diff --git a/EvtGenExternal/EvtPHOTOS.hh b/EvtGenExternal/EvtPHOTOS.hh
	--- a/EvtGenExternal/EvtPHOTOS.hh
	+++ b/EvtGenExternal/EvtPHOTOS.hh
	@@ -22,21 +22,64 @@
	#define EVTPHOTOS_HH

	#include "EvtGenBase/EvtAbsRadCorr.hh"
	+#include "EvtGenBase/EvtHepMCEvent.hh"
	+#include "EvtGenBase/EvtId.hh"
	+#include "EvtGenBase/EvtParticle.hh"
	+#include "EvtGenBase/EvtVector4R.hh"

	+#ifdef EVTGEN_PHOTOS
	+#ifdef EVTGEN_HEPMC3
	+#include "HepMC3/Units.h"
	+
	+#include "Photos/PhotosHepMC3Event.h"
	+#include "Photos/PhotosHepMC3Particle.h"
	+#else
	+#include "Photos/PhotosHepMCEvent.h"
	+#include "Photos/PhotosHepMCParticle.h"
	+#include "Photos/PhotosParticle.h"
	+#endif
	+#endif
	+
	+#include <mutex>
	#include <string>

	class EvtParticle;
	class EvtAbsExternalGen;

	-// Description: EvtGen's interface to PHOTOS for generation of
	-// QED final state radiation.
	+/*
	+ * Description: EvtGen's interface to PHOTOS for generation of
	+ * QED final-state radiation.
	+ */

	class EvtPHOTOS : public EvtAbsRadCorr {
	public:
	- void doRadCorr( EvtParticle* p ) override;
	+ EvtPHOTOS( const std::string& photonType = "gamma",
	+ const bool useEvtGenRandom = true,
	+ const double infraredCutOff = 1.0e-7,
	+ const double maxWtInterference = 64.0 );
	+
	+ void initialise() override;
	+
	+ void doRadCorr( EvtParticle* theParticle ) const override;

	private:
	- EvtAbsExternalGen* _photosEngine = nullptr;
	+#ifdef EVTGEN_PHOTOS
	+ GenParticlePtr createGenParticle( const EvtParticle& theParticle,
	+ bool incoming ) const;
	+
	+ int getNumberOfPhotons( const GenVertexPtr theVertex ) const;
	+
	+ // Default photon type, id, pdg number and mass
	+ std::string m_photonType = "gamma";
	+ EvtId m_gammaId = EvtId( -1, -1 );
	+ int m_gammaPDG = 22;
	+ double m_mPhoton = 0.0;
	+
	+ bool m_initialised = false;
	+
	+ static std::mutex photos_mutex;
	+
	+#endif
	};

	#endif
	diff --git a/EvtGenExternal/EvtPhotosEngine.hh b/EvtGenExternal/EvtPhotosEngine.hh
	deleted file mode 100644
	--- a/EvtGenExternal/EvtPhotosEngine.hh
	+++ /dev/null
	@@ -1,68 +0,0 @@
	-
	-/***********************************************************************
	-* 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_PHOTOS
	-#ifndef EVTPHOTOSENGINE_HH
	-#define EVTPHOTOSENGINE_HH
	-
	-#include "EvtGenBase/EvtHepMCEvent.hh"
	-#include "EvtGenBase/EvtId.hh"
	-#include "EvtGenBase/EvtParticle.hh"
	-#include "EvtGenBase/EvtVector4R.hh"
	-
	-#include "EvtGenModels/EvtAbsExternalGen.hh"
	-#ifdef EVTGEN_HEPMC3
	-#include "HepMC3/Units.h"
	-
	-#include "Photos/PhotosHepMC3Event.h"
	-#include "Photos/PhotosHepMC3Particle.h"
	-#else
	-#include "Photos/PhotosHepMCEvent.h"
	-#include "Photos/PhotosHepMCParticle.h"
	-#include "Photos/PhotosParticle.h"
	-#endif
	-
	-#include <string>
	-
	-// Description: Interface to the PHOTOS external generator
	-
	-class EvtPhotosEngine : public EvtAbsExternalGen {
	- public:
	- EvtPhotosEngine( std::string photonType = "gamma",
	- bool useEvtGenRandom = true );
	-
	- bool doDecay( EvtParticle* theMother ) override;
	-
	- void initialise() override;
	-
	- private:
	- std::string _photonType;
	- EvtId _gammaId;
	- int _gammaPDG;
	- double _mPhoton;
	- bool _initialised;
	-
	- GenParticlePtr createGenParticle( EvtParticle* theParticle, bool incoming );
	- int getNumberOfPhotons( const GenVertexPtr theVertex ) const;
	-};
	-
	-#endif
	-
	-#endif
	diff --git a/EvtGenModels/EvtNoRadCorr.hh b/EvtGenModels/EvtNoRadCorr.hh
	--- a/EvtGenModels/EvtNoRadCorr.hh
	+++ b/EvtGenModels/EvtNoRadCorr.hh
	@@ -36,7 +36,9 @@
	EvtNoRadCorr() { ; }
	virtual ~EvtNoRadCorr() { ; }

	- virtual void doRadCorr( EvtParticle* ) { ; }
	+ virtual void initialise() { ; }
	+
	+ virtual void doRadCorr( EvtParticle* ) const { ; }

	private:
	};
	diff --git a/History.md b/History.md
	--- a/History.md
	+++ b/History.md
	@@ -11,6 +11,11 @@
	===
	## R02-0X-00

	+12 Apr 2024 Fernando Abudinen
	+* D112: Removed EvtPhotosEngine and moved its functionality EvtPHOTOS class.
	+ Implemented const correctness in doRadCorr function of EvtAbsRadCorr and EvtNoRadCorr.
	+ This required moving initialisation of fsrEngine to setRadCorrEngine.
	+
	11 Apr 2024 Thomas Latham
	* D109: Remove broken or obsolete models
	- Models removed:
	diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
	--- a/src/CMakeLists.txt
	+++ b/src/CMakeLists.txt
	@@ -117,6 +117,9 @@
	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()
	@@ -150,6 +153,7 @@
	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})
	@@ -170,6 +174,7 @@
	# 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")
	@@ -206,6 +211,7 @@
	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})
	@@ -226,6 +232,7 @@
	# 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")
	diff --git a/src/EvtGenBase/EvtRadCorr.cpp b/src/EvtGenBase/EvtRadCorr.cpp
	--- a/src/EvtGenBase/EvtRadCorr.cpp
	+++ b/src/EvtGenBase/EvtRadCorr.cpp
	@@ -49,6 +49,10 @@
	void EvtRadCorr::setRadCorrEngine( EvtAbsRadCorr* fsrEngine )
	{
	_fsrEngine = fsrEngine;
	+
	+ if ( _fsrEngine ) {
	+ _fsrEngine->initialise();
	+ }
	}

	void EvtRadCorr::doRadCorr( EvtParticle* p )
	diff --git a/src/EvtGenExternal/EvtExternalGenFactory.cpp b/src/EvtGenExternal/EvtExternalGenFactory.cpp
	--- a/src/EvtGenExternal/EvtExternalGenFactory.cpp
	+++ b/src/EvtGenExternal/EvtExternalGenFactory.cpp
	@@ -27,10 +27,6 @@
	#include "EvtGenExternal/EvtPythiaEngine.hh"
	#endif

	-#ifdef EVTGEN_PHOTOS
	-#include "EvtGenExternal/EvtPhotosEngine.hh"
	-#endif
	-
	#ifdef EVTGEN_TAUOLA
	#include "EvtGenExternal/EvtTauolaEngine.hh"
	#endif
	@@ -99,25 +95,6 @@
	}
	#endif

	-#ifdef EVTGEN_PHOTOS
	-void EvtExternalGenFactory::definePhotosGenerator( std::string photonType,
	- bool useEvtGenRandom )
	-{
	- int genId = EvtExternalGenFactory::PhotosGenId;
	-
	- EvtGenReport( EVTGEN_INFO, "EvtGen" )
	- << "Defining EvtPhotosEngine using photonType = " << photonType << endl;
	-
	- EvtAbsExternalGen* photosGenerator = new EvtPhotosEngine( photonType,
	- useEvtGenRandom );
	- _extGenMap[genId] = photosGenerator;
	-}
	-#else
	-void EvtExternalGenFactory::definePhotosGenerator( std::string, bool )
	-{
	-}
	-#endif
	-
	#ifdef EVTGEN_TAUOLA
	void EvtExternalGenFactory::defineTauolaGenerator( bool useEvtGenRandom )
	{
	diff --git a/src/EvtGenExternal/EvtExternalGenList.cpp b/src/EvtGenExternal/EvtExternalGenList.cpp
	--- a/src/EvtGenExternal/EvtExternalGenList.cpp
	+++ b/src/EvtGenExternal/EvtExternalGenList.cpp
	@@ -28,14 +28,12 @@
	EvtExternalGenList::EvtExternalGenList( bool convertPythiaCodes,
	std::string pythiaXmlDir,
	std::string photonType,
	- bool useEvtGenRandom )
	+ bool useEvtGenRandom ) :
	+ m_photonType{ photonType }, m_useEvtGenRandom{ useEvtGenRandom }
	{
	// Instantiate the external generator factory
	EvtExternalGenFactory* extFactory = EvtExternalGenFactory::getInstance();

	- // Define the external generator "engines" here
	- extFactory->definePhotosGenerator( photonType, useEvtGenRandom );
	-
	if ( pythiaXmlDir.size() < 1 ) {
	// If we have no string defined, check the value of the
	// PYTHIA8DATA environment variable which should be set to the
	@@ -56,10 +54,12 @@
	{
	}

	-EvtAbsRadCorr* EvtExternalGenList::getPhotosModel()
	+EvtAbsRadCorr* EvtExternalGenList::getPhotosModel( const double infraredCutOff,
	+ const double maxWtInterference )
	{
	- // Define the Photos model, which uses the EvtPhotosEngine class.
	- EvtPHOTOS* photosModel = new EvtPHOTOS();
	+ // Define the Photos model, which uses the EvtPHOTOS class.
	+ EvtPHOTOS* photosModel = new EvtPHOTOS( m_photonType, m_useEvtGenRandom,
	+ infraredCutOff, maxWtInterference );
	return photosModel;
	}

	diff --git a/src/EvtGenExternal/EvtPHOTOS.cpp b/src/EvtGenExternal/EvtPHOTOS.cpp
	--- a/src/EvtGenExternal/EvtPHOTOS.cpp
	+++ b/src/EvtGenExternal/EvtPHOTOS.cpp
	@@ -20,18 +20,309 @@

	#include "EvtGenExternal/EvtPHOTOS.hh"

	-#include "EvtGenBase/EvtPatches.hh"
	+#include "EvtGenBase/EvtPDL.hh"
	+#include "EvtGenBase/EvtPhotonParticle.hh"
	+#include "EvtGenBase/EvtRandom.hh"
	+#include "EvtGenBase/EvtReport.hh"
	+#include "EvtGenBase/EvtVector4R.hh"

	-#include "EvtGenExternal/EvtExternalGenFactory.hh"
	+#include <iostream>
	+#include <sstream>
	+#include <vector>

	-void EvtPHOTOS::doRadCorr( EvtParticle* p )
	+using std::endl;
	+
	+#ifdef EVTGEN_PHOTOS
	+
	+// Mutex PHOTOS as it is not thread safe.
	+std::mutex EvtPHOTOS::photos_mutex;
	+
	+EvtPHOTOS::EvtPHOTOS( const std::string& photonType, const bool useEvtGenRandom,
	+ const double infraredCutOff,
	+ const double maxWtInterference ) :
	+ m_photonType{ photonType }
	+{
	+ photos_mutex.lock();
	+
	+ EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Setting up PHOTOS." << endl;
	+
	+ if ( useEvtGenRandom ) {
	+ EvtGenReport( EVTGEN_INFO, "EvtGen" )
	+ << "Using EvtGen random number engine also for Photos++" << endl;
	+
	+ Photospp::Photos::setRandomGenerator( EvtRandom::Flat );
	+ }
	+
	+ Photospp::Photos::initialize();
	+
	+ Photospp::Photos::setExponentiation( true );
	+
	+ /* Sets the infrared cutoff (default at 1e-7)
	+ * Reset the minimum photon energy, if required, in units of half of the decaying particle mass.
	+ * This must be done after exponentiation! Keep the cut at 1e-7, i.e. 0.1 keV at the 1 GeV scale,
	+ * which is appropriate for B decays
	+ */
	+ Photospp::Photos::setInfraredCutOff( infraredCutOff );
	+
	+ Photospp::Photos::setInterference( true );
	+
	+ /* Increase the maximum possible value of the interference weight
	+ * corresponding to 2^n, where n = number of charges (+,-).
	+ */
	+ Photospp::Photos::maxWtInterference( maxWtInterference );
	+
	+#ifdef EVTGEN_PHOTOS_NEWLIBS
	+ // This turns off/on virtual photons splitting into l+l-
	+ // It is false by default, but just to keep track of it here.
	+ Photospp::Photos::setPairEmission( false );
	+#endif
	+
	+ photos_mutex.unlock();
	+}
	+#else
	+EvtPHOTOS::EvtPHOTOS( const std::string& /photonType/,
	+ const bool /useEvtGenRandom/,
	+ const double /infraredCutOff/,
	+ const double /maxWtInterference/ )
	+{
	+ EvtGenReport( EVTGEN_WARNING, "EvtGen" )
	+ << " PHOTOS has been called for FSR simulation, but it was not switched on during compilation."
	+ << endl;
	+
	+ EvtGenReport( EVTGEN_WARNING, "EvtGen" )
	+ << " The simulation will be generated without FSR." << endl;
	+}
	+#endif
	+
	+#ifdef EVTGEN_PHOTOS
	+void EvtPHOTOS::doRadCorr( EvtParticle* theParticle ) const
	+{
	+ if ( !theParticle ) {
	+ return;
	+ }
	+
	+ /* Skip running Photos if the particle has no daughters, since we can't add FSR.
	+ * Also skip Photos if the particle has too many daughters (>= 10) to avoid a problem
	+ * with a hard coded upper limit in the PHOENE subroutine.
	+ */
	+ const int nDaughters( theParticle->getNDaug() );
	+ if ( nDaughters == 0 \|\| nDaughters >= 10 ) {
	+ return;
	+ }
	+
	+ /* Create a dummy HepMC GenEvent containing a single vertex, with the mother
	+ * assigned as the incoming particle and its daughters as outgoing particles.
	+ * We then pass this event to Photos for processing.
	+ * It will return a modified version of the event, updating the momentum of
	+ * the original particles and will contain any new photon particles.
	+ * We add these extra photons to the mother particle daughter list.
	+ * First, Create the dummy event.
	+ */
	+ auto theEvent = std::make_unique<GenEvent>( Units::GEV, Units::MM );
	+
	+ // Create the decay "vertex".
	+ GenVertexPtr theVertex = newGenVertexPtr();
	+ theEvent->add_vertex( theVertex );
	+
	+ // Add the mother particle as the incoming particle to the vertex.
	+ const GenParticlePtr hepMCMother = this->createGenParticle( *theParticle,
	+ true );
	+ theVertex->add_particle_in( hepMCMother );
	+
	+ /* Find all daughter particles and assign them as outgoing particles to the vertex.
	+ * Keep track of the number of photons already in the decay (e.g. we may have B -> K* gamma)
	+ */
	+ int iDaug( 0 ), nGamma( 0 );
	+ for ( iDaug = 0; iDaug < nDaughters; iDaug++ ) {
	+ const EvtParticle& theDaughter = *theParticle->getDaug( iDaug );
	+ const GenParticlePtr hepMCDaughter =
	+ this->createGenParticle( theDaughter, false );
	+ theVertex->add_particle_out( hepMCDaughter );
	+
	+ if ( theDaughter.getPDGId() == m_gammaPDG ) {
	+ nGamma++;
	+ }
	+ }
	+
	+ photos_mutex.lock();
	+
	+ /* Now pass the event to Photos for processing
	+ * Create a Photos event object */
	+#ifdef EVTGEN_HEPMC3
	+ Photospp::PhotosHepMC3Event photosEvent( theEvent.get() );
	+#else
	+ Photospp::PhotosHepMCEvent photosEvent( theEvent.get() );
	+#endif
	+
	+ // Run the Photos algorithm
	+ photosEvent.process();
	+
	+ photos_mutex.unlock();
	+
	+ // Find the number of (outgoing) photons in the event
	+ const int nPhotons = this->getNumberOfPhotons( theVertex );
	+
	+ int iLoop( 0 );
	+ // See whether Photos has created additional photons. If not, do nothing extra.
	+ if ( nPhotons > nGamma ) {
	+ // We have extra particles from Photos; these would have been appended
	+ // to the outgoing particle list
	+
	+ // Get the iterator of outgoing particles for this vertex
	+#ifdef EVTGEN_HEPMC3
	+ for ( const auto& outParticle : theVertex->particles_out() ) {
	+#else
	+ HepMC::GenVertex::particles_out_const_iterator outIter;
	+ for ( outIter = theVertex->particles_out_const_begin();
	+ outIter != theVertex->particles_out_const_end(); ++outIter ) {
	+ // Get the next HepMC GenParticle
	+ const HepMC::GenParticle* outParticle = *outIter;
	+#endif
	+ // Get the three-momentum Photos result for this particle, and the PDG id
	+ if ( outParticle ) {
	+ const FourVector HepMCP4 = outParticle->momentum();
	+ const double px = HepMCP4.px();
	+ const double py = HepMCP4.py();
	+ const double pz = HepMCP4.pz();
	+ const double energy = HepMCP4.e();
	+ const int pdgId = outParticle->pdg_id();
	+ const EvtVector4R newP4( energy, px, py, pz );
	+
	+ if ( iLoop < nDaughters ) {
	+ /* Original daughters: keep the original particle mass,
	+ * but set the three-momentum according to what Photos
	+ * has modified.
	+ */
	+ EvtParticle* daugParticle = theParticle->getDaug( iLoop );
	+ if ( daugParticle ) {
	+ daugParticle->setP4WithFSR( newP4 );
	+ }
	+
	+ } else if ( pdgId == m_gammaPDG ) {
	+ // Extra photon particle. Setup the four-momentum object
	+
	+ // Create a new photon particle and add it to the list of daughters
	+ EvtPhotonParticle* gamma = new EvtPhotonParticle();
	+ gamma->init( m_gammaId, newP4 );
	+ // Set the pre-FSR photon momentum to zero
	+ gamma->setFSRP4toZero();
	+ // Set its particle attribute to specify it is a FSR photon
	+ gamma->setAttribute( "FSR", 1 ); // it is a FSR photon
	+ gamma->setAttribute( "ISR", 0 ); // it is not an ISR photon
	+ // Let the mother know about this new photon
	+ gamma->addDaug( theParticle );
	+ }
	+
	+ // Increment the loop counter for detecting additional photon particles
	+ iLoop++;
	+ }
	+ }
	+ }
	+
	+ // Cleanup
	+ theEvent->clear();
	+
	+ return;
	+}
	+#else
	+void EvtPHOTOS::doRadCorr( EvtParticle* /theParticle/ ) const
	{
	- if ( !_photosEngine ) {
	- _photosEngine = EvtExternalGenFactory::getInstance()->getGenerator(
	- EvtExternalGenFactory::PhotosGenId );
	+}
	+#endif
	+
	+#ifdef EVTGEN_PHOTOS
	+void EvtPHOTOS::initialise()
	+{
	+ if ( m_initialised ) {
	+ return;
	}
	+ m_gammaId = EvtPDL::getId( m_photonType );

	- if ( _photosEngine ) {
	- _photosEngine->doDecay( p );
	+ if ( m_gammaId == EvtId( -1, -1 ) ) {
	+ EvtGenReport( EVTGEN_INFO, "EvtGen" )
	+ << "Error in EvtPHOTOS. Do not recognise the photon type "
	+ << m_photonType << ". Setting this to \"gamma\". " << endl;
	+ m_gammaId = EvtPDL::getId( "gamma" );
	}
	+
	+ m_gammaPDG = EvtPDL::getStdHep( m_gammaId );
	+ m_mPhoton = EvtPDL::getMeanMass( m_gammaId );
	+
	+ m_initialised = true;
	+}
	+#else
	+void EvtPHOTOS::initialise()
	+{
	}
	+#endif
	+
	+#ifdef EVTGEN_PHOTOS
	+GenParticlePtr EvtPHOTOS::createGenParticle( const EvtParticle& theParticle,
	+ bool incoming ) const
	+{
	+ /* Method to create an HepMC::GenParticle version of the given EvtParticle.
	+ */
	+
	+ // Get the 4-momentum (E, px, py, pz) for the EvtParticle
	+ const EvtVector4R p4 = incoming ? theParticle.getP4Restframe()
	+ : theParticle.getP4();
	+
	+ // Convert this to the HepMC 4-momentum
	+ const double E = p4.get( 0 );
	+ const double px = p4.get( 1 );
	+ const double py = p4.get( 2 );
	+ const double pz = p4.get( 3 );
	+
	+ const FourVector hepMC_p4( px, py, pz, E );
	+
	+ const int PDGInt = EvtPDL::getStdHep( theParticle.getId() );
	+
	+ // Set the status flag for the particle. This is required, otherwise Photos++
	+ // will crash from out-of-bounds array index problems.
	+ const int status = incoming ? Photospp::PhotosParticle::HISTORY
	+ : Photospp::PhotosParticle::STABLE;
	+
	+ const GenParticlePtr genParticle = newGenParticlePtr( hepMC_p4, PDGInt,
	+ status );
	+
	+ return genParticle;
	+}
	+
	+int EvtPHOTOS::getNumberOfPhotons( const GenVertexPtr theVertex ) const
	+{
	+ /* Find the number of photons from the outgoing particle list
	+ */
	+
	+ if ( !theVertex ) {
	+ return 0;
	+ }
	+
	+ int nPhotons( 0 );
	+
	+ // Get the iterator of outgoing particles for this vertex
	+#ifdef EVTGEN_HEPMC3
	+ for ( const auto& outParticle : theVertex->particles_out() ) {
	+#else
	+ HepMC::GenVertex::particles_out_const_iterator outIter;
	+ for ( outIter = theVertex->particles_out_const_begin();
	+ outIter != theVertex->particles_out_const_end(); ++outIter ) {
	+ // Get the next HepMC GenParticle
	+ const HepMC::GenParticle* outParticle = *outIter;
	+#endif
	+
	+ // Get the PDG id
	+ int pdgId( 0 );
	+ if ( outParticle != nullptr ) {
	+ pdgId = outParticle->pdg_id();
	+ }
	+
	+ // Keep track of how many photons there are
	+ if ( pdgId == m_gammaPDG ) {
	+ nPhotons++;
	+ }
	+ }
	+
	+ return nPhotons;
	+}
	+
	+#endif
	diff --git a/src/EvtGenExternal/EvtPhotosEngine.cpp b/src/EvtGenExternal/EvtPhotosEngine.cpp
	deleted file mode 100644
	--- a/src/EvtGenExternal/EvtPhotosEngine.cpp
	+++ /dev/null
	@@ -1,307 +0,0 @@
	-
	-/***********************************************************************
	-* 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_PHOTOS
	-
	-#include "EvtGenExternal/EvtPhotosEngine.hh"
	-
	-#include "EvtGenBase/EvtPDL.hh"
	-#include "EvtGenBase/EvtPhotonParticle.hh"
	-#include "EvtGenBase/EvtRandom.hh"
	-#include "EvtGenBase/EvtReport.hh"
	-#include "EvtGenBase/EvtVector4R.hh"
	-
	-#include <iostream>
	-#include <sstream>
	-#include <vector>
	-
	-using std::endl;
	-
	-EvtPhotosEngine::EvtPhotosEngine( std::string photonType, bool useEvtGenRandom )
	-{
	- _photonType = photonType;
	- _gammaId = EvtId( -1, -1 );
	- _gammaPDG = 22; // default photon pdg integer
	- _mPhoton = 0.0;
	-
	- EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Setting up PHOTOS." << endl;
	-
	- if ( useEvtGenRandom == true ) {
	- EvtGenReport( EVTGEN_INFO, "EvtGen" )
	- << "Using EvtGen random number engine also for Photos++" << endl;
	-
	- Photospp::Photos::setRandomGenerator( EvtRandom::Flat );
	- }
	-
	- Photospp::Photos::initialize();
	-
	- // Increase the maximum possible value of the interference weight
	- Photospp::Photos::maxWtInterference(
	- 64.0 ); // 2^n, where n = number of charges (+,-)
	- Photospp::Photos::setInterference( true );
	- Photospp::Photos::setExponentiation( true ); // Sets the infrared cutoff at 1e-7
	- // Reset the minimum photon energy, if required, in units of half of the decaying particle mass.
	- // This must be done after exponentiation! Keep the cut at 1e-7, i.e. 0.1 keV at the 1 GeV scale,
	- // which is appropriate for B decays
	- Photospp::Photos::setInfraredCutOff( 1.0e-7 );
	-
	- _initialised = false;
	-}
	-
	-void EvtPhotosEngine::initialise()
	-{
	- if ( _initialised == false ) {
	- _gammaId = EvtPDL::getId( _photonType );
	-
	- if ( _gammaId == EvtId( -1, -1 ) ) {
	- EvtGenReport( EVTGEN_INFO, "EvtGen" )
	- << "Error in EvtPhotosEngine. Do not recognise the photon type "
	- << _photonType << ". Setting this to \"gamma\". " << endl;
	- _gammaId = EvtPDL::getId( "gamma" );
	- }
	-
	- _gammaPDG = EvtPDL::getStdHep( _gammaId );
	- _mPhoton = EvtPDL::getMeanMass( _gammaId );
	-
	- _initialised = true;
	- }
	-}
	-
	-bool EvtPhotosEngine::doDecay( EvtParticle* theMother )
	-{
	- if ( _initialised == false ) {
	- this->initialise();
	- }
	-
	- if ( theMother == nullptr ) {
	- return false;
	- }
	-
	- // Create a dummy HepMC GenEvent containing a single vertex, with the mother
	- // assigned as the incoming particle and its daughters as outgoing particles.
	- // We then pass this event to Photos for processing.
	- // It will return a modified version of the event, updating the momentum of
	- // the original particles and will contain any new photon particles.
	- // We add these extra photons to the mother particle daughter list.
	-
	- // Skip running Photos if the particle has no daughters, since we can't add FSR.
	- // Also skip Photos if the particle has too many daughters (>= 10) to avoid a problem
	- // with a hard coded upper limit in the PHOENE subroutine.
	- int nDaug( theMother->getNDaug() );
	- if ( nDaug == 0 \|\| nDaug >= 10 ) {
	- return false;
	- }
	-
	- // Create the dummy event.
	- auto theEvent = std::make_unique<GenEvent>( Units::GEV, Units::MM );
	-
	- // Create the decay "vertex".
	- GenVertexPtr theVertex = newGenVertexPtr();
	- theEvent->add_vertex( theVertex );
	-
	- // Add the mother particle as the incoming particle to the vertex.
	- GenParticlePtr hepMCMother = this->createGenParticle( theMother, true );
	- theVertex->add_particle_in( hepMCMother );
	-
	- // Find all daughter particles and assign them as outgoing particles to the vertex.
	- // Keep track of the number of photons already in the decay (e.g. we may have B -> K* gamma)
	- int iDaug( 0 ), nGamma( 0 );
	- for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
	- EvtParticle* theDaughter = theMother->getDaug( iDaug );
	- GenParticlePtr hepMCDaughter = this->createGenParticle( theDaughter,
	- false );
	- theVertex->add_particle_out( hepMCDaughter );
	-
	- if ( theDaughter ) {
	- int daugId = theDaughter->getPDGId();
	- if ( daugId == _gammaPDG ) {
	- nGamma++;
	- }
	- }
	- }
	-
	- // Now pass the event to Photos for processing
	- // Create a Photos event object
	-#ifdef EVTGEN_HEPMC3
	- Photospp::PhotosHepMC3Event photosEvent( theEvent.get() );
	-#else
	- Photospp::PhotosHepMCEvent photosEvent( theEvent.get() );
	-#endif
	-
	- // Run the Photos algorithm
	- photosEvent.process();
	-
	- // Find the number of (outgoing) photons in the event
	- int nPhotons = this->getNumberOfPhotons( theVertex );
	-
	- // See if Photos has created additional photons. If not, do nothing extra
	- int nDiffPhotons = nPhotons - nGamma;
	- int iLoop( 0 );
	-
	- if ( nDiffPhotons > 0 ) {
	- // We have extra particles from Photos; these would have been appended
	- // to the outgoing particle list
	-
	- // Get the iterator of outgoing particles for this vertex
	-#ifdef EVTGEN_HEPMC3
	- for ( auto outParticle : theVertex->particles_out() ) {
	-#else
	- HepMC::GenVertex::particles_out_const_iterator outIter;
	- for ( outIter = theVertex->particles_out_const_begin();
	- outIter != theVertex->particles_out_const_end(); ++outIter ) {
	- // Get the next HepMC GenParticle
	- HepMC::GenParticle* outParticle = *outIter;
	-#endif
	-
	- // Get the three-momentum Photos result for this particle, and the PDG id
	- double px( 0.0 ), py( 0.0 ), pz( 0.0 );
	- int pdgId( 0 );
	-
	- if ( outParticle != nullptr ) {
	- FourVector HepMCP4 = outParticle->momentum();
	- px = HepMCP4.px();
	- py = HepMCP4.py();
	- pz = HepMCP4.pz();
	- pdgId = outParticle->pdg_id();
	- }
	-
	- // Create an empty 4-momentum vector for the new/modified daughters
	- EvtVector4R newP4;
	-
	- if ( iLoop < nDaug ) {
	- // Original daughters
	- EvtParticle* daugParticle = theMother->getDaug( iLoop );
	- if ( daugParticle != nullptr ) {
	- // Keep the original particle mass, but set the three-momentum
	- // according to what Photos has modified. However, this will
	- // violate energy conservation (from what Photos has provided).
	- double mass = daugParticle->mass();
	- double energy = sqrt( mass * mass + px * px + py * py +
	- pz * pz );
	- newP4.set( energy, px, py, pz );
	- // Set the new four-momentum (FSR applied)
	- daugParticle->setP4WithFSR( newP4 );
	- }
	-
	- } else if ( pdgId == _gammaPDG ) {
	- // Extra photon particle. Setup the four-momentum object
	- double energy = sqrt( _mPhoton * _mPhoton + px * px + py * py +
	- pz * pz );
	- newP4.set( energy, px, py, pz );
	-
	- // Create a new photon particle and add it to the list of daughters
	- EvtPhotonParticle* gamma = new EvtPhotonParticle();
	- gamma->init( _gammaId, newP4 );
	- // Set the pre-FSR photon momentum to zero
	- gamma->setFSRP4toZero();
	- // Let the mother know about this new photon
	- gamma->addDaug( theMother );
	- // Set its particle attribute to specify it is a FSR photon
	- gamma->setAttribute( "FSR", 1 ); // it is a FSR photon
	- gamma->setAttribute( "ISR", 0 ); // it is not an ISR photon
	- }
	-
	- // Increment the loop counter for detecting additional photon particles
	- iLoop++;
	- }
	- }
	-
	- // Cleanup
	- theEvent->clear();
	-
	- return true;
	-}
	-
	-GenParticlePtr EvtPhotosEngine::createGenParticle( EvtParticle* theParticle,
	- bool incoming )
	-{
	- // Method to create an HepMC::GenParticle version of the given EvtParticle.
	- if ( theParticle == nullptr ) {
	- return nullptr;
	- }
	-
	- // Get the 4-momentum (E, px, py, pz) for the EvtParticle
	- EvtVector4R p4( 0.0, 0.0, 0.0, 0.0 );
	-
	- if ( incoming == true ) {
	- p4 = theParticle->getP4Restframe();
	- } else {
	- p4 = theParticle->getP4();
	- }
	-
	- // Convert this to the HepMC 4-momentum
	- double E = p4.get( 0 );
	- double px = p4.get( 1 );
	- double py = p4.get( 2 );
	- double pz = p4.get( 3 );
	-
	- FourVector hepMC_p4( px, py, pz, E );
	-
	- int PDGInt = EvtPDL::getStdHep( theParticle->getId() );
	-
	- // Set the status flag for the particle. This is required, otherwise Photos++
	- // will crash from out-of-bounds array index problems.
	- int status = Photospp::PhotosParticle::HISTORY;
	- if ( incoming == false ) {
	- status = Photospp::PhotosParticle::STABLE;
	- }
	-
	- GenParticlePtr genParticle = newGenParticlePtr( hepMC_p4, PDGInt, status );
	-
	- return genParticle;
	-}
	-
	-int EvtPhotosEngine::getNumberOfPhotons( const GenVertexPtr theVertex ) const
	-{
	- // Find the number of photons from the outgoing particle list
	-
	- if ( !theVertex ) {
	- return 0;
	- }
	-
	- int nPhotons( 0 );
	-
	- // Get the iterator of outgoing particles for this vertex
	-#ifdef EVTGEN_HEPMC3
	- for ( auto outParticle : theVertex->particles_out() ) {
	-#else
	- HepMC::GenVertex::particles_out_const_iterator outIter;
	- for ( outIter = theVertex->particles_out_const_begin();
	- outIter != theVertex->particles_out_const_end(); ++outIter ) {
	- // Get the next HepMC GenParticle
	- HepMC::GenParticle* outParticle = *outIter;
	-#endif
	-
	- // Get the PDG id
	- int pdgId( 0 );
	- if ( outParticle != nullptr ) {
	- pdgId = outParticle->pdg_id();
	- }
	-
	- // Keep track of how many photons there are
	- if ( pdgId == _gammaPDG ) {
	- nPhotons++;
	- }
	- }
	-
	- return nPhotons;
	-}
	-
	-#endif

File Metadata

Mime Type: text/plain
Expires: Mon, Sep 29, 12:50 PM (19 h, 39 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6569215
Default Alt Text: D112.1759146623.diff (35 KB)

D112.1759146623.diffNo OneTemporaryActions

D112.1759146623.diffView Options

File Metadata

Event Timeline

D112.1759146623.diff
No OneTemporary
Actions

D112.1759146623.diff
View Options