diff --git a/EvtGenBase/EvtAbsRadCorr.hh b/EvtGenBase/EvtAbsRadCorr.hh --- a/EvtGenBase/EvtAbsRadCorr.hh +++ b/EvtGenBase/EvtAbsRadCorr.hh @@ -29,6 +29,7 @@ public: EvtAbsRadCorr(){}; virtual ~EvtAbsRadCorr(){}; + virtual void initialise() = 0; virtual void doRadCorr( EvtParticle* p ) = 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 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 #include 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 ) 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 . * -***********************************************************************/ - -#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 - -// 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,6 +36,8 @@ EvtNoRadCorr() { ; } virtual ~EvtNoRadCorr() { ; } + virtual void initialise() { ; } + virtual void doRadCorr( EvtParticle* ) { ; } private: diff --git a/History.md b/History.md --- a/History.md +++ b/History.md @@ -11,6 +11,10 @@ === ## R02-0X-00 +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: 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 +#include +#include -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 ) +{ + 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( 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*/ ) { - 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 . * -***********************************************************************/ - -#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 -#include -#include - -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( 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