Changeset View
Changeset View
Standalone View
Standalone View
src/EvtGenExternal/EvtPhotosEngine.cpp
Show All 25 Lines | |||||
#include "EvtGenBase/EvtPhotonParticle.hh" | #include "EvtGenBase/EvtPhotonParticle.hh" | ||||
#include "EvtGenBase/EvtRandom.hh" | #include "EvtGenBase/EvtRandom.hh" | ||||
#include "EvtGenBase/EvtReport.hh" | #include "EvtGenBase/EvtReport.hh" | ||||
#include "EvtGenBase/EvtVector4R.hh" | #include "EvtGenBase/EvtVector4R.hh" | ||||
#include <iostream> | #include <iostream> | ||||
#include <sstream> | #include <sstream> | ||||
#include <vector> | #include <vector> | ||||
#include <memory> | |||||
using std::endl; | using std::endl; | ||||
EvtPhotosEngine::EvtPhotosEngine( std::string photonType, bool useEvtGenRandom ) | EvtPhotosEngine::EvtPhotosEngine( std::string photonType, bool useEvtGenRandom ) | ||||
{ | { | ||||
_photonType = photonType; | m_photonType = photonType; | ||||
_gammaId = EvtId( -1, -1 ); | m_gammaId = EvtId( -1, -1 ); | ||||
_gammaPDG = 22; // default photon pdg integer | m_gammaPDG = 22; // default photon pdg integer | ||||
_mPhoton = 0.0; | m_mPhoton = 0.0; | ||||
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Setting up PHOTOS." << endl; | EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Setting up PHOTOS." << endl; | ||||
if ( useEvtGenRandom == true ) { | if ( useEvtGenRandom == true ) { | ||||
EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
<< "Using EvtGen random number engine also for Photos++" << endl; | << "Using EvtGen random number engine also for Photos++" << endl; | ||||
Photospp::Photos::setRandomGenerator( EvtRandom::Flat ); | Photospp::Photos::setRandomGenerator( EvtRandom::Flat ); | ||||
} | } | ||||
Photospp::Photos::initialize(); | Photospp::Photos::initialize(); | ||||
// Increase the maximum possible value of the interference weight | // Increase the maximum possible value of the interference weight | ||||
Photospp::Photos::maxWtInterference( | Photospp::Photos::maxWtInterference( | ||||
64.0 ); // 2^n, where n = number of charges (+,-) | 64.0 ); // 2^n, where n = number of charges (+,-) | ||||
Photospp::Photos::setInterference( true ); | Photospp::Photos::setInterference( true ); | ||||
Photospp::Photos::setExponentiation( true ); // Sets the infrared cutoff at 1e-7 | 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. | // 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, | // 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 | // which is appropriate for B decays | ||||
Photospp::Photos::setInfraredCutOff( 1.0e-7 ); | Photospp::Photos::setInfraredCutOff( 1.0e-7 ); | ||||
_initialised = false; | m_initialised = false; | ||||
} | } | ||||
void EvtPhotosEngine::initialise() | void EvtPhotosEngine::initialise() | ||||
{ | { | ||||
if ( _initialised == false ) { | if ( m_initialised == false ) { | ||||
_gammaId = EvtPDL::getId( _photonType ); | m_gammaId = EvtPDL::getId( m_photonType ); | ||||
if ( _gammaId == EvtId( -1, -1 ) ) { | if ( m_gammaId == EvtId( -1, -1 ) ) { | ||||
EvtGenReport( EVTGEN_INFO, "EvtGen" ) | EvtGenReport( EVTGEN_INFO, "EvtGen" ) | ||||
<< "Error in EvtPhotosEngine. Do not recognise the photon type " | << "Error in EvtPhotosEngine. Do not recognise the photon type " | ||||
<< _photonType << ". Setting this to \"gamma\". " << endl; | << m_photonType << ". Setting this to \"gamma\". " << endl; | ||||
_gammaId = EvtPDL::getId( "gamma" ); | m_gammaId = EvtPDL::getId( "gamma" ); | ||||
} | } | ||||
_gammaPDG = EvtPDL::getStdHep( _gammaId ); | m_gammaPDG = EvtPDL::getStdHep( m_gammaId ); | ||||
_mPhoton = EvtPDL::getMeanMass( _gammaId ); | m_mPhoton = EvtPDL::getMeanMass( m_gammaId ); | ||||
_initialised = true; | m_initialised = true; | ||||
} | } | ||||
} | } | ||||
bool EvtPhotosEngine::doDecay( EvtParticle* theMother ) | bool EvtPhotosEngine::doDecay( EvtParticle* theMother ) | ||||
{ | { | ||||
if ( _initialised == false ) { | if ( m_initialised == false ) { | ||||
this->initialise(); | this->initialise(); | ||||
} | } | ||||
if ( theMother == 0 ) { | if ( theMother == 0 ) { | ||||
return false; | return false; | ||||
} | } | ||||
// Create a dummy HepMC GenEvent containing a single vertex, with the mother | // Create a dummy HepMC GenEvent containing a single vertex, with the mother | ||||
Show All 28 Lines | bool EvtPhotosEngine::doDecay( EvtParticle* theMother ) | ||||
for ( iDaug = 0; iDaug < nDaug; iDaug++ ) { | for ( iDaug = 0; iDaug < nDaug; iDaug++ ) { | ||||
EvtParticle* theDaughter = theMother->getDaug( iDaug ); | EvtParticle* theDaughter = theMother->getDaug( iDaug ); | ||||
GenParticlePtr hepMCDaughter = this->createGenParticle( theDaughter, | GenParticlePtr hepMCDaughter = this->createGenParticle( theDaughter, | ||||
false ); | false ); | ||||
theVertex->add_particle_out( hepMCDaughter ); | theVertex->add_particle_out( hepMCDaughter ); | ||||
if ( theDaughter ) { | if ( theDaughter ) { | ||||
int daugId = theDaughter->getPDGId(); | int daugId = theDaughter->getPDGId(); | ||||
if ( daugId == _gammaPDG ) { | if ( daugId == m_gammaPDG ) { | ||||
nGamma++; | nGamma++; | ||||
} | } | ||||
} | } | ||||
} | } | ||||
// Now pass the event to Photos for processing | // Now pass the event to Photos for processing | ||||
// Create a Photos event object | // Create a Photos event object | ||||
#ifdef EVTGEN_HEPMC3 | #ifdef EVTGEN_HEPMC3 | ||||
▲ Show 20 Lines • Show All 52 Lines • ▼ Show 20 Lines | #endif | ||||
double mass = daugParticle->mass(); | double mass = daugParticle->mass(); | ||||
double energy = sqrt( mass * mass + px * px + py * py + | double energy = sqrt( mass * mass + px * px + py * py + | ||||
pz * pz ); | pz * pz ); | ||||
newP4.set( energy, px, py, pz ); | newP4.set( energy, px, py, pz ); | ||||
// Set the new four-momentum (FSR applied) | // Set the new four-momentum (FSR applied) | ||||
daugParticle->setP4WithFSR( newP4 ); | daugParticle->setP4WithFSR( newP4 ); | ||||
} | } | ||||
} else if ( pdgId == _gammaPDG ) { | } else if ( pdgId == m_gammaPDG ) { | ||||
// Extra photon particle. Setup the four-momentum object | // Extra photon particle. Setup the four-momentum object | ||||
double energy = sqrt( _mPhoton * _mPhoton + px * px + py * py + | double energy = sqrt( m_mPhoton * m_mPhoton + px * px + py * py + | ||||
pz * pz ); | pz * pz ); | ||||
newP4.set( energy, px, py, pz ); | newP4.set( energy, px, py, pz ); | ||||
// Create a new photon particle and add it to the list of daughters | // Create a new photon particle and add it to the list of daughters | ||||
EvtPhotonParticle* gamma = new EvtPhotonParticle(); | EvtPhotonParticle* gamma = new EvtPhotonParticle(); | ||||
gamma->init( _gammaId, newP4 ); | gamma->init( m_gammaId, newP4 ); | ||||
// Set the pre-FSR photon momentum to zero | // Set the pre-FSR photon momentum to zero | ||||
gamma->setFSRP4toZero(); | gamma->setFSRP4toZero(); | ||||
// Let the mother know about this new photon | // Let the mother know about this new photon | ||||
gamma->addDaug( theMother ); | gamma->addDaug( theMother ); | ||||
// Set its particle attribute to specify it is a FSR photon | // Set its particle attribute to specify it is a FSR photon | ||||
gamma->setAttribute( "FSR", 1 ); // it is a FSR photon | gamma->setAttribute( "FSR", 1 ); // it is a FSR photon | ||||
gamma->setAttribute( "ISR", 0 ); // it is not an ISR photon | gamma->setAttribute( "ISR", 0 ); // it is not an ISR photon | ||||
} | } | ||||
▲ Show 20 Lines • Show All 71 Lines • ▼ Show 20 Lines | #endif | ||||
// Get the PDG id | // Get the PDG id | ||||
int pdgId( 0 ); | int pdgId( 0 ); | ||||
if ( outParticle != 0 ) { | if ( outParticle != 0 ) { | ||||
pdgId = outParticle->pdg_id(); | pdgId = outParticle->pdg_id(); | ||||
} | } | ||||
// Keep track of how many photons there are | // Keep track of how many photons there are | ||||
if ( pdgId == _gammaPDG ) { | if ( pdgId == m_gammaPDG ) { | ||||
nPhotons++; | nPhotons++; | ||||
} | } | ||||
} | } | ||||
return nPhotons; | return nPhotons; | ||||
} | } | ||||
#endif | #endif |